Padé approximants of the Born series of electromagnetic scattering by a diffraction

We present the realization of a vectorial perturbation method based on the Born series applied to strong electromagnetic scattering problems. We present the general theoretical formalism and show a semianalytical implementation for scattering by diffraction gratings. We are particularly interested in the strong scattering regime, where the Born series is known to wildly (namely, exponentially) diverge. By applying Padé approximation to the vectorial Born series, we are able to obtain accurate results from divergent Born series. The method we present has the inherent beneﬁt of being close to the actual physical mechanism behind the formation of a scattered signal, as the solution is built step by step from a sequence of multiple-scattering events. This helps in the understanding of signal formation, which is a key element in inverse scattering problems that are relevant to optical metrology, among others.


I. INTRODUCTION
The analysis of periodic structures is important in a variety of fields: metrology in the semiconductor industry, photovoltaic applications, spectroscopy, and more [1][2][3].The diffractive properties of periodic structures, such as diffraction efficiencies, can be analyzed using tools that solve the forward electromagnetic scattering problem.There are many tools, ranging from general numerical methods, such as the finite-element method (FEM) or the finite-difference timedomain method, to specific numerical methods for periodic media such as the transfer-and scattering-matrix method [4,5] and rigorous coupled-wave analysis (RCWA) [6].Most of these methods have been extensively studied and they have been extended to handle many different cases.On top of that, numerical issues have effectively been addressed throughout the years, such as the convergence of Fourier series in RCWA [7,8].In general, however, these methods lack giving the user insight into how the underlying physical mechanisms generate the scattered signal.In a sense, they often are black-box solvers that yield the solution at once without giving an indication of what part of the structure would have contributed to what part of the signal.
A different way to approach scattering problems is by using perturbation methods.Perturbation methods such as the Born series give more insight in the scattering process because in building the approximation the relevance of different scattering orders is automatically obtained.Moreover, perturbation theory has the advantage that it turns a complicated problem, which requires solving a large system of coupled equations, into a sequence of much simpler problems which are straightforward to solve.In other words, one starts with the known solution to an approximate unperturbed problem and adds corrections to it to increasingly approach the solution of the actual problem.Each correction describes a certain type of interaction between the perturbation and the solution of the unperturbed problem, which for electromagnetic scattering is a type of scattering event.
The process of higher-order scattering, i.e., scattering beyond the first Born approximation, can be interpreted as equivalent to near-field illumination, where the near field is the scattered field due to the incident field or due to the scattered field of one order less.Such near-field illumination can have high spatial frequencies and if these spatial frequencies also have high amplitude this can lead to the scattered far field being sensitive to subwavelength features in the sample, i.e., to superresolution.
Regarding the computation of the Born series, the series rarely converges in cases of practical interest, which severely limits its application.Overcoming the (often wild) divergence has thus been a subject of research.For instance, Osnabrugge et al. developed a modified Born series with guaranteed convergence [9], which has further been applied to inverse problems and extended to the vector and the anisotropic case [10][11][12].Although convergence is guaranteed, it inherently needs a large computational domain extending beyond the scattering objects, which in practice can lead to high computer memory usage.Other approaches to deal with the divergence have been explored as well [13,14].For instance, the change of variables made in Ref. [13] effectively reshuffles the terms of the Born series to make it convergent.However, there is no clear procedure for choosing the parameter used in the change of variables.
In previous work of ours, we chose to retain the original Born series, but used Padé approximation to extract the solution from the series [15,16].Even in cases of divergence, the approximants were able to represent the solution accurately.We considered scalar scattering only, and thereby tested the method on two analytical problems for which the scalar description is rigorous.The method has been applied by other authors as well, for instance, to the problem of scattering by an inhomogeneous dielectric slab, where it was used to find transmission and reflection coefficients instead of the electric field on a grid [17,18].These investigations of the Born-Padé method have also been limited to scalar scattering problems.
It is important to note that, although the Padé approximant has to be computed for all points of a sufficiently fine mesh inside the scattering objects, these computations are simple and can be done easily in parallel and that there is no large system of linear equations to be solved, as is the case in all other rigorous methods for solving scattering problems.
In this work we present the fully vectorial formulation of the Born series applied to the vectorial Lippmann-Schwinger equation.More specifically, we consider the scattering from one-dimensional (1D) periodic structures for which we extend the Born-Padé method from a scalar to a vectorial formulation.By exploiting the periodicity of the problem, we derive a specific implementation of the vectorial recurrence relation for the terms in the Born series.Different formulations of the vector Born series have been presented before [19].
A forward model that incorporates multiple-scattering effects can be of interest for inverse problems.Experiments such as by Simonetti [20] have shown that subwavelength structures can be observed by accounting for multiply scattered waves.If multiple scattering is ignored, i.e., when using the single-scattering or first Born approximation, the resolution becomes diffraction limited.For example, the single-scattering approximation is used in ptychography where it is also called the multiplicative approximation.Examples where multiple scattering is used to achieve subwavelength resolution can be found in Refs.[21][22][23][24][25].A major reason for sticking to the first Born approximation is speed, because in an inverse problem often the forward problem has to be solved many times.
We formulate the vector Born series and show how to apply Padé approximation to the vector series in Sec.II.We develop a semianalytical implementation for 1D periodic structures in Sec.III.The method handles many cases: s or p polarization, classical or conical incidence, and dielectrics or metals.We present simulation results for two gratings in the strong-scattering regime in Sec.IV.We summarize and discuss our conclusions in Sec.V.

