## Abstract

We present an improved adaptive mesh process that allows the accurate control of the numerical solution of interest derived from the solution of the partial differential equation. In the cases of two-dimensional studies, such an adaptive meshing is applied to compute phenomenon involving high field gradients in near-field (electric intensity, Poynting’s vector, optical forces,…). We show, that this improved scheme permits to decrease drastically the computationnal time and the memory requirements.

© 2007 Optical Society of America

## 1. Introduction

Since the advent of the near-field optical microscopy and its use for the characterization of nano-structures, several numerical methods have been developed for the modeling of the interactions between the electromagnetic field and the nanostructures which exhibit high gradients in the electric field. These gradients are of interest for the prediction of Raman or fluorescence in near-field. In principle, all the methods which make possible the computation of the optical near-field around metallic nanostructures illuminated by an external source. Nevertheless, only a few of them are able to deal with high local variations of the electromagnetic field [1, 2]. This condition is directly connected with the numerical convergence of the algorithm and the discretization of the problem. Comparisons of various methods were reviewed in several papers [2, 3, 4, 5]. Among these different methods, the finite element method is a powerful tool that has proven its utility in the study of electromagnetic problems. In the present work, we compare the results to the Mie’s theory in the case of an infinite gold cylinder in order to improve the convergence of the numerical scheme to the solution of interest as the Poynting’s vector amplitude and the intensity of the electric field. We compare and discuss the advantages of the improved adaptive mesh scheme relatively to the classical adaptive remeshing scheme with *h*-method estimator. We show that the use of FEM with this improved adaptive meshing permits to decrease drastically the computational time and the memory requirement.

The organization of the paper is: Section 2 presents the formulation of the problem in near-field. Section 3 is devoted to the description of the improved adaptive mesh process, used to compute the total electromagnetic field and the Poynting’s vector amplitude. In Sec. 4, numerical results obtained with this improved mesh adaptation process in connection with FEM are presented and discussed before concluding in Sec. 5.

## 2. Formulation of the problem in near-field

We study the particular case of 2D structures (nanowires, nanotubes, nanoslits…) considered infinite along *z*-direction. The different media are supposed to be non-magnetic and a p-polarized illuminating wave is used to produce resonances and localized plasmons. Assuming a timeharmonic dependance of the form exp(*jωt*) for the fields [6], the partial differential equation (PDE) to be solved is the Helmholtz’ equation for the magnetic field *H _{z}* [7]:

where *c* is the velocity of light in vacuum and *ε _{r}* is the dimensionless relative complex permittivity of the considered medium (

*ε*=

*ε*

_{0}

*ε*) which is a function of the spatial coordinates (

_{r}*x*,

*y*). The electric field

**E**(

*x*,

*y*) is then deduced from the Maxwell’s equation:

and the Poynting’s vector **P**(*x*,*y*) depends on both fields:

where the ∗ denotes the complex conjugate. The curl operator (∇×.) is computed with finite elements instead of finite differences to avoid a loss in the accuracy. In conventional Finite-Elements approaches, the accuracy of the solution of the PDE (**H**) governs the remeshing process. Therefore a loss of accuracy is expected on **E**, especially if it exhibits high gradients. Moreover, in near-field computations, the discontinuity of the electric components normal to the material surfaces makes such gradients to appear in distinct zones where **H** is varying. In contrast to previous works where the *a posteriori* error estimator is achieved on the PDE solution [8, 9, 10, 11], the present work uses an *a posteriori* estimation of the interpolation error on the solution of interest (**E** or **P**) derivated from the PDE solution (*H*).

## 3. The improved adaptive mesh method

In this section we will describe the principles of the proposed adaptive-mesh scheme to be employed in our numerical experiments. The improvement of this strategy with respect to the basic adaptive meshing approaches lies on the simultaneous utilisation of an adaptive meshing process and an *a posteriori* estimation of the interpolation error in the electromagnetic computation. These two features allow an important reduction in the number of computed nodes.

In Fig. 1 we consider a one dimensional element *e* to outline the main differences between two widely employed adaptive meshing processes and our improved scheme. The basic adaptation process is shown in Fig. 1(a), it consists on dividing by two the length of segment *e*
^{1}
_{1} by adding only one new node and evaluating the error on the PDE solution. The *h*-method estimator is illustrated in Fig. 1(b), and permits to divide the length of the segment with respect to a threshold *δ* for the error estimator based on the PDE solution (*δ _{H}*). The meshing process in Fig. 1(c) corresponds to the improved scheme proposed in this work and permits to divide the length of the segment with respect to a threshold

*δ*for the error estimator based on the solution of interest (

*δ*or

_{E}*δ*,‥).

