## Abstract

An extension of the dual weighted residual (DWR) method to the analysis of electromagnetic waves in a periodic diffraction grating is presented. Using the *α*,0-quasi-periodic transformation, an upper bound for the *a posteriori* error estimate is derived. This is then used to solve adaptively the associated Helmholtz problem. The goal is to achieve an acceptable accuracy in the computed diffraction efficiency while keeping the computational mesh relatively coarse. Numerical results are presented to illustrate the advantage of using DWR over the global *a posteriori* error estimate approach. The application of the method in biomimetic, to address the complex diffraction geometry of the Morpho butterfly wing is also discussed.

## 1. Introduction

This paper investigates the interaction of electromagnetic waves with periodic diffraction gratings. This type of grating has been used recently in crystalline silicon solar cells [1], gas sensors [2] and medical X-ray imaging [3,4]. They are also used in acousto-optic devices to characterize the composition of drugs, and in spectroscopic sensing for detecting trace gases [5]. Such periodic structures also occur naturally in the Morpho butterfly wing. The iridescent scales of this butterfly display a different optical response to the presence of different vapours [2]. The information about the nature and concentration of the vapours provided by the reflectance spectra of the scales has therefore been proposed as a design for a gas sensor.

Theoretically, electromagnetic wave diffraction is based on solving Maxwell’s equations in the diffraction grating region [6,7]. If it is assumed that there is neither charge nor current, then the Maxwell equations reduce to the Helmholtz equation [7]. This boundary value problem has been investigated in two dimensions, for a periodic grating, using a finite-element method [8]. This work also introduced the concept of the *α*_{0}, 0-quasi-periodic transformation which is used in this paper. In order to achieve a reliable and efficient numerical method, it is necessary to know how the system parameters affect the error between the exact and approximate solution. These errors can be classified into three different sources; data errors, model errors and computation errors [9]. As the finite-element method is used here, then the focus will be on the error arising from the discretization (data and modelling). For more details on the finite-element modelling in periodic structure, see [10–12].

An *a priori* error estimate has already been derived in Lord & Mulholland [8] to guarantee that the discrete solution converges to the exact solution provided that the interpolation error is kept small with respect to the wavenumber *k*. This study found an upper bound on the error between the exact and the approximate solutions in terms of the exact solution and the solution to the dual problem. Hence, this upper bound depends on some stability constants whose dependency on the system parameters is unknown. This limits the numerical approach to using a uniform mesh, which is computationally expensive. In this paper, an *a posteriori* error estimate that arises when we discretize the Helmholtz problem for a periodic grating is derived. In contrast to the *a priori* error estimate, it will provide a computable upper bound so that local refinement of the mesh (h-version) or raising the degree of the polynomial basis (p-version) can be performed. By doing so, this will reduce and control both the numerical error and the computational cost. This process is known as the adaptive computation method [9,13,14]. There are two different approaches to estimate the approximation error (denoted *e*_{h}). The first approach consists of looking for upper and lower bounds on *e*_{h} (global error estimate) and the second one is to use a quantitative estimate of some local feature of the exact solution called the quantity of interest (strain, displacement, etc.) [9], and to look for an estimate of the error in this quantity of interest (goal-oriented error estimate). In many applications, the interest lies in the error that arises in some specific, real valued physical quantity of interest *Q* that depends on the solution *U*_{α,0}. The global error does not, however, provide useful bounds for this error in the quantity of interest. Also the sensitivity of the global error to local error sources is not properly represented when global stability constants are used [9,14]. This paper will employ a goal-oriented error estimate so that the error in the quantity of physical interest can be controlled and at the same time the efficiency of computing this quantity can be optimized. In particular, this paper focuses on the dual weighted residual (DWR) method [13]. Here, the dual problem is combined with the direct problem to derive an estimate of the error in the target quantity from each local residual on each of the mesh cells. In this way, the error is controlled locally in computing the target quantity. Given a grating profile, one of the main concerns is to know how much of the incident wave will be reflected and how much will be transmitted. This is achieved by computing the efficiency of the grating. Hence, rather than estimating the energy norm of the error in the solution to the Helmholtz problem, a particular linear functional of this error is estimated. This linear functional denoted by *Q* is chosen so that the error made by computing the efficiency of the grating is minimized.

