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 (ε = ε ₀ ε_r) which is a function of the spatial coordinates (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.

 

Fig. 1. Remeshing of one element e ¹ ₁ into e_k ² for (a) the basic adaptation (k = 1−2, with 3 nodes), (b) the adaptation with the h-method interpolating the error on the PDE solution (H_z) and (c) the improved adaptive scheme. The adaptives scheme (b) and the improved one (c) using h-method produce more than two new elements (k = 1−4, with 5 nodes). Moreover, with the basic adaptation (a) and the a posteriori h-method (b), the error is estimated on the computed PDE solution Hⁱ_k. This contrast to the improved scheme (c) for which S_ϕkⁱ denotes the interpolation error of the solution of interest (ϕ = E or ϕ = P,‥) on the remeshing with respect to the threshold δ_ϕi for the iterative step i.

Download Full Size | PDF

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 ¹ ₁ 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 (δ_E or δ_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:

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 ²+ν ² 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 + j1.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∣²/∣E₀∣² 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.

 

Fig. 2. Geometry of the study for the infinite circular-cylinder along the z-axis illuminated by a p-polarized incident wave of vector k.

Download Full Size | PDF

Table 1. Number of nodes and computational time (in s) for each mesh step for the classical remeshing and the improved adaptive remeshing process.

View Table

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∣²/∣H₀∣² 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).

 

Fig. 3. Normalized intensity of the electric field ∣E∣²/∣E₀∣² in the xy-plane computed by the FEM with the classical adaptive remeshing (a,c,e) and the improved adaptive scheme (b,d,f). The number of nodes are (a) N_nodes = 892, (b) N_nodes = 768, (c) N_nodes = 4818, (d) N_nodes = 5556, (e) N_nodes = 71260 and (f) N_nodes = 31669, respectively.

Download Full Size | PDF

 

Fig. 4. Computational mesh for the classical adaptive remeshing (a,c,e) and the improved adaptive scheme (b,d,f).

Download Full Size | PDF

 

Fig. 5. Last computational mesh for (a) the improved adaptive scheme and (b) the classical adaptive remeshing. The adaptation to the solution of interest (intensity of the electric field) clearly appears with the improved remeshing.

Download Full Size | PDF

 

Fig. 6. Normalized intensity of the magnetic field ∣H∣²/∣H₀∣² in the xy-plane computed by the FEM with classical adaptive remeshing (a,c,e) and improved adaptive scheme (b,d,f). The number of nodes are (a) N_nodes = 892, (b) N_nodes = 768, (c) N_nodes = 4818, (d) N_nodes = 5556, (e) N_nodes = 71260 and (f) N_nodes = 31669, respectively.

Download Full Size | PDF

 

Fig. 7. Evolution of the normalized intensity of the electric and magnetic field as a function of the mesh refinement (a,c) and (b,d) the errors, relatively to Mie computation, as a function of the distance from the center of the nano-object (radius a = 15 nm) along the x-axis for successive mesh refinements with the classical and the improved adaptive remeshing.

Download Full Size | PDF

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_nodes = 670) at the boundaries of the subdomains (nanoparticle, substrate).

 

Fig. 8. Geometry of the study for the infinite square-cylinder along the z-axis.

Download Full Size | PDF

The remeshing procedure consists in computing the H_z component, deducing the E_x and E_y and their associated P_x and P_y 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 N_nodes = 1580, N_nodes = 3371, N_nodes = 13205, N_nodes = 24906 and N_nodes = 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 ∣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.

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.

 

Fig. 9. Normalized amplitude of the Poynting’s vector ∣P∣/∣P₀∣ computed by the FEM with improved adaptive mesh scheme and the mesh view for (a,b) N_nodes = 3371, (c,d) N_nodes = 13205 and (e,f) N_nodes = 78175, respectively.

Download Full Size | PDF

 

Fig. 10. Evolution of the intensity as a function of the distance along (a) the x-axis and (b) the y-axis, for various mesh refinements.

Download Full Size | PDF

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.