_{P}The strategy of the global adaptive process is governed by a map mesh element size that ensures the conformity with the underlying geometry domain and of the accuracy of the solution of interest. This size map results from an *a posteriori* estimation of the interpolation error on the field of interest (**E**, *P*,‥).

The main steps of the improved general adaptive remeshing scheme in a domain Ω defined from its boundary Γ can be summarized as follows:

- 0 generation of a first coarse mesh
*ℱ*(Ω) of Ω, - 1 Computation of the field
*H*(PDE solution) by FEM on_{z}*ℱ*(Ω), - 2 Derivation of the solution of interest
*S*(i.e._{ϕ}**E**or**P**,‥) on*ℱ*(Ω), - 3 Estimation of the physical error (i.e. computation of the interpolation error of the solution of interest
*S*(_{ϕ}*ϕ*being**E**or**P**,‥) based on the estimation of the discrete Hessian of*S*; definition of a physical size map_{ϕ}*ℋ*(Ω) allowing to bound this error with a given threshold*δ*),_{ϕ} - 4 remeshing of the domain conforming with the size map
*ℋ*(Ω), - 5 loop to step 1 when threshold (
*δ*) is not satisfied._{ϕ}

To bound the deviation *η _{ℱ}* between the exact solution

*S*and the physical solution

*S*associated with the mesh

_{ℱ}*ℱ*(Ω) (obtained by finite element computation or derivation), an indirect approach is used. This consists in defining

*S*˜

_{ℱ}an interpolation of

*S*on mesh

*ℱ*(which can be a piecewise linear or quadratic function, depending on the degree of the elements of

*ℱ*) and

*η*˜

_{ℱ}the deviation between

*S*and

*S*˜

_{ℱ}. Such an interpolation error

*e*˜

_{ℱ}can be used to bound the gap between the solution obtained numerically and the unknown exact solution:

where ∥·∥ is a norm and *C* a constant independent from *ℱ*. In other words, the finite element error is bounded by the interpolation error of the solution [9, 10, 11, 12, 14, 15]. Therefore, the problem is reduced to the characterization of those meshes on which the interpolation error is limited by a threshold [13]. Also, the examination of a ‘measurement’ of the interpolation error makes it possible to obtain constraints on the size of the mesh element and the interpolation error gives the way to quantify the degree of ‘regularity’ of the exact solution *S*. The remeshing process consists in building a new mesh of Ω by using the interpolation error achieved on the solution of interest *S _{ϕ}*.

The resolution method, based on the local deviation of the solution surface, gives a bound locally in the vicinity of each mesh vertex. This deviation can be quantified simply by considering the Hessian along the normal to the surface (i.e. the second fundamental form of the surface). Let *P* be a vertex of the solution surface. Locally, in the vicinity of *P*, this surface has a parametric representation *ψ*(*u*,*ν*), (*u*,*ν*) being the parameters, with *P*=*ψ*(0,0). Applying the Taylor expansion of order two to *ψ* in the vicinity of *P*, we obtain:

where the derivative of *ψ* are evaluated at (0,0) and *ê* = (1,1,1). If *n*(*P*) denotes the normal to the surface at *P*, then the quantity 〈*n*(*p*)∣(*ψ*(*u*,ν)−*ψ*(0,0))〉 (〈·∣·〉 being the scalar product) represents the gap point *ψ*(*u*,*ν*) to the tangent plane at *P*, and can be written as:

which is proportional to the second fundamental form of the surface for *u*
^{2}+*ν*
^{2} small enough. Hence, to measure the local deviation of the surface, one only has to consider, for each node of the surface, the maximal deviation from the adjacent nodes to the tangent plane at this node. This size at each node of *ℱ*(Ω) is then defined as being proportional to the inverse of this deviation [12]. Formally, the size associated with a node *w* of *ℱ*(Ω) can be written as:

where η(*w*,*S*(*w*)) denotes this new deviation. The size map *ℋ*
_{ϕ} associated with the vertices of *ℱ*(Ω) is deduced.

Due to the optimization of the position of the new vertex in respect to the *a posteriori* interpolation error achieved on the field of interest, the improved adaptive remeshing procedure permits to reduce the number of the iterations and to control the accuracy of the solution. That contrasts to the basic adaptive process where two loop sequences were necessary; the first one on the error on the PDE solution and the second one on the error on the field of interest (**E**, **P**,…). In near-field physics, the physical phenomenon of interest are those implying a strong variation of the electric field **E** (or the Poynting’s vector amplitude ∣**P**∣ or the electromagnetic force **F**), while the variations of the magnetic field **H** (the PDE solution) are smooth. In contrast with the classical remeshing procedure and/or the classical adaptive meshing scheme with *h*-method, the proposed adaptive scheme optimizes the mesh only in the regions where the solution of interest (**E**, or ∣**P**∣, or **F**,…) presents high local variations. Moreover, it permits to control *a priori* the error on the solution of interest (**E**, or ∣**P**∣, or **F**,…) and not on the PDE solution (**H**), as it is done with the classical error estimator.