There are two types of grating, the perfectly conducting and the dielectric transmitting grating [8]. There are also two fundamental polarizations, the transverse electric (TE) mode and the transverse magnetic (TM) mode [7,15,16]. Hence, there are four cases to investigate: Case 1A/B: perfectly conducting; and Case 2A/B: transmitting dielectric where A (B) denotes the TE (TM) wave. In §2, the mathematical statement of the problem is posed both in its continuous and discrete forms. In §3, the dual problem is introduced, which then facilitates the derivation of the *a posteriori* error estimate in §4. The final section provides some numerical results and applications.

## 2. Problem statement

In this section, the geometry of the problem, the associated function spaces, the form of the Helmholtz problem under the *α*_{0},0-quasi-periodic transformation and the boundary conditions are detailed.

### (a) Boundary value problem

Using the *α*_{0},0-quasi-periodic transformation with the finite-element method, it has been shown in Lord & Mulholland [8] that the periodic diffraction grating problem can be solved numerically inside the truncated single vertical strip, *Ω*=[0,*d*]×[−*B*,*B*] as shown in figure 1. Define *Γ*_{+}={(*x*,*y*):0≤*x*≤*d*,*y*=*B*} and *Γ*_{−}={(*x*,*y*): 0≤*x*≤*d*,*y*=−*B*}. Also define the wavenumber *k* to be
2.1where *Ω*_{1}=[0,*d*]×[*b*,*B*], *Ω*_{0}=[0,*d*]×[−*b*,*b*] and *Ω*_{2}=[0,*d*]×[−*B*,−*b*] with *d* the grating period, *b* and *B* positive real numbers. The incident wave is given by , where , and *θ* is the angle of incidence of the wave as shown in figure 1. The scattered and diffracted waves are composed of bounded outgoing waves by demanding that >0 and ≥0, where () denotes the real (imaginary) part of *k*_{j} for *j*∈{0,1,2}.

The notation below will allow us to define the various function spaces arising on the boundaries of the domains in figure 1 which are
and the function spaces used inside *Ω* which are
with , *H*^{s}(*Ω*) and *H*^{s}([0,*d*]) are Sobolev spaces [17,18].

### Definition 2.1