II. VECTOR BORN-PADÉ METHOD
We are interested in scattering by an object which has isotropic relative permittivity ε r (r) and which can be spatially inhomogeneous.The object is spatially bounded, though it can also be infinitely periodic in one or more dimensions.The scatterer is surrounded by a homogeneous background, which we choose to be vacuum or air (ε r = 1); k 0 is the wave number in the homogeneous background.Furthermore, no free charge sources are present.We consider time-harmonic electromagnetic fields with time dependence e −iωt (which we will not explicitly write in all formulas from now on), where ω > 0 is the angular frequency.We allow ε r to be complex with a positive imaginary part.Also, we assume negligible magnetization at optical frequencies so that μ r = 1.
In the scalar forward problem, (a component of) the electric field U satisfies the Helmholtz equation where ε r (r) = ε r (r) − 1 is the relative permittivity contrast with respect to the homogeneous background (so outside the scatterer it vanishes).The scalar description is rigorous in certain problems only, e.g., scattering by an infinitely long cylinder [26].For others, its use is often motivated by Born and Wolf's argumentation [27].It starts by considering which can be derived from Faraday's and Gauss's laws.If the permittivity varies slowly enough, the last term on the left is negligible.However, as soon as the permittivity makes a jump at an interface between two materials, this reasoning cannot hold, as the gradient of the permittivity will generate a δ function.Only when the electric field is perpendicular to the normal of the interface is the term zero.Thus, as long as there is a jump in the permittivity and the possibility of oblique incidence, the vectorial formulation is required.For scattering by a diffraction grating, this term is zero for a TE-polarized incident plane wave in a classical mount, for which the electric field is polarized along the grating grooves.In a conical mount, this term is not generally zero and thus a rigorous description is necessarily vectorial.We can derive an equation similar in structure to Eq. (1) for the vector case, containing the electric field E and the vector potential A as unknowns.If we were to apply perturbation theory in the scalar case, we would expand the field U (r) as a power series.For the vector case, we will expand both E and A as a power series.

