1. Introduction

Since their introduction by Bérenger in [1] for the time dependent Maxwell’s equations, Perfectly Matched Layers (PMLs) have become a widely used technique in computational physics. The idea is to enclose the area of interest by surrounded layers which are absorbing and perfectly reflectionless. Subsequent formulations [2, 3] allow to consider PMLs as a complex change of co-ordinates, which implies equivalent material properties [4, 5]. A properly designed change of variables preserves the solution in the region of interest while introducing exponential decay at infinity. It is then obvious to truncate the problem to a finite domain, set a convenient boundary condition on the boundary of this domain and apply the FEM. In a scattering problem, the “traditional co-ordinate stretching” is frequency dependent, yielding a constant attenuation of propagating waves. However, these kinds of PMLs are inefficient for periodic problems when dealing with grazing angles of diffracted orders, i.e. when the frequency is near a Wood’s anomaly [6, 7], leading to spurious reflexions and thus numerical pollution of the results. An important question in designing absorbing layers is thus the choice of their parameters : the PML thickness and the absorption coefficient. To this aim, adaptive formulation have been set up, most of them employing a posteriori error estimate [8, 9, 10]. In this paper, we propose Adaptive PMLs (APMLs) with a suitable co-ordinate stretching, depending both on incidence and grating parameters, capable of efficiently absorbing propagating waves with nearly grazing angles. In the first section, we detail our FEM formulation and expose the problems arising when dealing with Wood’s anomalies. Section 2 is dedicated to the mathematical formulation used to determine PML parameters adapted to any diffraction orders. In the last section, we provide a numerical example of a dielectric slit grating showing the efficiency of our approach in comparison with classical PMLs.

2. Setup of the problem and notations

2.1. Governing equations and description of the structure

We denote by x, y and z the unit vectors of an orthogonal co-ordinate system Oxyz. We deal with time-harmonic fields, so that the electric and magnetic fields are represented by complex vector fields E and H with a time-dependence in exp(−iωt), which will be discarded in the sequel. Moreover, we denote k₀ = ω/c.

The materials in this paper may be z-anisotropic, so the tensor fields of relative permittivity ε͇ and relative permeability μ͇ are of the following form :

The grating is illuminated by an incident plane wave of wave vector defined by the angle θ₀ : k⁺ = α⁺x+β⁺y = k⁺(sinθ₀x − cosθ₀y). Its electric (resp. magnetic) field is linearly polarized along the z-axis, this is the so-called transverse electric or s-polarization case (resp. transverse magnetic or p-polarization case) :

The diffraction problem we are dealing with is to solve Maxwell’s equations in harmonic regime, i.e. to find the unique solution (E, H) of :

Under the aforementioned assumptions, the diffraction problem in non conical mounting can be separated in two fundamental scalar cases TE and TM. Thus we search for a z-linearly polarized electric (resp.magnetic) field E = e(x,y)z (resp. H = h(x,y)z). Denoting ε͇̃ and μ͇̃ the 2 × 2 matrices extracted from ε͇ and μ͇:

The gratings we are dealing with are made of three regions (see Fig. 1) :

 

Fig. 1 Set up of the problem and notations

Download Full Size | PDF

2.2. Implementation of the PMLs

Transformation optics have recently unified various techniques in computational electromagnetics such as the treatment of open problems, helicoidal geometries or the design of invisibility cloaks [4]. These apparently different problems share the same concept of geometrical transformation, leading to equivalent material properties. A very simple and practical rule can be set up [5] : when changing the co-ordinate system, all you have to do is to replace the initial materials properties ε͇ and μ͇ by equivalent material properties ε͇_s and μ͇_s given by the following rule :

In this framework, the most natural way to define PMLs is to consider them as maps on a complex space ℂ³, which co-ordinate change leads to equivalent permittivity and permeability tensors. We detail here the different co-ordinates used in this paper.

In this paper, we use rectangular PMLs [12] absorbing in the y-direction and we choose a diagonal matrix J = diag(1, s_y(y), 1), where s_y(y) is a complex-valued function defined by :