Let . We define the Dirichlet to Neumann (DtN) maps, , where , and
2.2for *n*∈*Z*, *α*_{n}=*α*_{0}+2*πn*/*d* and , such that , for *j*∈{1,2} with . The case where corresponds to the resonance phenomenon and is not considered here [7]. It can be shown that is a continuous function from [11, p. 29].

The solution to the scattering problem is *α*_{0}-quasi-periodic, that is it can be written as the product of a periodic function with e^{iα0x}. By denoting this periodic function by *U*_{α,0}, the *α*_{0}, 0-quasi-periodic transformation can be used [8]. Hence, the scattering problem when an incident wave interacts with a periodic grating can be reformulated as *AU*_{α,0}=*f*, where
2.3for Cases 1A, 1B and 2A and
2.4for Case 2B. ∇_{α,0} is given by ∇+i(*α*_{0},0), *I*_{d} is the identity operator, *f* is some given data, is a linear operator on , and
2.5For all cases, *U*_{α,0} also satisfies the following boundary conditions:
2.6and
2.7where *n* represent the outward unit normal and in addition for Case 1
2.8

The exact solution of the boundary value problem is denoted by *U*_{α,0} and *U*_{α,0h}∈*X*^{α}, such that is the corresponding numerical solution.

### (b) Continuous problem

The continuous variational formulations given for the four cases are represented as a single problem as follows. Find such that
2.9for all with
2.10and
2.11where
2.12such that *k* is given by equation (2.1) and *n*_{x} is the outward unit normal with respect to the *x*-axis. For Case 2B, because both *U*_{α,0} and *v* are periodic, then
2.13

### (c) Discretized problem

If *X*^{α} is a finite-dimensional subspace of [17–19], then the discrete problem is to find *U*_{α,0h}∈*X*^{α} such that
2.14for all *v*_{h}∈*X*^{α}, with *a*(⋅,⋅) and (*f*,⋅) as given, respectively, by equations (2.10) and (2.11).

Where *q* and *F* are given, respectively, by equations (2.12) and (2.5). The discretization error is denoted by
2.15In order to establish an *a posteriori* error estimate, the dual problem is studied in §3.

## 3. Dual problem

Here, the DWR method is used to provide a goal-oriented error estimate [13]. A quantity of interest is defined and shown to be a continuous linear functional, as this is a necessary condition for the dual problem to be well posed. The dual problem is then formulated and used to study the error in the quantity of interest.

### (a) Quantity of interest

Let , and define the map *Q*, where , as
3.1such that *c*_{n}=0 or *c*_{n}=1 and the Fourier coefficient *f*^{(n)}(±*B*) of *f*(*x*,±*B*) is given by
Relating these to the diffraction efficiency [8, §§ 4.1.4 and 4.2.3] gives
because , for |*α*_{n}|<|*k*_{2}|, where *c*_{n}=1 was used in equation (3.1) to keep notation simple. and are complex scalars called Rayleigh coefficients used to compute the diffraction efficiency of order *n* [8]. Hence, the quantity of interest is chosen so that it is related to the computation of the diffraction efficiencies corresponding to the propagating waves, given the condition |*α*_{n}|<|*k*|. In addition, the constants *c*_{n} have been introduced so that a particular efficiency order *m* can be isolated. Hence, *c*_{n}=1 if *m*=*n* and *c*_{n}=0 otherwise.

### Lemma 3.1

*Let* , *k*_{ref}>0 *such that* |*k*|>*k*_{ref}, *C a positive constant and let Q be given by equation* (3.1), *then Q is a continuous linear functional and*

### Proof.

The linearity of *Q* follows from its definition. From Schwarz’s inequality
As the sum is finite and bounded by a constant. We can estimate this constant as follows:
with *C* a positive constant, then the change of variable *x*=*α*_{n}−*α*_{0} gives
Evaluating the integrals gives
Hence,
from the definition of the -norm [20, p. 150] for . □

### (b) Strong form and variational form of the dual problem

A dual function is introduced in order to estimate *Q*(*U*_{α,0})−*Q*(*U*_{α,0h})=*Q*(*e*_{h}), and it satisfies
3.2where *a* is given by equation (2.10) and *A* is given by equation (2.3) or (2.4). Hence, the strong form of the dual problem given by equation (3.2) is provided by lemma 3.2.

### Lemma 3.2

*Let* *be the dual solution corresponding to the dual problem of equation* (2.10) *then z satisfies*
3.3*The functional J*_{±} *is defined as*
3.4*with j*=1,2, *respectively, for J*_{+} (*on Γ*_{+}) *and J*_{−} (*on Γ*_{−}), *and* *is the dual of* [21, p. 476].

### Proof.

The divergence theorem and applying integration by parts to *a*(*v*,*z*) in equation (3.2) gives
and
where *n*_{i} is the outward unit normal in the direction of *x*_{i}. This leads to
from equation (2.13). Hence, for Cases 1A and 2,
for all . Let , then
and . Note that
3.5such that *J*_{±} is given by (3.4). For Case 1B, (*vn*_{x},*z*)_{∂F}=(*v*,*n*_{x}*z*)_{∂F}. □

As the exact dual solution *z* is not known, a weak formulation is used in the numerical implementation to approximate *z*. From lemma 3.2, the variational formulation corresponding to the dual problem is given as follows.

Let (weak solution) then, for any , for Cases 1A, 2A and 2B 3.6and for Case 1B Having formulated the direct and the dual problem, they can be used to establish a goal-oriented error estimate, using the DWR method, where the quantity of interest is related to the grating efficiency.

## 4. *Aposteriori* error estimation

Let *z* be the dual solution associated with the dual problem given by lemma 3.2. Then, for any ,
such that
where *A** is the dual operator of *A* defined in equation (2.3) or (2.4), [21, p. 476]. Using equation (3.5) leads to the following problem: find such that
4.1for any

Let *X*^{α} be a finite-element subspace of order *p* of and let *ζ*_{h} be any regular partition of *X*^{α} [17,18,22,23]. Denote by *h* the maximum mesh size of the triangular elements in this partition. The finite-element approximation of the dual problem given by equation (4.1) is to find *z*_{h}∈*X*^{α} such that , for any *ϕ*_{h}∈*X*^{α}.

The estimate of the linear functional of the error *Q*(*e*_{h}) is given by the following theorem.

### Theorem 4.1

*Let U*_{α,0h} *be the solution to equation (*2.14*), e*_{h} *be given by equation (*2.15*), ζ*_{h}*={K} be a partition of X*^{α}*, and p*_{K} *and h*_{K} *be the polynomial order and the mesh size associated with the element K. Denote the field equation residual by R*_{h}*(K), then*
4.2*The flux residual r*_{h}*(E) for Case 2 is given by*
4.3*and for Case 1 by*
4.4*with [∂*_{n}*U*_{α,0h}] *denoting the jump of the normal derivatives of U*_{α,0h}, *E an edge of the element K and q as defined by equation* (2.12). *Then*
4.5*where the cell residuals ρ*_{k} *and weights w*_{K} *are given by*
4.6*and*
4.7*with* .

### Proof.

From equations (3.2) and (2.14), *Q*(*U*_{α,0})−*Q*(*U*_{α,0h})=*a*(*U*_{α,0}−*U*_{α,0h},*z*), and using equations (2.9) and (3.2) *Q*(*e*_{h})=(*f*,*z*)_{Γ+}−*a*(*U*_{α,0h},*z*)=(*f*,*z*)_{Γ+}−(*AU*_{α,0h},*z*). Let *ζ*_{h} be a partition of the domain *F* into mesh cells *K*. From cellwise integration by parts [24, p. 28], [25, p. 12] and by using Galerkin orthogonality [17, p. 58] with equation (2.3) for Cases 1A, 1B and 2A give
for all *ϕ*∈*X*^{α} and . Using a similar argument for Case 2B gives
Using equations (2.6)–(2.8) for Case 1 gives equations (4.2)–(4.4). Then
From the Cauchy–Schwarz inequality [17, p. 50]
The proof is completed by choosing *ϕ*=*z*_{h} and by using equations (4.6) and (4.7). □

## 5. Numerical methods and numerical results

For completeness a pseudo-code algorithm for the implementation of the method proposed in this paper is detailed in §5*a*. After that, three numerical examples are used to illustrate the computational benefits of this adaptive method.

### (a) Algorithm for the adaptive mesh

The operations are done elementwise by looping over all elements of a given triangulation. Hence, each element is mapped to a reference element through an affine transformation.

— step 0. Choose a tolerance TOL;

— step 1. Assemble the mass matrix for where

*ψ*_{i}is some basis in the corresponding finite-element subspace*X*^{α}.Using a reference element,

*ψ*_{i}becomes . Write*a*(*ψ*_{i},*ψ*_{j})=*K*_{ij}+*B*_{ij}and use equation (2.10) to give Numerical integration is used to calculate the surface and element integrals to get*K*_{ij}[26–28];— step 2. Assemble the load vector.

Similar to step 1, equation (2.11) becomes which is also computed using numerical integration;

— step 3. Assemble the DtN operators.

From step 1, , where 5.1 can be computed as follows using equation (2.2): 5.2where

*x*_{j}and*x*_{j+1}delimit the element edge where is different from zero and hence . Since is a polynomial of order*N*then Hence, equation (5.2) becomes 5.3where It can be shown by integrating by parts and by induction that 5.4for*q*≤*l*−1. Equations (5.1), (5.3) and (5.4) can then be used to compute*B*_{ij};— step 4. Apply the periodicity boundary conditions for all cases and the Dirichlet boundary conditions for Case 1A.

The problem is find

*U*_{α,0}such that . By denoting and by noting that the solution satisfies 5.5for

*i*,*j*=1,…*N*. For the periodic boundary condition, let*N*_{p}denote the set of nodes*i*such that*ϕ*_{i}belongs to the boundary*x*=0 or*x*=*d*and*N*_{O}denote the set of nodes*j*such that*ϕ*_{j}belongs to ∂*Ω*_{3}for Case 1A. The techniques described in Ainsworth [29] are used to apply the periodicity constraints and the Dirichlet constraints to the nodes*i*∈*N*_{p}and*j*∈*N*_{O};— step 5. Solve the system.

Solve system (5.5) and find the coefficients of the approximate solution to the Helmholtz problem;

— step 6. Repeat steps 1–5 to solve the dual problem to find

*z*_{h}using equation (3.6).As the variational form of the dual problem as given in equation (3.6) has a form similar to the direct problem, the same technique is used to solve the dual problem;

— step 7. Repeat steps 1–5 to solve the dual problem using equation (3.6) with a higher polynomial order or a finer mesh to give an approximation.

— step 8. Compute the upper bound given by equation (4.5) and denoted this upper bound by 5.6Standard techniques can be used (reference element and numerical integration) [26,28,30,31] to compute equation (4.2). For the flux residual given by equations (4.3) and (4.4), the are computed using definition 2.1 and equation (5.3). For the jump derivative [∇

*ψ*_{i}.*n*], let us denote by*K*and*K*^{1}two triangles which share an edge*E*, and let the affine transformation which maps*K*to the reference element (and*ψ*_{i}to ), be described by where*J*is a 2×2 matrix. Let the affine transformation which maps*K*^{1}to the reference element (and to ), be described by where*J*^{1}is a 2×2 matrix. We then have because*n*^{1}=−*n*. Hence, using*y*=*g*(*x*) to represent the curve , gives where the edge*E*is delimited by*g*(*x*_{0}) and*x*_{0}.— step 9. Let

*N*_{T}be the list of elements*K*where*ρ*(*K*)*w*_{K}is in decreasing order and such that . Refine the mesh in*N*_{T};— step 10. Check that the mesh is periodic if it is not defined

*Γ*_{L}={*y*:(0,*y*)∈*Ω*} and*Γ*_{R}={*y*:(*d*,*y*)∈*Ω*}. Refine the edge elements on these boundaries until*Γ*_{L}=*Γ*_{R}so that the mesh is periodic; and— repeat steps 1–10 until

*I*_{err}≤TOL.

### (b) Numerical results

Having briefly described the numerical technique used to solve numerically the Helmholtz problem, the benefits that can be obtained by using an adaptive grid, and in particular where this adaptivity is driven by the DWR method is illustrated below.

### Example 5.1

In the following example, the transmitting dielectric lamellar grating as shown in figure 2 is considered. This type of grating, used in modelling multiscale phenomena grating problems, has been studied in Bao *et al.* [32] using a hybrid approach that combines a perfectly matching layer technique and an adaptive finite-element method driven by a global *a posteriori* error estimate and it can be applied in solving optimal design problems. Here, the wavenumbers are fixed as *k*_{1}=2*π* and *k*_{2}=(0.22+6.71i)2*π*, the angle of incidence is *θ*=*π*/6 and the period is *d*=1. A TM field (Case 2B) is considered as described by equations (2.10) and (2.11), and the reflection efficiency of order zero is computed [8], sections 4.2.3. To use the DWR method, equations (3.6) and (4.5) are employed to find the dual solution and the error bounds. The algorithm in §5*a* is followed with TOL=10^{−4}. Let d.f. denote the degrees of freedom, that is the total number of nodes needed to define all triangular elements on the domain, and all the polynomial bases on each of these elements. To provide a basis for a relative error, the global method in Bao *et al.* [32] is used with 201205 d.f. which gives . The relative error using the global *a posteriori* error estimate in Bao *et al.* [32] and the DWR method can be compared by defining . Figure 3*a* shows that the DWR method converges faster than the global *a posteriori* error method studied in Bao *et al.* [32]. Also note that the indicative computed error *I*_{err} given in equation (5.6) decreases monotonically which shows the convergence of the DWR method. A relative error of 10^{−6} has been chosen as the stopping criteria because the mesh becomes very irregular after this and will lead to a numerical instability [33,34]. In addition, as the focus is on the local error in the quantity of interest, there is a pollution from the neglected global error. In practice, this problem can be avoided by refining the mesh only in the neighbourhood of the high stress gradients; that is, at the singularities in the geometry. However, the areas of high stress gradients cannot be guaranteed to coincide with the areas of interest and this may lead to pollution error. Hence, the goal-oriented error estimation method could be improved by using a balance between the local error and the global error so that the mesh is refined to assure a high level of accuracy of the quantity of interest and at the same time to not underestimate the effect of the global error [14].