A. Vector Born series
First, we define the vector potential for the H field instead of the B field, which is possible because we neglect magnetization: For a homogeneous medium, it can be shown that each of the Cartesian components of A(r) satisfies the homogeneous Helmholtz equation [28,29].In the inhomogeneous case of an object placed in vacuum or air, we can derive that The reader can refer to Appendix A for details of the derivation of Eq. ( 4).We observe that the components in Eq. ( 4) are decoupled from each other while the right-hand side acts as a source.Just as the scalar Helmholtz equation (1) has an integral solution, namely, the Lippmann-Schwinger equation, Eq. ( 4) also has an integral solution.This solution is in terms of a convolution of the scalar free-space Green's function G 0 (r; r ) with the induced current density where A (i) is the vector potential of the incident field [a homogeneous solution of Eq. ( 4)] and G 0 (r; r ) is The second term on the right-hand side of Eq. ( 5) is the vector potential A (s) of the field scattered by the object.Thus, the vector potential of the total field A(r) is the sum of that of the incident and the scattered field.We see that A (s) depends on A itself through E under the integral.When using perturbation theory to solve Eq. ( 5), one uses the incident field as an initial solution and adds corrections to it to increasingly approach the exact solution.The series that is formed by the initial solution plus the corrections is the Born series.Mathematically, the series can be obtained by introducing a dummy (unitless) parameter σ .First, we insert σ in Eq. ( 4) by multiplying the source on the right-hand side by σ .The perturbation in our problem is this right-hand side term whose strength can be tuned by varying σ .If we set σ = 0, we get the unperturbed problem of the homogeneous background, while σ = 1 gives the problem that we are eventually interested in.Second, we simultaneously expand E(r) and A(r) as power series in σ : This is different from the scalar case, where we would only expand in series one field U (r).We substitute these representations of E(r) and A(r) in Eq. ( 5).Collecting equal powers of σ , we obtain a set of equations for the terms in each series.The zeroth-order terms are E (i) (r) and A (i) (r), which are related through the relation For the higher-order terms ( 0) we find the two-step recurrence relation We can thus compute the ( + 1)st order by first calculating A +1 from E and then calculating E +1 from A +1 .We start the recurrence with E (i) , which we recall is the electric field of the incident plane wave.Finding A (i) is not necessary, because all higher-order terms can be computed starting from E (i) .
We see that A +1 can be calculated componentwise from E with Eq. ( 9) because no mixing between different Cartesian components occurs.When calculating E from A with Eq. ( 10), components mix due to the double derivative.The use of both E and A in the Born series has the advantage, as in the scalar case, that the integration is a convolution with the scalar Green's function, which has a 1/r singularity.Our derivation thus does not need the dyadic Green's function, which has a 1/r 3 singularity (although that singularity could be integrated analytically [30]).On the other hand, the double derivative still can cause trouble by creating singularities when computing the derivatives.
We note that the recurrence relation in Eqs. ( 9) and ( 10) can be combined into a single calculation step, such that E +1 can directly be calculated from E , where I is the unit dyad.Still, we prefer to explicitly show the step through the vector potential A +1 because it better clarifies at which phase of the scattering process the mixing of the components of the electric field (namely, depolarization effects) actually takes place, given a certain input polarization of E (i) .Some notes on the properties and the validity of the presented Born series follow.Objects can be either bounded or infinitely periodic along one or several Cartesian coordinates.Perfect electric conductors cannot be handled, as it would not be possible to perform the integration in Eq. ( 9) for infinite permittivity, except when the perfect conductor can be taken into account in the Green's tensor.Furthermore, the derived formalism satisfies the Lorenz gauge, which makes it automatically Lorentz invariant [31].

