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 exact 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₃ [20], LiNbO₃ [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 }$ _a (resp. µ̄_a) represents the conjugate complex of ε_a (resp. µ_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 ₀:=ω/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^{+}≔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^{-}≔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.

 

Fig. 1. Schematics of the gratings studied in this paper.

Download Full Size | PDF

This grating is illuminated by an incident plane wave of wave vector k ⁺=αx+β+y=k ⁺ (sinθ0x-cosθ₀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 ⁰ _e (resp. A ⁰ _m) is an arbitrary complex number. The magnetic (resp. electric) field derived from E ⁰ _e(resp. H ⁰ _m) is denoted H ⁰ _e(resp. E ⁰ _m) and the electromagnetic field associated with the incident field is therefore denoted (E ⁰,H ⁰) which is equal to (E ⁰ _e,H 0 _e) (resp. (E ⁰ _m,H ⁰ _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=δ_zy=0). Moreover, the left upper matrix extracted from δ͇ is denoted $\underset{¯}{\underset{¯}{\widehat{\delta }}}$ , namely:

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 ₁ and u ^d₂, 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₂ 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 u0 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 ₁, is defined as the unique solution of the equation:

The function u ₁ 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₂ which is simply defined as the difference between u and u ₁:

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₂ 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 𝒮₁ is thus localized within the groove region as expected. It remains to compute more explicitly the source term 𝒮₁. Making use of the linearity of the operator 𝓛 and the equality u ₁=u ^d ₁+u0, the source term can be splitted into two terms

where

and

Now, bearing in mind that u0 is nothing but a plane wave u ₀=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):

The same holds for the term associated with the diffracted field (u^d ₁=ρ 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 ²(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 ²(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 Γ_r the lines parallel to the 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 ²(curl,d,k), the boundary term for Γ_l∪Γ_r is

${\int }_{{\Gamma }_l\cup {\Gamma }_r}\overline{u\prime }(\underset{¯}{\underset{¯}{\xi }}\nabla u)·\mathbf{n}dS={\int }_{{\Gamma }_l\cup {\Gamma }_r}\overline{u_{\#}^{\text{'}}}{e}^{-\mathrm{ikx}}(\underset{¯}{\underset{¯}{\xi }}\nabla \left(u_{\#}{e}^{+\mathrm{ikx}}\right))·\mathbf{n}dS=$${\int }_{{\Gamma }_l\cup {\Gamma }_r}\overline{u_{\#}^{\text{'}}}\left(\underset{¯}{\underset{¯}{\xi }}\left(\nabla u_{\#}+\mathrm{ik}{u}_{\#}\mathbf{x}\right)\right)·\mathbf{n}dS=0$

because the integrand $\overline{u_{\#}^{\text{'}}}\left(\underset{¯}{\underset{¯}{\xi }}\left(\nabla u_{\#}+\mathit{ik}{u}_{\#}\mathbf{x}\right)\right)·\mathbf{n}$ is periodic along x and the normal n has opposite directions on Γ_l and Γ_r 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.

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 Γ_l (i.e. shifted by a quantity -d along x) up to the factor e^ikd. The practical implementation in the finite element method is described in details in [25, 26]

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{¯}{\underset{¯}{{\epsilon }_s}}(x,y)$ and $\underset{¯}{\underset{¯}{{\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_s=diag (1, 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{¯}{\underset{¯}{{\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)^T which is solution of Eq. (3) except that we have replaced ε͇ and µ by $\underset{¯}{\underset{¯}{{\epsilon }_s}}$ and $\underset{¯}{\underset{¯}{{\mu }_s}}$. The main feature of this latest field is the remarkable correspondence with the first field F ; whatever the function sy provided that it equals 1 for Y+∗<y<Y-∗ (cf. Fig. 2a), the two fields F and F _s are identical in the region Y+∗<y<Y-∗ [27]. In other words, the PML is completely reflectionless. In addition, for complex valued functions s_y (∑m{s_y} strictly positive in PML), the field F_s does converge exponentially towards zero (as y tends to ±∞, cf. Fig. 2c and 2d) although its physical counterpart F does not. As a consequence, F_s 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.

 

Fig. 2. PML adapted to substratum and superstratum (TM case). Note (Fig. 2d) that the radiated field on each extreme boundary of the PML is at least 10^-8 weaker than in the region of interest.

Download Full Size | PDF

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}{y}_c>Y_{*}^{+}\\ y_c& \mathrm{for}{Y}_{*}^{-}<y_c<Y_{*}^{+}\\ Y_{*}^{-}+{\zeta }^{-}\left(y_c-Y_{*}^{-}\right)& \mathrm{for}{y}_c<Y_{*}^{-}.\end{array}$.

Now, let us consider a propagating plane wave in the substratum u_n(x,y) :=exp(i(αx-β-+,_ny)) can be reexpressed in stretched coordinates in the PML 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^sc(yc) is given by exp ((β′-+,_nζ″-+β″-+,_nζ′-)yc). Keeping in mind that β′-+,_n and/or β″-+,_n are positive numbers, the function U^sc decreases exponentially towards zero as y_c tends to -∞ (Fig. 2d) provided that ζ- belongs to ℂ+ :={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₂ 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 χ+_s depending on the equivalent anisotropic properties of the PML given by Eq. (7), Eq. (8), Eq. (27) and Eq. (28).

• 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 χ^g 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).

• 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_n=0. Combining Eq. (36) and Eq. (37) at a fixed y ₀ 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 ₀. 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 ofC∞ class. From a numerical point of view, it appears that the interpolated approximation of the unknown function, namely u^d₂ 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_n and R_n) defined by :

This calculation is performed at twenty different y ₀ 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.9524i. We present the R ₀ 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.

 

Fig. 3. Rectangular groove grating : This pattern is repeatedly set up with a period d=1 µm. This grating has been studied by [23] and is one of our points of reference.

Download Full Size | PDF

Table 1. Reflected efficiencies versus mesh refinement. Note that the efficiencies are properly computed (two significant digits) even for a rather coarse mesh.

View Table | View all tables in this article

4.3. Experiment set up based on existing materials

Let us now consider a trapezoidal (cf. Fig. 4) anisotropic grating made of aragonite (CaCO₃) deposited on an isotropic substratum (SiO₂, εS_iO₂=2.25). Along the anisotropic crystal axis, its dielectric tensor can be written as follows [20] :

 

Fig. 4. Diffractive element pattern. This element is made of aragonite for which the dielectric tensor is given by Eq.(41) and is deposited on a silica substrate (d=600nm).

Download Full Size | PDF

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 λ ₀=633nm 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 λ ₀/10√ε. Efficiencies are still post-processed in accordance with the calculation presented section 4.1.

4.4. Numerical results

 

Fig. 5. Real part of the total calculated field depending on θ₀ (TM and TE cases).

Download Full Size | PDF

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 ^TE1 (cf. Table 2: 0.322510 versus 0.124722), whereas T ^TM-1=T ^TM1=0.20313. Note that the material asymmetry ratio $\frac{0.251}{2.592}$ between the diagonal term and off-diagonal term of ${\underset{¯}{\underset{¯}{\epsilon }}}_{\mathrm{Ca}{\mathrm{CO}}_3}^{+45°}$ is much smaller than the optical asymmetry ratio at normal incidence, i.e. $\frac{T_1-T_{-1}}{T_1}$.

Table 2. Transmitted and reflected efficiencies of propagative orders deduced from field maps shown Fig.5

View Table | View all tables in this article

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.