### Example 5.2

Here, the DWR method is compared with the use of a uniform mesh. As above, the relative error using the uniform mesh and the DWR method is defined by . Figure 3*b* shows that the error associated with the DWR method does not decrease monotonically unlike that associated with the uniform mesh when we increase d.f. However, when d.f.>10^{4}, the DWR converges faster and requires fewer degrees of freedom than using the uniform mesh.

### Example 5.3

One of the main reasons for using a finite-element approach is of course its ability to tackle any prescribed geometry. The examples discussed so far have dealt with simple geometries and this of course has enabled a comparison with other results in the literature. In order to develop more sophisticated grating structures that can help push this technology forward, it is necessary to be able to investigate complex geometries. A recent experimental exploration of what could be achieved by freeing up these geometrical constraints has been inspired by a naturally occurring periodic diffraction grating. When scattered by light, the Morpho butterfly wing produces colour and this allows a dynamical control of light flow and wavelength interaction [35] (figure 4). The scattering of the Morpho butterfly wing has inspired applications in biomimetics such as gas sensors, electronic display screens and paints for cars [2,36–38]. This transmitting dielectric grating type has been studied experimentally in Vukusic *et al.* [35] using a focused laser technique to examine the total reflectivity and transmissivity. In this paper, an adaptive finite-element method driven by the DWR is investigated to study the Morpho butterfly using the experimentally measured geometry. Other numerical attempts have been conducted to study the complex geometry of the Morpho butterfly wings [39] but they have had to compromise on the geometry. In those studies, the butterfly geometry was approximated by a series of rectangular blocks stacked periodically in the vertical direction. In addition, a poor agreement of the computed transmitted energy with the experimental result was reported for larger wavelengths [39]. In this paper, one period of the butterfly wing image of the Morpho didius in Vukusic *et al.* [35] was digitized, the coordinates of the grating interface extracted and a Savitzky Golay filter [40], pp. 183, 644–645 used to smooth the resulting geometry (figure 4*b*).