## 4. Applications to high gradient of the electric field computation

At this stage, we use an example to illustrate the computation of the the electromagnetic field with the improved adaptive mesh scheme. We compare the result of our numerical experiments that we have obtained with a classical remeshing process with *h*-method (*a posteriori* error estimation on the PDE solution *H _{z}*). In the following, the classical remeshing process with

*h*-method is denoted as classical remeshing process.

#### 4.1. Case of the gold nano-cylinder

It consists in the computation of the electric and magnetic field intensities; **E** and **H**, around an infinite circular-cylinder of gold with radius *a* = 15 nm, illuminated with a p-polarized wave at *λ* = 660 nm, as depicted in Fig. 2. The harmonic solution of the problem (Eq. 1) enables the use of the experimental values of the relative permittivities of the materials (*ε _{r}* = −13.6969 +

*j*1.0356) and a non-regular non-Cartesian adaptive mesh to reproduce the curvature of the nanostructure. The FEM adaptive mesh is realized employing the BL2D-V2 bidimensional mesher [16] to match closely the geometry of the object. This simple case is studied through a comparison of the results obtained employing FEM with a classical adaptive meshing, with the improved adaptive mesh scheme and with the analytical solution computed with the Mie’s Theory [17, 18].

Figures 3 show the evolution of the normalized intensity of the electric field ∣**E**∣^{2}/∣**E _{0}**∣

^{2}for the starting coarse mesh (see Fig. 3(a,b) for the classical remeshing and the improved remeshing, respectively), the first and last remeshing (see Fig. 3(c) and Fig. 3(e) for the classical remeshing, Fig. 3(d) and Fig. 3(f) for the improved adaptive mesh scheme). The computational times and the number of nodes at each iterative step for the two remeshing procedures are given in Table 1. This computation time can strongly vary, depending on the resonant behavior of the diffracting structure [19]. Furthermore, the computation time of FEM strongly depends on the number of mesh refinements and of the method used to inverse the matrix obtained from the variational scheme. Therefore the present CPU times are indicative instead of absolute values. In this paper, we use an optimized LU decomposition for sparse matrix [20] or, for larger matrices, an iterative biconjugate stabilized gradient with hermitian preconditioning.

The evolution of the computational mesh as a function of the adaptive remeshing scheme (see Fig. 4(a,c,e) for the classical remeshing and Fig. 4(b,d,f) for the improved remeshing). In contrast with a classical mesher, the refinement and the ‘regularity’ is obvious in the case of the improved meshing (Fig. 4(b,d,f)), which permits to describe accurately the shape of the field intensity. Moreover, the adaptation of the mesh to the solution of interest (i.e. intensity of the electric field) clearly appears in Fig. 5(a) with the improved remeshing scheme in contrast to a classical remeshing process (Fig. 5(b)).

We also show in Fig. 6(a-f) the evolution of the normalized intensity of the magnetic field ∣**H**∣^{2}/∣**H _{0}**∣

^{2}as function of the remeshing for the both the classical and the improved adaptive process. The starting coarse mesh and two successive remeshing are presented in Fig. 6(a,c,e) for the classical remeshing and in Fig. 6(b,d,f) for the improved remeshing. In contrast to the strong variations of the intensity and shape of the electric field, the intensity and shape of the magnetic field are smoothly varying. Therefore an optimization of the mesh seems to be less critical when only the magnetic field is the solution of interest. However, this is rarely the case and that clearly shows the adavantage of an adaptive meshing process based on the solution of interest (

**E**,

**P**) instead of the computational variable

**H**(the PDE solution).

In order to test the remeshing, we focus on the region close to the nano-object where high gradient of the electromagnetic field occurs. We compute the variation of the normalized electromagnetic intensity at short distances of detection from the center of the cylinder along the x-axis where the contribution of the scattered field is supposed to be important. In this region, a fine meshing is necessary for FEM to accurately describe the intensity variations (as illustrated in Fig. 3). Consequently, the remeshing is adapted to the electric field computation (Poynting’s vector amplitude, force vector amplitude,…) instead of the *H _{z}* amplitude. As shown in Table 1, the remeshing process permits to decrease drastically the total number of nodes. We can also note that a distinction must be done between: ‘convergence of the numerical method’ and ‘convergence to the physical solution’. We consider that the result of the computation has converged to the physical solution when the relative error between the result and the analytical solution (when it exists) is smaller than a maximal tolerance (e.g.∥

*δ*∥

_{ϕ}_{∞}≤ 0.5%). when no analytical solution exists, the relative error can be computed by comparing the solutions obtained with two successive mesh refinements or by the interpolation error estimator of the remeshing procedure.

