Field-based recovery technique for improved adaptive finite element analysis of photonic devices

A new adaptive finite element (FE) technique is proposed to enhance the accuracy of the solution of wave propagation in integrated photonic devices. The suggested recovery technique improves significantly the finite element convergence rate compared to refinement techniques. The method refines the elements based on the field distribution. This recovery technique estimates the error at each element by comparing the linear FE solution and the recovery solution using least-squares method. The numerical aspects and accuracy of the method are investigated using the light propagation in both one-dimensional and two-dimensional photonic devices.


Introduction
Modern photonic devices have complex structures that require efficient numerical techniques capable of modelling such designs.In this regard, the finite element method (FEM) attracts the interest of many researchers (Hoole et al. 1988).The conventional FEM (Obayya et al. 2000a(Obayya et al. , b, 2002) ) is a robust technique for analyzing various optical devices.However, getting accurate results requires fine meshing.This in turn demands high computational resources due to the produced large number of degrees of freedom (DOF) (Fernandez et al. 1993).A possible way to reduce the DOF is to use higher order elements.Although higher order solution reduces the DOF, it produces less sparse matrices, reduces the flexibility of placing nodes, and complicates the problem (Koshiba et al. 2000).Another way to reduce the DOF is to use first order elements with improved accuracy as in smoothed FEM method (Atia et al. 2015).This method relies on smoothing the derivatives across the boundaries between elements to fix numerical discontinuities inherited in the conventional FEM shape functions derivatives.Alternatively, one can reduce the DOF by selectively meshing areas where most of the field exists with small element sizes.However, large element sizes are used for areas where the field vanishes.This technique requires a proper knowledge of the solution where the field is confined.Therefore, adaptive meshing techniques have been proposed (Hoole et al. 1988;Fernandez et al. 1993;Tsuji and Koshiba 2000).The main challenge of any adaptive algorithm is the selection criteria of regions that need more meshing.The node points are distributed based on the errors in these regions.Hence, various error estimators have been developed in order to distribute the node points.An adaptive technique for solving Poisson equation that uses the node perturbation technique as an error estimator was developed (Hoole et al. 1988).Later, a predefined user function to calculate the error over the elements was proposed (Fernandez et al. 1993).Another error estimator was introduced based on power distribution inside the waveguide (Tsuji and Koshiba 2000).However, the field variation in the computational domain reported in Hoole et al. (1988), Fernandez et al. (1993), Morimoto et al. (2021) was not considered.Therefore, the previous adaptive algorithms fail to follow rapid field variations.This is very important for instance to photonic cavities because accurate representation of field decay significantly affects the cavity quality factor.In Schmidt (1993), a beam propagation technique based on new error estimator that compares the linear FEM solution with the local improved quadratic solution was developed.This technique needs an extra effort at each propagation step due to the required iterative calculations to obtain the local improved quadratic solution.
Here, in order to select the refinement elements, the local errors are estimated for each element.The local error is the difference between the linear FEM solution and the recovered solution calculated using the least-squares technique.Such technique is called "recovery-type" estimator (Kikuchi and Saito 2007).Further, the recovery technique is considered as a post processing method to reconstruct and improve the FEM solution (Yan and Zhou 2001).The new adaptive FEM technique is proposed to enhance the numerical accuracy of the solution of the light propagation in integrated photonic devices.The adaptive technique can refine the elements based on the initial field distribution.Our approach can be used in time and frequency domain simulations.However, here we only apply it on frequency domain method (Tsuji and Koshiba 2002).Three different photonic structures are tested and analyzed using the suggested approach.First, the transmission through reflectionless slab waveguide is calculated.Then, 2D-photonic crystal (PhC) cavity is studied which ensures the high accuracy of the quality factor and resonance frequency calculated by the reported algorithm.Further, 2D-PhC multiplexer-demultiplexer is analyzed for wavelengths of 1.3 µm and 1.55 µm where each wavelength is propagated through different PhC channels.It is shown that the reported technique has a node distribution allocated at the field location with fast convergence rate compared to conventional random refinement techniques.

Finite Element Formulation
For a 2D computational domain ( Ω ) where the variation is in xy/ plane, the weak form of the Helmholtz wave equation in the frequency domain can be expressed as with, where k o is the wave number, n is the refractive index, and ψ is the weight function which is expanded using the same basis as ϕ according to Galerkin formulation.The E z and H z are respectively the electric and magnetic field components in z-direction.The s x and s y are the PML parameters in x-and y-directions as shown in Fig. 1 that can be found in Tsuji and Koshiba (2002).
After dividing into first order triangular elements, Eq. 1 is rewritten as where (1) Here, {N} is the FEM shape function vector, and ∑ e extends over all elements.To excite the structure, we define the source plane Γ which divides the analysis region into Ω 1 and Ω 2 subregions.Thus, due to the continuity condition across Γ , Eq. 2 can be rewritten as Tsuji and Koshiba (2002) where ∑ Γ extends over the elements attached toΓ , {N} Γ is the 1D shape functions vector, ∫ l refers to a line integral over Γ, and ϕ 1 and ϕ 2 are the incident fields in the sub-regions , Ω 1 andΩ 2 , respectively.The fields ϕ 1 and ϕ 2 are expressed using the eigenmodes.The lin- ear system of equation presented in Eq. 3 is then solved for the field distribution over the computational domain.