The wavenumber inside the butterfly wing was set as *k*_{2}=(1.57+0.06*i*)2*π*/0.455, the refractive index is *n*_{2}=(1.57+0.06*i*) and outside the butterfly wing as *k*_{1}=2*π*/0.455; these belong to the range of values given in Vukusic *et al.* [35]. A TM field (Case 2B) was imposed, as described by equations (2.10) and (2.11), and the DWR method was used. Therefore, equations (3.6) and (4.5) were used to find the dual solution and the error bounds. The degrees of freedom used was d.f.=11557, with *M*=15 Fourier terms and the degree of the polynomial basis *p*=3. A TOL of 0.001 was used in the algorithm in §5*a* to reproduce the imaginary and real parts of the magnetic field (figure 5). The final indicator computed error *I*_{err} as given in equation (5.6) shows that the discretization error in solving the Helmholtz problem was kept under *I*_{err}=3×10^{−4}. The spectra of the transmission and reflection coefficients, in air (solid line) and in isopropyl alcohol (dashed line), corresponding to one layer (a single period is shown in figure 4*b*), are shown, respectively, in figure 6*a*,*b*. A good quantitative comparison with the experimental results (dots) in Vukusic *et al.* [35] corresponding to single iridescent ground scales in air is obtained. The reduced absolute reflection in isopropyl alcohol compared to that in air, as reported in Vukusic *et al.* [35], is also observed. Unfortunately, there were no measured values given for this latter case in order for a quantitative comparison with our theoretical predictions to be made.