We show in Fig. 7(a,c) the evolution of the normalized (a) electric and (c) magnetic intensity computed with FEM, for various distance from the center of the nano-object (radius *a* = 15 nm) along the x-axis, for successive mesh refinements with classical and improved adaptive remeshing. Figure 7(b,d) shows the error relatively to the Mie’s results as a function of the distance and the mesh refinements. Mie and FEM computations are performed using the same non-Cartesian mesh. In contrast to the classical remeshing, which need four iterations (0-step and four iterative steps), only two improved adaptive remeshing are necessary to stabilize the solution and to assure a convergence of the numerical results to the physical solution (with an relative error *δ _{ϕ}* around 0.5%, see Fig. 7(b)). The smallest cell size is less than 0.03 nm in the region where high gradient of field occurs.

#### 4.2. Case of the gold nano-square

The last illustrative example concerns the study of a more complex diffracting object, as depicted in Fig. 8. It consists of an infinite gold square-cylinder of semi-length *a* = 15 nm deposited on a glass substrate (*ε _{r}* = 2.25) and illuminated by a p-polarized incident field under an angle

*θ*=60

*°*along the y-axis. Strong confinement of the energy (Poynting’s vector amplitude) can be observed at the edges of the square section. Like in the previous example, the use of an optimized mesh refinement processing is crucial to assure the convergence of the numerical result. The used computational domain Ω consists in a square-box of 600 nm × 600 nm and the first step consists in defining a coarse mesh (

*N*= 670) at the boundaries of the subdomains (nanoparticle, substrate).

_{nodes}The remeshing procedure consists in computing the *H _{z}* component, deducing the

*E*and

_{x}*E*and their associated

_{y}*P*and

_{x}*P*components (see Eq. 3). The improved adaptive meshing uses directly the interpolation error on the the Poynting’s amplitude. In the present computations, only five remeshing steps are necessary to converge and stabilize the numerical solution. These remeshing steps (1, 2, 3, 4 and 5) involve

_{y}*N*= 1580,

_{nodes}*N*= 3371,

_{nodes}*N*= 13205,

_{nodes}*N*= 24906 and

_{nodes}*N*= 78175, respectively. Their associated computation times are 2s, 4s, 10s, 46s, 130s and 940s for the coarse mesh (step 0) and the 5 remeshing steps (1 to 5), respectively. We show in Fig. 9(a-f) the maps of the Poynting’s vector amplitude ∣

_{nodes}**P**∣/∣

**P**∣ computed by FEM for three of the five remeshing steps (i.e. steps 2, 3 and 5). The values of the peak amplitude at the corners of the object are directly related to the mesh refinement used and reveal the difficulty to describe a strong gradient of the field near a material geometric singularity.

_{0}Finally, we present the evolution of the amplitude around the nanostructure. Figure 10(a) shows the amplitude along the x-axis, at the interface between the substrate and the nanostructure, and Fig. 10(b) is along the y-axis (i.e. on the lateral face of the square-cylinder), for various mesh refinements. The topic is to compute the value of the Poynting’s vector amplitude (or intensity, force,‥) and to describe accurately the shape of the physical field in the vicinity of the edges of the object. The convergence and the stabilization to the solution as a function of the remeshing clearly appears.

## 5. Conclusion

In this paper we have presented an improved adaptive mesh process based on the *a posteriori* error indicator estimation on the solution. This approach is used in connexion with the numerical method FEM for the computation of the intensity of the total electromagnetic field and/or the Poynting’s vector amplitude in context of near-field enhancement. We have shown that a sufficient grid size refinement must be considered to assure the convergence and stabilization of the solution, particularly in the regions where a strong confinement of the light around the nanostructure takes place. The *a posteriori* error estimator is achieved on the field of interest (**E**, **P**). The main advantage of this *a posteriori* error estimator is to decrease the remeshing number on the PDE solution and to control the error on the solution. We have shown that an adaptive mesh scheme in connection with FEM is suitable to describe accurately the complex geometries and to represent accurate map of the intensity. This improved procedure accelerates drastically the convergence of the solution and minimizes both the memory requirement and the computational time. Moreover, its extension and its application to problems in three dimensions is direct.