We are now in position to define the so-called substituted field F_s = (E_s, H_s), solution of Eqs. (3) with ξ͇ = ξ͇_s and χ = χ_s. F_s equals the field F in the region y^b < y < y^t (with y^b = −h⁻ and y^t = h^g + h⁺), provided that s_y(y) = 1 in this region (cf. Fig. 2). The complex co-ordinate mapping y(y_c) is simply defined as the derivative of the stretching coefficient s_y(y) with respect to y_c. With simple stretching functions, we can obtain a reliable criterion on the decaying of the fields. A classical choice is :

 

Fig. 2 The basic cell used for the FEM computation of the diffracted field $u_2^d$.

Download Full Size | PDF

In that case, the complex valued function y(y_c) defined by Eq. (9) is explicitly given by :

Let’s now study the behaviour of the field in the PMLs. In the substrate and the superstrate, according to Bloch’s theorem, the diffracted field u^d(x,y) can be written as :

We focus on the bottom PML, as similar considerations can be conducted for the top PML. We consider a diffraction order propagating in the substrate, which can be expressed in the stretched co-ordinate system as :

The non oscillating part of this function is given by :

It then remains to derive the skin depth, $l_n^{-}$, associated with the propagating order n. This characteristic depth is defined as the length for which the field is divided by e :

So $l_n^{-}={\left({\beta }_n^{\prime ,-}{\zeta }^{″,-}+{\beta }_n^{″,-}{\zeta }^{\prime ,-}\right)}^{-1}$, and we take the larger value among the $l_n^{-}$ :

We set the height of the bottom PML region ĥ⁻ = 10l⁻.

2.3. The FEM formulation

Our FEM formulation is that described in [13, 14]. It relies on the fact that the diffraction problem can be rigorously treated as an equivalent radiation problem with sources inside the diffractive object. The expression of the source terms depends on the field u₁ solution of the annex problem of a simple interface, and is known in closed form. In this context, the role of the PMLs is to damp the field radiated by the sources. We set quasi-periodicity conditions on lateral boundaries and Neumann homogeneous conditions are imposed on the outward boundary of each PML. By so doing, the electromagnetic wave is not forced to vanish and the values of the computed field on these boundaries give valuable information on the absorbing efficiency of the PMLs. The cell defined in Fig. 2 is meshed using 2^nd order Lagrange elements [15]. In the numerical examples in the sequel, the maximum element size is set to ${\lambda }_0/\left(N_m\sqrt{\Re 𝔢\left({\epsilon }_r\right)}\right)$, where N_m is an integer (between 6 and 10 is usually a good choice). The final algebraic system is solved using a direct solver (PARDISO).

Note that we use this FEM formulation in this paper but the results of the sequel about PMLs are general and can be applied to any numerical method used for the modelling of diffraction gratings.

2.4. Weakness of the classical PML for grazing diffracted angles

Resuming part 2.2, we consider only the bottom PML, as analogous conclusions holds for the top PML. The efficiency of the classical PML fails for grazing diffracted angles, in other words when a given order appears/vanishes : this is the so-called Wood’s anomaly, well known in the theory of gratings. In mathematical words, there exists n₀ such that ${\beta }_{n_0}^{-}\simeq 0$. The skin depth of the PML then becomes very large. To compensate it, one should increase the value of ζ^″,−, but this gives rise to spurious numerical reflections due to an overdamping. For a fixed value of ĥ⁻, if ζ^″,− is too weak, the absorption in the PMLs is insufficient and the wave reflects on the outward boundary of the PML. To illustrate these typical behaviours (cf. Fig. 3), we compute the diffracted field by a grating with a rectangular cross section of height h^g = 1.5μm and width L^g = 3μm with ε^g = 11.7, deposited on a substrate with permittivity ε⁻ = 2.25. The structure is illuminated by a p-polarized plane wave of wavelength λ₀ = 10μm and of angle of incidence θ₀ = 10° in the air (ε⁺ = 1). All materials are non magnetic (μ_r = 1) and the periodicity of the grating is d = 4μm. We set ${\widehat{h}}^{-}=10l_0^{-}$ and ζ^′,− = 1.

 

Fig. 3 Zero^th transmitted order by a grating with a rectangular cross section (see parameters in text, part 2.4) for different values of ζ^″,− : blue line, ζ^″,− = 1, correct damping; green line, ζ^″,− = 0.1, underdamping; red line, ζ^″,− = 20, overdamping.