## 6. Summary and conclusion

This paper considers the use of finite elements to numerically solve the Helmholtz equation in periodic diffraction gratings. There are two ways of establishing an *a posteriori* error estimate which are the global error estimate and the goal-oriented error estimate. This paper considers the latter approach and uses a DWR method. A quantity of interest *Q*, was defined and shown to be a linear continuous functional. This afforded the formulation of the dual problem and, when combined with the direct problem, was used to establish an upper bound for the error estimate in theorem 4.1. The evaluation of the error in the functional *Q* represents the primary output from the model, namely the diffraction efficiency.

In the study of Lord & Mulholland [8], the *a priori* error estimate was studied to guarantee that the discrete solution will converge to the exact solution provided that we keep the interpolation error small with respect to the wavenumber *k*. However, as the upper bounds provided in the *a priori* error estimate depend on the exact solution, they cannot be computed. Hence, numerical implementations are limited to using uniform meshes which are computationally expensive. By contrast, for the *a posteriori* error estimate, the upper bounds given by equation (4.5) can be evaluated. The discrete dual problem *z*_{h} is found in the same way as that used to find the discrete direct solution *U*_{α,0h} in equation (3.6). For the exact solution of the dual problem *z*, because it is not known analytically, it can be approximated either by solving the dual problem numerically in a very fine mesh or by increasing the polynomial order using the same mesh as *z*_{h}. Hence, the error in the targeted quantity *Q*(*e*_{h}) can be estimated from the local contribution of each error indicator *ρ*_{K} and *w*_{K} as defined in equations (4.6) and (4.7). In fact, these error indicators are the cell residuals *ρ*_{K} multiplied by the weights *w*_{K} taken from the computed solution. The *ρ*_{K} in turn consists of the field equation residual (*R*_{h}(*K*)) and the flux residual (*r*_{h}(*E*)) which indicate the smoothness of the discrete solution. The weights *w*_{K} capture the influence of the cell residuals on the targeted error *Q*(*e*_{h}) since, differentiating *Q*(*e*_{h}) with respect to *ρ*_{K}, gives *w*_{K}.