Field-Based Recovery Technique
The proposed adaptive technique is considered as a posteriori error estimate method where the exact solution is replaced, in the expression of the true error, with an approximation recovered from the numerical solution.Figure 2 shows a flow chart of the presented technique.At each iteration, the problem is solved and an error calculation for each element is obtained.Then, the mesh is refined in zones where high errors exist.The error estimation (ε i ) at each element is calculated by the formula where k is the number of elements in the computational domain, ‖ ⋯ ‖ is the L 2 norm, ũ is the linear FEM approximation to the exact solution u and R ũ is the recovered solution of ũ .Once the estimation of the local is calculated at each element, the refinements take place in the elements that satisfy ε i > αε max , where ε max is the maximum local error in the domain and α(0 < α < 1) is a parameter that controls the number of refined elements.When α is close to 1, the refinement takes place only in the elements with an error close to the maximum.On the other hand, when α is close to 0, a large number of elements is refined, even those with a small error.The global error ε g is given by where the adaptation process iterates until ε g is less than a certain tolerance (ϵ) (Kikuchi and Saito 2007).In order to obtain the recovered solution R ũ at a certain mesh node z , a patch Ω z of mesh elements directly attached to z is selected as shown in Fig. 3a.Then, the linear FEM solution is extended to the mid-edges by interpolation (Zhang and Naga 2005;Naga and Zhang 2004).
Next, a quadratic polynomial P z that best fits the FEM solution ũ is constructed, using the least-squares method, at the nodes shown in Fig. 3b.The recovered solution R ũ at a certain interior node z is defined as the fitting polynomial, i.e.
At any boundary node ẑ , the recovered solution is equal to the FEM solution.