B. Padé approximation
When applied to actual cases of interest, the Born series often diverges.In that case, one can still obtain a valid solution by resorting to Padé approximants, as presented in Ref. [15].In short, instead of simply adding up the terms of the series one by one, one computes Padé approximants of the Born series.The definition of a Padé approximant of order M/N is where the Padé coefficients A m and B n are in total M + N + 1 unknowns.The idea is to match two representations of the solution, the Born series and the Padé approximant, inside the sufficiently small disk in the complex plane of σ where both the Born series and the Padé approximants converge.This yields a system of equations for computing the A m and B n coefficients from the first M + N + 1 terms of the Born series.In other words, we compute the A m and B n coefficients of P N M such that it is identical to the first M + N + 1 terms of the Born series.For more background information on Padé approximation, see Ref. [32].For details on how we compute the A m and B n coefficients, refer to our previous paper [15].We only consider symmetrical Padé approximants, those with a numerator and denominator of the same order [32], i.e., we choose M = N.In the end we have to set σ = 1 in the approximants to obtain the approximate solution to our problem.
In the present problem, the terms of the Born series are vectors.Still, the series is a componentwise summation, so we calculate Padé approximants for every component of E individually.The computational load thus increases by a factor of 3 compared to the scalar case.If some component of the incident E field is zero, it is possible that all the other terms of this component in the Born series are also zero.We then skip the Padé approximation of that component to reduce the computational load.Moreover, it is important to note that, although the Padé approximants have to be computed for sufficiently many points inside the scattering objects, each of these computations is simple and no large system of linear equations has to be solved.
For the convergence of Padé approximants there is no general rigorous theory, like there is for Taylor series.So-called Stieltjes series are the exception, for which it is guaranteed that the sequences of approximants P N N and P N N+1 have a limit for N → ∞.Moreover, they decrease and increase monotonically, thus giving upper and lower bounds to the approximation, respectively [32].If the limits are equal, the limit is a unique Stieltjes function.The coefficients of a Stieltjes series must be real, but since the Born series is generally complex valued, one should look at the real and imaginary parts of the Born terms separately to determine whether the series is of the Stieltjes type.In general, the Born series is not a Stieltjes series, since it is possible to show a counterexample for scattering by a 1D slab, which has an analytical solution and where the real part of the Born terms is not a Stieltjes series.

III. APPLICATION TO A 1D PERIODIC GRATING
The vector Born series can be further developed for the important case of 1D periodic structures.We will derive a specific implementation of the recurrence relation in Eqs. ( 9) and ( 10) by exploiting the periodicity of the problem, by which the permittivity and the electric field can be represented as a Fourier series.Because of the periodicity in x and the invariance of the problem in y, the 3D integral in Eq. ( 9) reduces to a sum of 1D integrals, as we will show.Multiple steps in the computation are done analytically, such as the derivatives in Eq. ( 10), which prevents potential numerical problems.