This error estimate affords an automatic mesh adaptation based on the local error indicators *ρ*_{K} and *w*_{K}. The DWR method uses these error indicators to maximize the accuracy of the computed diffraction efficiency by choosing a TOL and demanding that as the problem is solved. This can be used to optimize the computational efficiency. So, rather than having a uniform mesh, where all elements of the mesh are refined at each step, the refinement is only performed where the error indicators are large; keeping the coarse mesh where the error indicators are small.

A pseudo-code implementation of this adaptive finite-element method driven by the DWR for solving the periodic grating Helmholtz problem was also described. The merits of using such method were then illustrated for particular examples. In §5*b*, it was shown that the DWR method achieved an acceptable accuracy, converged faster and required fewer degrees of freedom than the global *a posteriori* error estimate proposed in Bao *et al.* [32].

Finally, the geometrical freedom that the finite-element method allows was explored to study a naturally occurring, periodic diffraction grating in the form of a butterfly wing. Its diffraction properties have previously been experimentally measured but this is the first attempt to mirror those results using a finite-element approach. The TM field was successfully calculated and showed the complex interaction between the scatterer and the field. Good quantitative comparison of the experimental results in Vukusic *et al.* [35] with the numerical results presented using here the DWR method for the reflectivity and transmissivity spectrum were also obtained.

There are still some open questions regarding this approach of numerically solving these diffraction problems. It would be interesting to investigate numerically the sensitivity of the grating efficiency to small changes in the geometry of the grating profile. This would be important when considering the inverse problem where the goal is to achieve a desired grating efficiency spectrum by constructing the grating profile that would give rise to this spectrum. It would also be useful to derive an analytical relation to justify the choice of parameters such as the number of the Fourier terms, the size of truncated domain *B*, the order of the polynomial basis and the mesh size *h*. It could be shown for example that the accuracy of the numerical results is not significantly affected by the truncation of the DtN operators.

- Received March 12, 2013.
- Accepted August 29, 2013.

- © 2013 The Author(s) Published by the Royal Society. All rights reserved.