Download Full Size | PDF

3. Construction of an adaptive PML

To overcome the problems pointed out in the preceding section, we propose a co-ordinate stretching that rigorously treats the problem of Wood’s anomalies. The wavelengths “seen” by the system are very different depending on the order at stake :

Thus if a classical PML is adapted to one diffracted order, it will not be for another, and vice versa. The idea behind the APML is to deal with each and every order when progressing in the absorbing medium.

Once again the development will be conducted only for the PML adapted to the substrate. We consider a real-valued co-ordinate mapping y_d(y), the final complex-valued mapping is then y_c(y) = ζ⁻y_d(y), with the complex constant ζ⁻, with ζ^′,− > 0 and ζ^″,− > 0, accounting for the damping of the PML medium. We start by transforming Eq. (15), $n\mapsto {\beta }_n^{-}$ a function with integer argument into a function with real argument continuously interpolated between the imposed integer values. Indeed, the geometric transformations associated to the PML has to be continuous and differentiable in order to compute the Jacobian. For this purpose, we choose the parametrization :

The function y(y_d) has two poles, denoted $y_{d,\pm }^{★}=d\left(\pm \sqrt{{\epsilon }^{-}}-\text{sin}{\theta }_0\right)$. When $y_{d,\pm }^{★}=\pm n{\lambda }_0$ with n ∈ ℕ^★, ${\beta }^{-}\left(y_{d,\pm }^{★}\right)={\beta }^{-}\left(\pm n{\lambda }_0\right)={\beta }_{\pm }^{-}=0$, i.e. we are on a Wood’s anomaly associated with the appearance/disappearance of the ±n^th transmitted order. We now search for the nearest point to $y_{d,\pm }^{*}$ associated with a Wood’s anomaly, denoting :

To avoid the singular behaviour at $y_d=y_{d,\pm }^{★}$, we continue the graph of the function y_d(y) by a straight line tangent at $y_d^0$, which equation is $t_0\left(y_d\right)=s\left(y_d^0\right)\left(y_d-y_d^0\right)+y\left(y_d^0\right)$, where $s\left(y_d\right)=\frac{\partial y}{\partial y_d}\left(y_d\right)$ is the so-called stretching coefficient. The final change of co-ordinate is then given by :

 

Fig. 4 Example of a co-ordinate mapping ỹ(y_d) used for the APML (black solid line). The graph of y_d(y) (blue solid line) is continued by a straight line t₀(y_d) tangent at $y_d^0$ (red dashed line) to avoid the singular behaviour at $y_d=y_d^{★}$.

Download Full Size | PDF

Eventually, the complex stretch s_y used in Eq. (10) is given by :

Equipped with this mathematical formulation, we can tailor a layer that is doubly perfectly matched :

4. Numerical example

We now apply the method described in the preceding parts to design an adapted bottom PML for the same example as in part 2.4. The parameters are the same, and we choose the wavelength of the incident plane wave close to the Wood’s anomaly related to the +1 transmitted order ( ${\lambda }_0=0.999y_{d,+}^{★}$). Moreover, we set the length of the PML ${\widehat{h}}^{-}=1.1y_{d,+}^{★}$ and choose absorption coefficients ζ⁺ = ζ⁻ = 1 + i. For both cases (PML and APML), parameters are alike, the only difference being the complex stretch s_y.

The field maps of the norm of H_z, E_x and E_y are plotted in logarithmic scale on Fig. 5, for the case of a classical PML and our APML. We can observe that the field H_z that is effectively computed is clearly damped in the bottom APML (leftmost on Fig. 5(b)) whereas it is not in the standard case (leftmost on Fig. 5(a)), causing spurious reflections on the outer boundary. The fields E_x and E_y are deduced from H_z thanks to Maxwell’s equations. The high values of E_y at the tip of the APML (rightmost on Fig. 5(b)) are due to very high values of the optical equivalent properties of the APML medium (due to high values of s_y), which does not affect the accuracy of the computed field within the domain of interest.

 