Numerical Results
Figure 4a shows a silica slab with refractive index n = 1.44 surrounded by air.The length of the silica slab d = mλ 2n (where m = 1, 2, 3, … and λ is the operating wavelength in air) is chosen to obtain complete transmission with no reflection.During this study, the slab thickness is taken as 1615 nm which corresponds to m = 3.The purpose of this example is to show the convergence of the proposed technique to the exact solution with transmission of unity.An initial mesh is first used then the proposed adaptive technique is applied.Figure 4b, d show the distributed field E y of the travelling wave through the slab wave- guide at wavelengths of λ = 1.55 μm and λ = 0.775 μm .Figure 4c, e illustrate the mesh distribution after applying the proposed technique at the studied wavelengths after applying the mesh refinement technique.It may be noticed that the mesh density follows the electric field pattern produced in the computational domain.This coincides with the idea of the technique where the mesh follows the field unlike the previous work done in this area where the refinement is based on the gradient of the field (Babuvška and Rheinboldt 1978).Figure 4f shows the mesh distribution after applying the proposed refinement technique with DOF = 82,300 at λ = 1.55 µm.The solid green triangles represent the largest element with size 1.25 × 10 −7 μm , while the dark red triangles represent the smallest elements with size 3.88 × 10 −9 μm .It is clear from the figure that the mesh density increases at the peak regions of the electric field shown in Fig. 4b regardless it is a maximum or minimum value.This is due to at the peak regions the gradient of the field component increases that required the mesh to be denser at these regions in order to mimic the rapid changes in the field values.
Figure 5 shows the variation of the absolute error versus the number of degrees of freedom (DOF).The absolute error is computed using The problem is designed to have T exact = 1.0 while T r is the transmission obtained by the suggested technique after mesh refinement process i.e.T r = | | S 21 | | .It may be noticed that the absolute error drops down rapidly when the DOF is between 8 × 10 4 and 1 × 10 5 .Then the error stabilizes afterward.
In order to test the advantages of the proposed method in modeling complex photonic structures, a 5 × 5 PhC cavity shown in Fig. 6 is analyzed.The PhC devices demand numerous numbers of elements for accurate analysis.The selected square lattice PhC cavity (Esquerre et al. 2005) consists of dielectric wires with radius r = 0.2 × a with a lattice constant of a = 0.58652 μm .The dielectric wires are embedded in air and have a refractive index n wire = 3.4 .The central wire is removed in order to form the required cavity.Further, perfectly matched layer boundary condition is applied around the studied PhC to truncate the computational domain.Here, not only the field inside the cavity that dictates the cavity quality factor and normalized frequency but also the field decay and since our approach cares about the error in field no matter how weak it is, the computations of the quality factor is more precise.
In the proposed adaptive FEM, an initial conventional FEM solution is obtained where the computational domain is discretized into 35744 elements.After applying the recov- ery technique, the recovered solution is obtained at each iteration.Based on the initial and the recovered solutions, the local error given by Eq. ( 1) can be estimated.It is required to refine the elements based on a certain ratio α .In this example, the tolerance ϵ is equal to 10 −6 .Figure 7a, b show the initial discretized computational domain and after five itera- tions, respectively.The field distribution inside the cavity is shown in Fig. 6b.It may be seen that the mesh is refined in the subdomains where the field exists without sacrificing the mesh where the field decays.Therefore, this method is called the field-based mesh refinement technique.
Figure 8 shows the normalized frequency and the quality factor variations with the DOF using the adaptive proposed refinement.It is revealed from these figures that the normalized frequency and the quality factor converge to 0.3788 and 177.88 , respectively at DOF of 696,29 .On the other side, Fig. 9 illustrates the studied parameters with the DOF using random refinement.It is evident from these figures that neither the normalized frequency nor the quality factor could converge even for DOF higher than 7 × 10 4 .Table 1 shows the calculated values of the quality factor Q , and resonance frequency f 0 and DOF by the sug- gested algorithm and that obtained in Esquerre et al. (2005).It may be seen from this table that an excellent agreement is achieved.Further, the DOF is reduced by 35% using the proposed adaptive technique.
Next, photonic crystal multiplexer-demultiplexer (MUX/DEMUX) is studied for wavelengths of 1310/1550 nm (Bhatia and Gupta 2014).Figure 10 shows the schematic diagram of the studied PhC demultiplexer that consists of rectangular lattice of silicon pillars embedded in air.The lattice constant is a = 570nm and the rod diameter is d = 0.36a while the refractive index of the rods is n = 3.4 (Yan and Zhou 2001).Two filters are formed through the T-junction by adjusting the pillars diameter.Filter 1 has diameter defect of d 1 = 0.5a which allows wavelength 1310 nm goes upward towards port 3 and filter 2 has  diameter defect of d 2 = 0.72a which allows wavelength 1550 nm goes downward towards port 2. Figure 11 shows mesh refinement generated from the proposed technique at (a) λ = 1.31 μm and (b) λ = 1.55 μm while DOF is fixed to 8.5 × 10 4 .It can be easily noticed that mesh refinement occurs in accordance with propagation of the field inside the device.
The computed electromagnetic field distribution E z for the studied demultiplexer is shown in Fig. 12 for two transmitted signals with wavelength 1.31 μm (Fig. 12a) and 1.55 μm (Fig. 12b).It can be seen that wavelength of 1.31 μm is guided in the upper channel while the lower channel guides the wavelength of 1.55 μm as received by monitor 2 which ensures the high accuracy of the proposed model.
It should be noted that the polynomial-preserving recovery (PPR) technique has been extended to recover solution or continuous gradient from a finite element solution of an arbitrary order in 2D and 3D problems (Naga and Zhang 2005).Therefore, it is believed that the proposed field based recovery technique is applicable to three-dimensional (3D) problems and will be considered for our future work.Typically, the 3D FEM is based on tetrahedron element with four faces which is analogous to the triangle shape in the 2D implementation.Therefore, each FEM node will have three translational degrees of freedom.Consequently, the 3D models are more complex with more variables than 2D models where more computational resources will be needed with long execution time.

Conclusion
A novel adaptive mesh generation is proposed and tested.The technique consists of posterior error estimation calculated based on a previous approximate solution, to guide the generation of a new mesh.This new mesh is built starting from a minimal number of triangular elements, which are then, in several sweeps, repeatedly refined according to the field locations.This method allows optimum use of the available degree of freedom in order to obtain accurate solution.The high accuracy of the reported technique is tested using the transmission through photonic structure, a 2D photonic crystal cavity and 2D photonic crystal multiplexer-demultiplexer.

Fig. 1
Fig. 1 Schematic diagram of FEM discretization of 2D optical waveguide problem

Fig. 2 Fig
Fig. 2 Flowchart of the suggested adaptive FEM

Fig. 4 a
Fig. 4 a Silica slab embedded in Air b, d the distributed field E y of the travelling wave through the studied photonic structure at wavelengths of 1.55 µm and 0.755 µm c, e mesh distribution after applying the proposed technique at the studied wavelengths.f shows the mesh distribution after applying the proposed refinement technique with DOF = 82,300 at λ = 1.55 µm

Fig. 5 Fig. 6 aFig. 7
Fig. 5 The change of the absolute error as a function of the number of degrees of freedom

Fig. 10
Fig. 10 Schematic diagram of the optical MUX /DEMUX based on a PhC waveguide

Fig. 11
Fig. 11 The corresponding mesh refinement generated from the proposed algorithm at a λ = 1.31 μm and b λ = 1.55 μm

Table 1
Esquerre et al. (2005)of the quality factor , resonance frequency 0 and DOF by the suggested algorithm and that obtained inEsquerre et al. (2005)