A. Scope of the 1D periodic problem
Consider a structure which, with respect to a Cartesian coordinate system (x, y, z), is infinitely periodic along x with period p and invariant along y.Along z, the structure is bounded; it is contained within z min < z < z max .The structure can have an arbitrary, finite-value complex relative permittivity profile ε r (x, z) within the unit cell, so besides homogeneous stratified and grating layers also gradient index layers are allowed.In any case, the structure should be surrounded by the same lossless background medium above and below.If the background medium were different above and below the grating, one could accommodate the Green's function for it, but that is outside the scope of this paper.Figure 1 schematically shows an example of a grating consisting of multiple layers and the defined coordinate system.The z coordinate increases downward (as we consider the incident field to be propagating downward), and because the x coordinate increases to the right, the y coordinate increases out of the paper.
A plane wave of wavelength λ in vacuum is scattered by the grating.It has spatial frequency q 0 = 1/λ and real wave FIG. 1. Unit cell of an example of a diffraction grating with the definition of the coordinate system (x, y, z).The +y direction points out of the paper.Although this example is piecewise homogeneous, gradient-index media are also allowed vector q (i) = (q (i) x , q (i) y , q (i) z ), where with q ⊥ = √ q 2 x + q 2 y .The plane wave has a polarization given by the polarization unit vector ê(i) (q (i)  x , q (i) y ).For an spolarized plane wave we define the electric field perpendicular to the plane of incidence, i.e., perpendicular to the (x, z) plane.The electric field of a p-polarized plane wave lies in the plane of incidence.Their polarization unit vectors are, respectively, ês (q x , q y ) = q y q ⊥ , −q x q ⊥ , 0 for q ⊥ = 0, ( êp (q x , q y ) = λ −q x q z q ⊥ , −q y q z q ⊥ , q ⊥ for q ⊥ = 0. (15) For normal incidence (q ⊥ = 0), we define the two as ês = (0, −1, 0) and êp = (−1, 0, 0).For reference, for a classical mount (q y = 0) we have ês = (0, −1, 0) and êp = λ(−q z , 0, q x ).All in all, the electric field of an s-or a ppolarized plane wave with amplitude W (i) is written as

B. Periodicity and diffraction orders
Next we look at how the periodicity of the problem in x determines the form of the solution.The Floquet-Bloch theorem states that the solution is pseudoperiodic in x, i.e., it is periodic up to a phase factor e iαx , where α = 2π q (i) x is determined by the incident wave vector [28].The Floquet-Bloch theorem holds more generally for an incident field that is pseudoperiodic, i.e., periodic except for a phase factor.Furthermore, the periodic permittivity contrast can be represented by a Fourier series with Fourier coefficients ε (m) r (z): It is easy to see that also every term in the Born series is pseudoperiodic in x, i.e., also satisfying the Floquet-Bloch theorem.This reduces our problem of computing the Born series to finding Fourier coefficients E (m) +1 (z) given E (m) (z) for every 1.Higher-order Born terms (higher than the two lowest terms) can be considered to be generated by an illuminating field given by the highest preceding Born term.The latter contains high spatial frequencies in the near field because of the presence of the scattering object.If these high spatial frequencies have high amplitude as in the case of a (near) resonance, then this local illumination can cause that highspatial-frequency information about the scatterer to influence the scattered far field, i.e., it can cause superresolution.

C. Born recurrence relation
We will now derive an expression for computing the Fourier coefficients E (m)  +1 (z) given E (m) (z).Starting from the general recurrence relation ( 9) and (10), we observe that the integral in Eq. ( 9) is a convolution of the Green's function with the product ε r E .In our method, we will evaluate the integral in discrete Fourier space, so we need to find the Fourier coefficients of ε r E first.We can either compute the Fourier coefficients of ε r and E separately and get the Fourier coefficients of the product by convolution in Fourier space or multiply the two in real space and then directly get the Fourier coefficients of the product.The former option would require the appropriate convolution rule as we deal with discontinuous functions [8].We therefore choose the latter option and directly expand ε r E at the beginning of each Born iteration.Since for every z, x → ε r (x, z)E (r) is pseudoperiodic, there exist coefficients c (m) (z) such that x +m/p)x e i2πq (i) y y .(19) The coefficients are computed numerically at the beginning of every Born iteration given the output E (r) of the previous iteration step.For the incident plane wave we have c (0) 0 (z) = W (i) e (i) and c (m)  0 (z) = 0 for |m| > 0. To evaluate the integral over the Fourier components, we use the Weyl expansion to expand the Green's function in plane waves [33]: This will reduce the 3D integration in Eq. ( 9) to a 1D integral over z, a Fourier sum in x, and no summation in y.The rest of the derivation of the Born recurrence relation is done in Appendix B; we will only state the final result here.Before we do that, we define some wave-vector quantities for every diffraction order m: q xm = q (i) x + m/p, q y = q (i) y , q ⊥m = q xm x + q y ŷ, and Furthermore, part of the derivation is the integration over z, which can be split in two parts, corresponding to whether z < z or z > z in Eq. ( 20): These two integrals depend on z by the exponential factor e ±i2πq zm z (which can be pulled out of the integral), but z also appears in either of the integration limits.Notice that for every m the pair of integrals in Eqs. ( 22) and ( 23) is equivalent to a Born iteration in a 1D (aperiodic) scattering problem [15].The final result for the Fourier coefficients of E +1 is Equation (24) shows that we can immediately calculate the next Born order of E without calculating A first.We see here that the dependence of the terms on z by e ±i2πq zm z corresponds to up-or downward traveling waves and their amplitudes depend on the z position within the structure.For z < z min , C (m) ,< (z) = 0 and hence there are only waves with e −i2πq zm z , while for z > z max we have C (m)  ,> (z) = 0 and then all waves depend on z by e i2πq zm z .This makes sense as scattered waves only travel away from the structure.Furthermore, we see that we can derive the double derivative in Eq. ( 10) analytically, removing the need to do it numerically.

IV. RESULTS
The formalism we developed in Sec.III is general: It is valid for any incident plane wave and any number of layers.We show results for two cases of s polarization: one where scattering is in the moderate-scattering regime and one for strong scattering.For p polarization, slow convergence of Fourier series occurs, which requires special methods to accelerate; however, that is beyond the scope of the present work.We compare the results of each case with the FEM solution in COMSOL MULTIPHYSICS.

A. Details of the numerical implementation
As the problem is invariant in y up to a phase factor, we solve the problem in the (x, z) plane for y = 0. We discretize one period of the grating on a 2D Cartesian grid, which is chosen uniform in both x and z.In the z direction, the grid extends 250 nm above and below the structure, although strictly we only need to integrate over the structure itself to compute the Born series.The computation of the Born series starts at the zeroth-order field: a plane wave with amplitude W (i) = 1.Thus, the iterative procedure starts with E (i) (x, z) of the plane wave and the relative permittivity contrast ε r (x, z) of the grating as the input.

Born iterations and the Gibbs phenomenon
At every iteration step, we compute the truncated Fourier representation of the product ε r E .We use the maximum number of Fourier components as determined by the grid spacing per the Nyquist theorem: N Fourier < p/2 x.Then Eq. ( 24) is evaluated to obtain the Fourier coefficients of the next Born order.Before the start of the next iteration, we have to sum the Fourier series to return to real space.To counteract the Gibbs phenomenon at this point, we sum the Fourier series using the Fejér kernel [34].Essentially, we take the average of the sequence of partial Fourier sums, which is a Cesàro sum.Denoting the partial Fourier sum up to order m = ±M by E (M )  +1 , the Cesàro sum is In the limit of M max → ∞, it converges uniformly to E +1 (x, z).We use Eq. ( 25) as input to the next iteration.

Quantifying the divergence of the Born series
When the scattering is weak, the Born series converges.In the strong-scattering regime, the terms of the Born series grow faster the stronger the scattering is.We can get an indication of the rate of convergence or divergence and therefore the scattering strength by looking at the maximum of the modulus of every Born term E ,y : If we suppose that A grows exponentially, i.e., A = Born A 0 , with Born the growth rate, then we can estimate Born from the Born series up to the Nth Born term by We have weak and strong scattering (convergence and divergence) for Born 1 and Born 1, respectively.Around Born = 1 this parameter is not sufficient to discriminate between a converging and a diverging Born series.

Padé approximation
After computing enough Born terms, we calculate the Padé approximants of the Born series for each point on the 2D grid.To find the approximant P N N , one needs 2N + 1 terms of the Born series computed by performing 2N iterations.Finding the coefficients of the Padé approximant of order N/N then consists of solving a linear system of 2N + 1 equations.Note that the Padé coefficients have to be calculated separately for each Padé approximant of a different order.

Optimal Padé approximant
Beforehand, it is not known which order Padé approximant is necessary to get a desired accuracy.So we could benefit from knowing at what rate the approximants are converging towards a result and at what point the improvement starts to level off.To get an indication of how quickly the approximants are converging as N increases, we can look at the relative difference of an approximant of order N compared to the preceding approximant of order N − 1: There is one major limiting factor for computing Padé approximants: machine precision.The higher the order of the Padé approximant, the larger the system of equations needed to compute the Padé coefficients.The condition number of this system increases, requiring higher precision.Beyond a certain Padé order N, the improvement of the approximation levels off, given that we use double-precision complex floating point numbers.This is the maximum available precision in scientific computing software packages like NUMPY.If instead we use single precision, the approximation levels off at a lower Padé order N and consequently a less accurate result can be obtained.

Comparing the results with COMSOL
For the finite-element solution in COMSOL, we simulate one period of the grating.We use periodic boundary conditions in the x direction and perfectly matching layers in the z direction.The COMSOL mesh has a maximum element size of 20 nm outside the scatterer and 10 nm inside.To make the comparison, we linearly interpolate the result from the FEM mesh to the rectangular grid used in our method.We quantify the quality of each approximant of E y relative to the COM-SOL result E COMSOL by determining the relative error over the scatterer, defined as

C. Example 1: Moderate scattering
First we consider a grating of SiO 2 , which has a low relative permittivity ε r,SiO 2 = 2.18 at λ = 633 nm [35].The grating has a pitch of 600 nm and consists of a rectangular grating layer on top of a planar layer.The grating layer has a thickness of 200 nm and a linewidth of 300 nm, while the planar layer is 150 nm thick.See Fig. 2 for | ε r (x, z)| for this example.
The incident field is an s-polarized plane wave in a classical mount, so only E y is nonzero and E x and E z are zero in both the incident and the scattered field.The plane wave comes in obliquely with q (i) x = 0.5q 0 , so q (i) z = 0.87q 0 .Its polarization unit vector is ê(i) s (0.5, 0) = (0, −1, 0).Furthermore, we use a grid spacing of x = z = 10 nm.According to the Nyquist theorem, we can use up to 60 Fourier components.The first step is to compute the Born series.The modulus of E y of the = 1, 2, 3 terms of the Born series are shown in Fig. 3 using a normalization factor displayed in the bottomleft corner.This factor indicates that the series diverges.A better indication is given by Eq. ( 27), which yields Born = 1.55.While this shows the divergence, it however does not satisfy either Born 1 or Born 1, so we cannot say with certainty that it is a case of weak or strong scattering.
From the Born series, we computed Padé approximants up to order 30/30, which requires 61 Born terms.Figure 6 shows how the error (29) of P N N relative to the COMSOL solution decreases with increasing N, leveling off after P 5  5 .However, if we choose our optimal approximant by looking at how the difference between subsequent approximants develops according to Eq. ( 28) (without knowledge of the COMSOL solution), we find that the convergence of the Padé approximants only levels off after P 13  13 (see Fig. 6).We show this optimal approximant in Figs. 4 and 5, as well as the COMSOL solution and the difference between the two, in terms of amplitude and phase.Padé approximation achieves an error of about 2% in terms of the modulus of the solution as seen in Fig. 4.  The error (29) of the approximants relative to the COMSOL solution levels off earlier as a function of N than the error (28) relative to the preceding approximant.This could be caused by multiple factors, such as the difference in the discretization of the grating in our simulation and in COMSOL.Decreasing x and z in our simulation decreases this difference such that the relative error levels off later.Furthermore, as seen in Fig. 6, the error ( 28) levels off at a lower N when using singleprecision floating numbers instead of double precision.This illustrates that the precision used determines the maximum achievable accuracy with Padé approximation.

D. Example 2: Strong scattering
Now we consider a case that exhibits stronger scattering: a grating with high-permittivity materials.The grating has period p = 450 nm and has four layers.The top layer is a rectangular grating layer of silicon (80 nm thick and 180 nm wide).Underneath it are planar layers of SiN (60 nm), gold (40 nm), and SiO 2 (80 nm).
Again, we first compute the Born series; the modulus of E y for the first Born orders = 1, 2, 3 is shown in Fig. 8.The plots show the normalized modulus of every term, with the normalization factor displayed in the bottom-left corner of each plot.It can be seen that the series diverges.The parameter defined by Eq. ( 27) now is Born = 6.74, which means this case is far into the strong-scattering regime.
While the Born series wildly diverges, Padé approximation can retrieve the complex electric field from the series.The optimal approximant is P 13  13 , shown in Figs. 9 and 10 in a comparison with the COMSOL simulation.To compute P 13  13 , 27 Born orders are needed, of which the highest Born order has amplitudes which are 10 21 times larger than the incident field.
The approximants are more accurate further away from the grating than inside or near it.Some artifacts are also visible, e.g., noiselike artifacts which are signs of the limits of machine precision when solving for the Padé coefficients.We observe them in higher-order approximants only.
The error relative to the COMSOL solution is shown in Fig. 11 for Padé approximants of E y up to P 30  30 , showing that the approximation does not improve further after P 8  8 .Looking at the relative difference between P N N and P N−1 N−1 according  to Eq. ( 28), the approximants level off after P 9 9 , as seen in Fig. 11.The minimum of Eq. ( 28) is at N = 13, P 13  13 is the optimal approximant.For higher-order approximants, this relative difference starts increasing slightly again.In Fig. 11 we also show the error (28) when using singleprecision floating numbers, denoted by c64.Again, the error levels off at a lower order N, now not even nearing the lowest error when compared to COMSOL.

V. CONCLUSION
We have developed a vectorial perturbation method based on the Born series to solve the forward problem of electromagnetic scattering by a diffraction grating.We presented a general formalism, as well as an implementation of the Born recurrence relation for the grating problem.To handle cases in the strong-scattering regime, where the Born series diverges, we computed Padé approximants of the Born series.This is done for each electric-field component separately, for every point on the coordinate grid.
The method requires solving relatively small systems of equations for the Padé coefficients, which moreover can be done in parallel at each step of the iteration.There is no need to solve a large system of linear equations as is the case in other rigorous methods; hence the requirements on computer memory are low.
We demonstrated the method in two examples of scattering of an s-polarized plane wave by a grating: one exhibiting moderately strong scattering by a dielectric and one exhibiting strong scattering by a grating with dielectric and metallic layers.The method is general: It allows for a classical and conical mount, s and p polarization, and any configuration of grating and planar layers with finite relative permittivity.For p-polarized scattering, attention must be paid to the slow convergence of the Fourier harmonics, a known issue in other methods as well.We leave consideration of p polarization for future work.Machine precision limits the method such that higher-order Padé approximants do not further converge towards the exact solution beyond a certain Padé order, when noiselike artifacts begin to appear in the approximants.
The method can be extended to 2D periodic structures, which would involve a second Fourier summation, in y.This does not affect the Born series formalism presented Sec.II, but it requires modifications in Sec.III, taking into account the second summation.Furthermore, although the problem is only solved for a single incident plane wave, a solution could also be obtained for other incident fields.This would require decomposing it in Bloch waves, repeating the calculation for every Bloch wave, and coherently summing the results.
The presented Born-Padé method can provide insight into how the underlying physics forms the solution because the Born series builds the solution step by step, in terms of singleand multiple-scattering effects.Also when computing Padé approximants, one exactly knows which Born orders are used in the calculation, thus tracing which scattering effects are contributing.A perturbation method like the one presented here could be of interest for the design of periodic struc-tures, as it contributes to the understanding of how the parts of a structure affect the generated signal.Furthermore, the Born-Padé method provides insight into the importance of multiple-scattering events.If multiple scattering is strong, this indicates that the scattered far field may be sensitive to subwavelength details.For these reasons, the presented method is very interesting for use in inverse problems as well, with potential applications in, for example, (optical) metrology.
By computing the integrals over q x and q y we obtain A +1 (r) = −iωε 0 ∞ m=−∞ i 4π q zm e i2π (q (i) x +m/p)x e i2πq where q zm has been defined in Eq. (21).The integral over z can be split in two parts: one where z < z and one where z > z, i.e., C (m)  ,< (z) = where q ⊥m = (q (i) x + m/p)x + q (i) y ŷ is the transverse component of the wave vector of the mth order and r ⊥ = x x + yŷ.

FIG. 2 .
FIG. 2. Absolute value of the relative permittivity contrast | ε r (x, z)| for one period of the simulated grating in Example 1.Both layers are SiO 2 with relative permittivity ε r,SiO 2 = 2.179.The pitch is 600 nm.

FIG. 3 .
FIG. 3. Modulus of E y for the Born terms (a) = 1, (b) = 2, and (c) = 3 for Example 1.The terms are normalized so that the same color bar can be used, with the normalization factor displayed in the bottom-left corner of every plot.The grating outline is shown by a white dashed line.

FIG. 6 .
FIG. 6. Convergence of the Padé approximants for Example 1 shown in two ways: the relative error of P N N compared to the COMSOL simulation [see Eq. (29)] and to P N−1 N−1 [see Eq. (28)].The latter is also shown for single-precision computation, denoted by c64.

FIG. 7 .
FIG. 7. Absolute value of the relative permittivity contrast | ε r (x, z)| for one period of the simulated grating in Example 2. The layers are (from top to bottom) silicon, SiN, gold, and SiO 2 .The pitch is 450 nm.

FIG. 8 .
FIG. 8. Modulus of E y for the Born terms (a) = 1, (b) = 2, and (c) = 3 for Example 2. The terms are normalized so that the same color bar can be used, with the normalization factor displayed in the bottom-left corner of every plot.The grating outlines are shown by white dashed lines.