## Acknowledgments

Authors gratefully thank Dr D. Macìas for fruitful discussions on the improvement of the paper.

## References and links

**1. **J.P. Kottmann and O.J.F. Martin, “Accurate solution of the volume integral equation for high-permittivity scatterers,” IEEE Trans. Antennas Propag. **48**,1719–1726 (2000). [CrossRef]

**2. **D. Barchiesi, B. Guizal, and T. Grosges, “Accuracy of local field enhancement models: toward predictive models?,” Appl. Phys. B ,**84**,55–60 (2006). [CrossRef]

**3. **D. Barchiesi, C. Girard, O.J.F. Martin, D. Van Labeke, and D. Courjon, “Computing the optical near-field distributions around complex subwavelength surface structures: A comparative study of different methods,” Phys. Rev. E **54**,4285–4292 (1996). [CrossRef]

**4. **B. Guizal, D. Barchiesi, and D. Felbacq, “Electromagnetic beam diffraction by a finite lamellar structure,” J. Opt. Soc. Am. A **20**,2274–2280 (2003). [CrossRef]

**5. **T. Grosges, A. Vial, and D. Barchiesi, “Models of near-field spectroscopic studies: comparison
between Finite-Element and Finite-Difference methods,” Opt. Express **13**,8483–8497 (2005), http://www.opticsexpress.org/abstract.cfm?URI=OPEX-13-21-8483. [CrossRef] [PubMed]

**6. **M. Born and E. Wolf, *Principle of Optics* (Pergamon Press, Oxford, 1993).

**7. **J. Jin, *The Finite Element Method in Electromagnetics* (John Wiley and Sons, New York,1993).

**8. **P. Ingelström and A. Bondeson, “Goal-Oriented error estimation and h-adaptivity for Maxwell’s equations,” Comput. Methods Appl. Mech. Eng. **192**,2597’2616 (2003). [CrossRef]

**9. **P. Houston, I. Perugia, and D. Schotzau, “Energy norm a posteriori error estimation for mixed discontinuous Galerkin approximations of the Maxwell operator,” Comput. Methods Appl. Mech. Eng. **194**,499–510 (2005). [CrossRef]

**10. **D. Pardo, L. Demkowicz, C. Torre-Verdìn, and L. Tabarovsky, “A goal-oriented hp-adaptive nite element method
with electromagnetic applications. Part I: Electrostatics,” Int. J. Numer. Methods Eng. **65**,1269–1309 (2005). [CrossRef]

**11. **D. Xue and L. Demkowicz, “Modeling of electromagnetic absorption/scattering problems on curvilinear geometries using hp finite/infinite element method,” Finite Elem. Anal. Design **42**,570–579 (2006). [CrossRef]

**12. **H. Borouchaki, P. Lang, A. Cherouat, and K. Saanouni, “Adaptive remeshing in large plastic strain with damage,” Int. J. Numer. Methods Eng. **63**,1–36 (2005). [CrossRef]

**13. **M. Berzins, “Mesh quality: a function of geometry, error estimates or both?,” Eng. Comput. **15**,236–247 (1999). [CrossRef]

**14. **M. Ainsworth and J.T. Oden, “A posteriori error estimation in finite element analysis,” Comput. Methods Appl. Mech. Eng. **142**,1–88 (1997). [CrossRef]

**15. **R. Radovitzky and M. Ortiz, “Error estimation and adaptive meshing in strongly non-linear dynamic problems,” Comput. Methods Appl. Mech. Eng. **172**,203–240 (1999). [CrossRef]

**16. **P Laug and H Borouchaki 2003 “BL2D-V2: mailleur bidimensionnel adaptatif,” *Report INRIA RT-0275*http://www-rocq1.inria.fr/gamma/cdrom/www/bl2d-v2/INDEX.html

**17. **G. Mie, “Beiträge zur Optik trüber Medien, speziell kolloidaler Metallösungen,” Ann. Phys. **25**,377–445 (1908). [CrossRef]

**18. **H. Du, “Mie-scattering calculation,” Appl. Opt. **43**,1951–1956 (2004). [CrossRef] [PubMed]

**19. **C. Ropers, D.J. Park, G. Stibenz, G. Steinmeyer, J. Kim, D.S. Kim, and C. Lienau, “Femtosecond Light Transmission and Subradiant Damping in Plasmonic Crystals,” Phys. Rev. Lett. **94**,113901–4 (2005). [CrossRef] [PubMed]

**20. **T.A. Davis and I.S. Duff, “A combined unifrontal multifrontal method for unsymmetric sparse matrices,” ACM T. Math Software **25**,1–20 (1999). [CrossRef]