Fig. 5 Field maps of the logarithm of the norm of H_z, E_x and E_y for the dielectric slit grating at ${\lambda }_0=0.999y_{d,+}^{★}$ (same parameters as in part 2.4). (a) : classical PML with inefficient damping of H_z in the bottom PML. (b) : APML where the H_z field is correctly damped in the bottom PML. For both cases the thickness of the PML is ${\widehat{h}}^{-}=1.1y_{d,+}^{★}$.

Download Full Size | PDF

Another feature of our approach is that it efficiently absorbs the grazing diffraction order, as illustrated on Fig. 6 : the +1 transmitted order does not decrease in the standard PML (blue solid line), and reaches a high value at y = −ĥ⁻, whereas the same order tends to zero as y → −ĥ⁻ in the case of the adapted PML (blue dashed line).

 

Fig. 6 Modulus of the u_n for the three propagating orders with adapted (dashed lines) and classical PMLs (solid lines). Note that the classical PMLs are efficient for all orders except for the grazing one (n = 1) as expected. This drawback is bypassed when using the adaptive PML.

Download Full Size | PDF

To further validate the accuracy of the method, we compare the diffraction efficiencies computed by our FEM formulation with PML and APML to those obtained by another method. We choose the Rigorous Coupled Wave Analysis (RCWA), also known as the Fourier Modal Method (FMM, [16]). For the chosen parameters, only the 0^th order is propagative in reflexion and the orders −1, 0 and +1 are non evanescent in transmission. We can also check the energy balance B = R₀ + T₋₁ + T₀ + T₊₁ since there is no lossy medium in our example. Results are reported in Table 1, and show a good agreement of the FEM with APML with the results from RCWA. On the contrary, if classical PML are used, the diffraction efficiencies are less accurate compared to those computed with RCWA. Checking the energy balance leads the same conclusions : the numerical result is perturbed by the reflection of the waves at the end of the PML if it is not adapted to the situation of nearly grazing diffracted orders.

Table 1. Diffraction efficiencies R₀, T₋₁, T₀ and T₊₁ of the four propagating orders, and energy balance B = R₀ + T₋₁ + T₀ + T₊₁, computed by three methods : RCWA (line 1), FEM formulation with APML (line 2), FEM formulation with classical PML (line 3).

View Table

Eventually, to illustrate the behaviour of the adaptive PML when the incident wavelength gets closer to a given Wood’s anomaly, we computed the mean value of the norm of H_z along the outer boundary of the bottom PML γ = 〈|H_z(−ĥ⁻)|〉_x, when ${\lambda }_0=\left(1+{10}^{-n}\right)y_{d,+}^{★}$ and ${\lambda }_0=\left(1-{10}^{-n}\right)y_{d,+}^{★}$, for n = 1, 2,...10. The results are shown in Fig. 7. As the wavelength gets closer to $y_{d,+}^{★}$, γ first increases but for n > 3, it decreases exponentially. However, in all cases, the value of γ remains small enough to ensure the efficiency of the PMLs.

 

Fig. 7 Mean value of the norm of H_z along the outer boundary of the bottom PML γ = 〈|H_z(−ĥ⁻)|〉_x, for λ₀ approaching the Wood’s anomaly $y_{d,+}^{★}$ by inferior values ( ${\lambda }_0=\left(1-{10}^{-n}\right)y_{d,+}^{★}$, red squares) and by superior value ( ${\lambda }_0=\left(1+{10}^{-n}\right)y_{d,+}^{★}$, blue circles) as a function of n.

Download Full Size | PDF

5. Conclusion

In this paper, we have proposed an adaptive PML that can treat rigorously Wood’s anomalies in numerical analysis of diffraction gratings. It is based on a complex-valued co-ordinate stretching that deals with grazing diffracted orders, yielding an efficient absorption of the field inside the PML. We provided an example in the TM polarization case (but similar results hold for the TE case), illustrating the efficiency of our method. The value of the magnetic field on the outward boundary of the PML remains small enough to consider there is no spurious reflection. The formulation is used with the FEM but can be applied to others numerical methods. Moreover, the generalization to the vectorial three-dimensional case is straightforward : the recipes given in part 3 do work irrespective of the dimension and whether the problem is vectorial. In addition, although the designed co-ordinate stretching is specific to the context of Wood’s anomalies, one can adapt this kind of non uniform PML to others wave problems by following the methodology exposed here.