{'fraction_non_alphanumeric': 0.032359389493211904, 'fraction_numerical': 0.012142676448267139, 'mean_word_length': 4.346821460775473, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 11, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Improving accuracy of the numerical solution of Maxwell's equations by processing edge singularities of the electromagnetic field

In this paper we present a methodology for increasing the accuracy and accelerating the convergence of numerical methods for solution of Maxwell's equations in the frequency domain by taking into account the be-havior of the electromagnetic field near the geometric edges of wedge-shaped structures. Several algorithms for incorporating treatment of singularities into methods for solving Maxwell's equations in two-dimensional structures by the examples of the analytical modal method and the spectral element method are discussed. In test calculations, for which we use diffraction gratings, the significant accuracy improvement and convergence ac-celeration were demonstrated. In the considered cases of spectral methods an enhancement of convergence from algebraic to exponential or close to exponential is observed. Diffraction efficiencies of the gratings, for which the conventional methods fail to converge due to the special values of permittivities, were calculated.


Introduction
It is well known that the electromagnetic field near geometric edges, both conductors and dielectrics, may possess singularities [1]. This leads to a slow convergence and sometimes even to the inability to obtain a solution in the seemingly simplest cases by conventional methods. For example, not so long ago the authors of [2] found that the calculation of the scattering of a TM polarized electromagnetic field on a lamellar diffraction grating ( Fig. 1a) with certain values of the permittivities is not possible even by means of such proven methods as the Fourier Modal Method (FMM, RCWA) [3] and the Analytical Modal Method (AMM) [4,5,6]. It was shown in [7] that the same problem arises under certain conditions in the case of triangular diffraction gratings and generally in the case of any objects where external electromagnetic radiation is scattered at the intersection of faces with certain values of dielectric permittivity and with a certain value of the refractive index of the external environment. In [7,8] the reason of this was investigated and the conditions under which conventional calculation methods cease to converge were deduced. The singularities primarily occur in the components of the electromagnetic field lying in a plane perpendicular to the edge of wedge [1]. In particular, for TM polarization the components of the electric field near the edge of diffraction gratings asymptotically change proportionally to 1 ( , ) ρ ϕ , the constants a, the function ( ) b ϕ and the value of 1 τ characterizing the singularity depend on the grating parameters and the refractive index of the external medium.
In number of simple cases the value of 1 τ can be obtained analytically [9,10], in the rest it is necessary to solve the transcendental equation numerically. In general, 1 τ is a complex value and τ was considered in [2]. Hereinafter, by the term 0 + we denote an infinitesimal positive value. Some modifications of FMM were proposed in [2,11] to solve the problem in the case of lamellar gratings. However, the authors of [8] expressed doubt about the reliability of the obtained results; besides, that did not solve the problem in the whole. The authors of [8] also noted that they do not know the methods currently in use that could give the correct result for 1 Re( ) 0 τ = + . Of course, the case under consideration is specific, however, very often even in ordinary calculations it is necessary to obtain greater speed and accuracy than existing methods allow, which as a rule do not take into account edge singularities of the field explicitly.
To date, various approaches have been developed for processing edge singularities, including: Adaptive Spatial Resolution (ASR) in modal methods [12,13], mesh refinement toward the singularities and various adaptive strategies in the Finite Element Method [14,15], decomposition of the solution into a regular part and a singular part as in the Singular Field Method [16] and in the Singular Complement Method [17,18], the Weighted Regularization Method [19] and the method using of Perfectly Matched Layers at corners [20]. The simplest strategies, such as ASR, do not use additional information about the behavior of the field near the singularity. While more advanced ones, such as the Singular Complement Method, use the explicit knowledge of the expression of the singularities near the edge.
It must be said that the behavior of the electromagnetic field near the edge of wedge as in Fig.  1b) and near the edge of wedge assembly in Fig. 1a) and Fig. 1c) was studied quite a long time ago [21,22], since understanding the law of the field behavior, we can take it into account in the numerical calculation. For example, the convergence of spectral methods which are based on the expansion of the solution in smooth functions such as orthogonal polynomials or trigonometric functions, to which FMM and AMM belong, directly depends on the smoothness of the solution [23]. If in the region under consideration the solution and all its derivatives are continuous, a very fast, so-called exponential or spectral convergence is realized, when the higher expansion coefficients decrease exponentially. If the solution has singularities in the expansion domain, then the convergence will be, at best, algebraic, when the coefficients decrease according to a power law. If the solution in the region near the singularity is known, one can use the expansion in smooth functions outside this region, and then join these two solutions at the boundary, thus very fast convergence can be achieved.
In the original work [24] the field near the singularity was obtained in the form of a power series in ρ , so-called Meixner series, with coefficients calculated from a recurrent system of differential equations. It was later clarified [25] that in the case of angles that are multiples of / 2 π this series also contains the logarithms of ρ , and that the logarithms appear in a more general case when angles are rational numbers multiple of π [25,26].
Based on the results obtained earlier by other authors, we somewhat reformulate the method of calculating the field near the singularity for more convenient integration with solution methods used in the rest part of the domain. We implement the transition from the static case, where the problem re-duces to solving the eigenvalue problem similar to that arising in modal methods, to the case of finite wavelength. As examples of practical application the solution near the edges is added to several wellknown methods. In the first case it is the AMM, where a system of linear equations for the coefficients of the transmitted and the reflected modes is modified as a result of taking into account the field singularity. In the second case it is the spectral element method, where the region containing the singularity acts as one of the elements. In addition, an algorithm for constructing the matrix of Dirichlet-to-Neumann map operator (DtN) is presented for the region containing the singularity. This allows adding singularity processing to any multi-domain method for solving Maxwell's equations in the twodimensional case when the DtN algorithm can be applied to join computational domains [27]. In particular, we briefly touch on the boundary integral method [28] by which a verification of reference values is carried out for the examples presented in this work. In all cases considered, the singularity handling allows to significantly increase the accuracy and the computation speed, as well as achieve exponential or close to exponential convergence. Additionally, we performed calculations of the electromagnetic field in the cases 1 Re( ) 0 τ = + analyzed in [7,8], for which it was not possible to obtain the result by conventional methods and the author knows only one method described in Ref. [20] which have allowed to tackle such problems.

Solution near the edge of a wedge assembly in polar coordinates
The general solution of Maxwell's equations for nonconical diffraction in the two-dimensional (2-D) case can be represented as a linear combination of two fundamental polarizations: the transverse electric polarization (TE) and the transverse magnetic (TM). In this paper, we shall focus on the case of TM polarization as on a more complex one; in the case of TE polarization, everything can be done in a similar way. For TM polarization only the z H component of magnetic field is not zero and the equation in the frequency domain for this component in polar coordinates is: where ( ) ε ϕ is the relative dielectric constant, 0 2 / k π λ = , λ denotes the wavelength of the incident radiation. Hereafter, for simplicity we set the relative magnetic permeability to unity. We shall solve Eq. (1) in a circle 0 ρ ρ ≤ , see Fig. 1). For convenience we introduce a new variable: 0 / r ρ ρ = and equation (1) takes the form: where ( ) 0 2 / γ π ρ λ = is the dimensionless constant much less than unity.
We shall seek the field z H in the form of a series:  (2).
It is solved by the separation of variables: where ( ) j ϕ Φ are the eigenfunctions of the one-dimensional Helmholtz equation: corresponding to the eigenvalues 2 j τ . For definiteness, we shall chose the sign of j τ so that Re( ) 0 j τ ≥ .
Equation (6) is completely analogous to the eigenvalue problem that arises in modal methods, for example, in case of a lamellar grating with replacing the x coordinate by the azimuth angle ϕ and the wavelength λ → ∞ 0 ( 0) k = , see the next paragraph. The solutions of Eq. (6) in some simple cases can be found analytically (see Appendix 1). The main properties of the eigenfunctions ( ) , including their linear independence and completeness, which is important for representing through them an arbitrary condition on the boundary 1 r = , as well as the orthogonality with weight 1 ( ) ε ϕ − to the adjoint eigenfunctions ( ) j ϕ + Φ corresponding to the differential operator adjoint to the operator in Eq. (6) are given in [6], and the general theory of such non-self-adjoint boundary eigenvalue problems is presented, for example, in [29].
Substituting Eq. (5) into Eq. (4) and taking into account Eq. (6), we obtain the equation for 0 ( ) j R r : with solution: For Re( ) 0 j τ > the term with a negative sign in exponent is rejected due to physical reason (the requirement of a finite energy in a finite volume) 1 . We set the coefficients 1 j a = so that on the boundary 1 r = the functions 0 ( ) j R r are equal to unity. In case of 0 γ > we shall seek the functions ( , ) j G r ϕ in the form: with the condition that upon passing to the static case 0 where matrix elements , and we take into account the orthogonality of ( ) The functions , ( ) k j R r can be found by the iteration method, taking the static case 0 ( ) j R r as the zeroth-order approximation (see Appendix 2.1). The solution is obtained as a series in even powers of γ : where , , ( )  (14) absolutely converges for any finite values of γ when the wedge angles are the product of π and rational numbers, and we briefly discuss the case of angles that are multiples of an irrational factor of π ; there are also formulas for practical calculations of the coefficients of the series, see Eq.(110)-(112). Since we know the behavior of the functions , ( ) k j R r for 0 r → , this functions can also be found by a direct solution of the Cauchy problem for a system of ordinary differential equations (Appendix 2.3).
Summarizing all the above, the solution near the edge of a wedge assembly inside the circle 0 ρ can be written as follows: where j τ are the square roots of the eigenvalues such that Re( ) 0 j τ ≥ , and ( ) ϕ Φ are corresponding eigenfunctions of the problem (6) arising in a static case ( λ = ∞ ), , 0 ( / ) k j R ρ ρ are the solutions of the system of ordinary differential equations (12), and the coefficients j h are found from the boundary conditions at 0 ρ ρ = .

Incorporating the solution near the edge of a wedge assembly into the Analytical Modal Method
In this section we shall apply the above-described technique for calculating the electromagnetic field near edges to accelerate the convergence of the Analytical Modal Method which is the base representative of the class of modal methods. Within these methods the simulation domain, for example, a diffraction grating, is divided into parallel layers, in each of which the permittivity changes only in the plane of the layer and is assumed to be constant along the perpendicular direction allowing the separation of spatial variables, see Fig.1a). Within each layer the eigenmodes of the electromagnetic field are calculated and the general solution is then expressed by means of an eigenmode expansion. The expansion coefficients can be found by applying the proper boundary conditions. In the case of TM polariza-tion the magnetic field within each l-th layer (l=0..L, where l=0 corresponds to the superstrate layer and l=L to the substrate) with relative permittivity ( ) l x ε satisfies the following equation: where p is the period of grating, n 0 is the refractive index of the superstrate and θ is the incidence angle of external radiation in the form of a plane wave ( Fig. 1): After separating the variables in Eq. (16), the solution within each l-th layer can be written as [4,6]: where , The sign of square root in (19) is chosen so that its imaginary part would be positive, whereas for real positive 2 , l n κ the positive square root should be used. Thus, the directions of propagation and decaying of modes coincide. Various methods for solving Eq. (20) can be found in [4,6,30,31] and the same methods are applicable for solving Eq. (6).
In the edges regions bounded by circles of radius 0 ρ and marked in Fig. 1a) by a dashed line the solution is written in the form (15). In numerical calculations, we restrict the sum (15) over j to the first w N terms. The number of terms in the sum over k is chosen to be larger than w N , so that with their further increase the final calculation results do not change within the required accuracy. For simplicity the numbers w N will be taken the same for all wedge assemblies w. The equations for the coefficients c in Eq. (19) for each layer and the coefficients h in Eq. (15) for each region near the edge are obtained by equating the values of the field and its normal derivative multiplied by 1 ε − on both sides of the interface Γ indicated in Fig. 1a) by the dashed line: with a subsequent integration of these equations over a system of test functions. . More details on the preparation of equations for the coefficients can be found in Appendix 3.
The radius 0 ρ of the circle around the edge of a wedge assembly is chosen so that the ratio of the number of functions (9) to the length of the projection of the circle boundary onto the x axis, see Fig. 2b), was approximately equal to the ratio of the number of modes to the period. The radius in this case is 0 / (4 ) At lower values of 0 ρ the convergence slows down; at large values, in addition, the condition number of the matrix of the resulting system of linear equations grows rapidly. The condition number increases for the same reason as in the Rayleigh method [32], due to the fact that exponentially decaying modes in Eq. (19), when integrated over the radius of the circle, give a very small contribution, which, for a finite length of the mantissa of numbers, is practically not distinguishable for different modes against the background of the remaining coefficients of the linear system.  Table 1. Parameters of the lamellar diffraction gratings in Fig. 1a) used as examples. The TM-polarized electromagnetic wave is incident at an angle θ .
As practical examples we calculated the diffraction efficiencies of several lamellar gratings in Fig. 1a). Table 1 shows the parameters of these gratings and references to some works in which they were previously used for such tests. In this paper we do not touch on hypersingularity [37]   in the PFMM method is selected from the condition of best convergence. The value 0 0.8484816789046437 R = , calculated using multiple-precision arithmetic and rounded to 16 significant digits, was used as a reference. The reference values used in the work were verified by several methods, including the boundary integral method with incorporated edge singularity processing; for more detail, see the Section 4 of the article. From Fig. 3 it can be seen that the FMM, PFMM methods and the regular version of AMM ( 0 w N = ) significantly lose in convergence rate to the modified AMM method, which explicitly takes into account edge singularities. In addition, in the case of a proportional increasing the number of functions w N in expansion (15) simultaneously with the number of modes M the modified AMM method shows exponential or close to exponential convergence (Fig. 3a), while the rest of the presented spectral methods are limited by algebraic convergence due to the presence of edge singularities.  Table 1 versus the number of modes M calculated by various modal methods; b) the results of calculation using the modified AMM method with a fixed number of functions N w used for representing the solution near the edges.
In the case of a fixed number of functions w N in the expansion (15) of the solution near the edges the convergence of the method is algebraic, Fig. 3b). Besides, with increasing w N the rate of convergence grows significantly. It can be noted that if earlier, in order to achieve relative accuracy 12 10 − by the usual AMM method M= 5

10
× modes were required in this example [31], now, to obtain the same accuracy using the modified AMM method, only a few hundred modes are enough. Fig. 4 shows the calculation of the first-order transmitted diffraction efficiency 1 T of a dielectric lamellar grating No. 2 of Table 1. As in the case of the metal diffraction grating the processing of edge singularities significantly accelerates convergence. And with a simultaneous proportional increase in the number of functions w N together with the number of modes M, an exponential or close to exponential convergence is also observed in Fig. 4a).  Table 1 versus the number of modes M calculated by various modal methods; b) the results of calculation using the modified AMM method with a fixed number of functions N w used for representing the solution near the edges.
In the following example, Fig. 5, the reflectivity of a lamellar diffraction grating No 3 of Table 1, featured the real part 1 τ of all four edge singularities equal to +0, is calculated. It was previously studied in the papers [2,8,11], which we talked about in the introduction. Conventional methods such as FMM, PFMM, and AMM do not work in this case, and only the results of the modified AMM method are shown in the figures. If the singularity is explicitly taken into account, the value     Table 1, depending on the required calculation accuracy and the value of 1 Im( ) ε . In the area highlighted in blue it is more advantageous to use the PFMM method, in the green area the modified AMM method is preferable.
Until now, we have compared the convergence of the methods depending on the number of modes M . However, in real applications the time of calculation of the result with the required accuracy is primarily important. The calculation of the field near the corners in the form (15) takes a certain amount of time and, on the whole, the modified AMM method, which explicitly takes into account the edges singularities, works slower than the usual AMM method with the same M . The number of operations when calculating by formula (15) is proportional to 4 w N and in the case when w N M ∼ the total number of operations in the modified AMM method grows as 4 M . While in most modal methods the number of operations is proportional to 3 M , moreover, in the iterative implementation of the AMM method it is proportional to 2 M [31]. However, the convergence of the conventional methods is only algebraic, that is, the calculation errors is proportional to q M − with a certain constant q . At the same time, the convergence of the modified AMM method is exponential or close to exponential, and the calculation error is proportional to e c M − ⋅ with some constant c . Therefore, with an increase in the result accuracy requirements, the modified method will ultimately win in terms of time consumption. Also, with decreasing of the values of 1 Re( ) τ , conventional methods need an exponential increase in the number of modes to achieve the required accuracy, whereas the convergence of the modified method practically does not change. Thus, conventional methods, which do not take into account the edges singularities explicitly, can gain in computation time only if the required accuracy is moderate and the value of 1 Re( ) τ is not too small. As an example, we provide a diagram approximately showing when one method or another is preferable in the case of diffraction grating No. 3 of Table 1, depending on the required calculation accuracy and the value of 1 Im( ) ε . The modified AMM method is compared with the PFMM method in Fig. 7). This test was carried out using only one processor core without any parallelization.

Incorporating the solution near the edges into the Spectral Element Method
In this section, we shall incorporate the solution near the edges into the Spectral Element Method (SEM) [38], which is a universal method and which is very close to the finite element method. In this method the computational domain is divided into subdomains, in each of which the permittivity remains uniform or vary smoothly. The solution in the subdomains is represented in the form of an expansion over known basis functions, for which the orthogonal polynomials are usually used. In a subdomain containing an edge we shall represent the solution in the form of the series (15). The equations for the expansion coefficients in SEM are wrote in such a way that within each subdomain the solution satisfies Maxwell's equations, and at the boundaries with other subdomains it satisfies the corresponding conditions for the electromagnetic field. The coefficients can be calculated all at once by solving one large system of linear equations, which is usually done by some iterative method or, by implementing the Dirichlet-to-Neumann mapping (DtN) technique, where firstly the boundary conditions for the subdomains are found, and after that the Dirichlet problem is solved separately for each subdomain [27,39,40].
Speaking about the spectral element method, we adhere to the terminology of the book [23], referring to methods that use the spectral approach inside the subdomains, not only limiting by the Galerkin method in a weak form, as in the original proposal [38]. A formalism based on the application of the Dirichlet-to-Neumann (DtN) map, which we shall use, was previously discussed by us in [41], here we only briefly outline the main points.
In the case of TM polarization the magnetic field in each subdomain with uniform permittivity ε satisfies the following equation: It can be solved both by the Galerkin method and by the pseudo-spectral method [23,39,40]. If a subdomain has not the rectangular form, then mapping into a rectangular area is applied by introducing new coordinates. For joining solutions in neighboring regions, it is convenient to use the Dirichlet-to-Neumann map [27,39,40]. To obtain the DtN map operator matrix for a subdomain, on the boundary Γ of this subdomain a set of N linearly independent boundary conditions at N points of the Gauss or Gauss-Lobatto quadrature is defined and, by solving the Dirichlet problem for each of them, the corresponding normal derivatives at these points of the boundary are calculated. After that we compose two square matrices N N × : in the first matrix F there will be N boundary conditions as columns, and in the second matrix ∂ n F there will be columns of the corresponding normal derivatives at the boundary. Since the boundary condition vectors are chosen linearly independent, an arbitrary boundary condition Γ f defined at these points can be represented as their linear combination with some coefficients с, and the corresponding normal derivative Γ ∂ n f can be easily calculated: where D is the DtN map operator matrix. Boundary conditions, as well as normal derivatives, can be specified not only in the form of values at N points, but also in the form of N expansion coefficients for some fixed set of basis functions at the boundary, which is convenient if we are calculating a solution inside the subdomain using the Galerkin method. In this case, we shall say that DtN map acts in the space of these basis functions. The transition from one representation to another, as well as the calculation of the DtN map matrix for the subdomain containing the edge in various representations is given in Appendix 4.
Knowing the DtN map operator matrices 1 D and 2 D for two subdomains with a shared border ′ Γ : Conversely, if the conditions on the boundaries 1 2 , Γ Γ are known, the conditions on the common boundary ′ Γ will be: If we want to connect two adjacent boundaries 2 Γ and 3 Γ of the region characterized by a matrix of DtN map operator: with conditions on the boundaries: where , α β are diagonal matrices as, for example, in the case of quasiperiodic boundary conditions (17), then, substituting Eq. (28) into Eq. (27), we obtain the DtN map for the remaining part of the boundary 1 Γ : And vice versa, after the conditions on the boundary 1 Γ are found, the conditions on the boundaries 2 3 , Γ Γ will be: Thus, a DtN map operator matrix of the entire domain can be constructed from individual elements. After such a matrix is calculated, the unknown coefficients 1 − c of the reflected wave in superstrate and 2 + c of the transmitted wave in substrate can be found using expansion (19) of the field in the superstrate and in the substrate, and the definition of the DtN map at the upper and lower boundaries of the domain: where the square matrices l Ψ consist of the columns , are known and specify the zero order incident wave in the superstrate and the upward going wave in the substrate, respectively. Instead of equating the field at grid points, the Galerkin method can be implemented by equating field integrals with test functions on both sides of the boundary, for more details see [41]. After the magnetic field at the boundary of the entire domain is obtained, we can, using Eq. (26) and Eq. (30), and acting in the opposite order, find the boundary conditions for each of the subdomains. Then, the field inside the subdomains can be calculated by solving the Dirichlet problem. In particular, the calculation of the coefficients j h in the expansion (15) of the magnetic field in the region near the edge is given in Appendix 4. It may be noticed that the DtN mapping approach is very similar to the S-matrix technique used in modal methods [42], with the difference that instead of reflection and transition coefficients, the field values at the subdomains boundaries are searched.
It should be noted that in the considered approach the DtN map matrices of individual subdomains can be calculated by using an arbitrary method, not necessarily by the spectral method. In the case of a uniform dielectric constant inside each element, it is very convenient to calculate the DtN map matrix by using the boundary integral method [28,36]. Since in this method the calculated subdomain does not need to be mapped into a rectangle, we get considerable freedom in choosing the partition into subdomains. In particular, we can represent regions of uniform permittivity as these subdomains. The subdomains containing the edges, as before, are calculated separately, and the whole technique of combining the subdomains remains the same. The boundary integral method scales much better with the number of collocation points increasing, and we employed its multiple-precision arithmetic implementation with incorporated edge singularity processing to verify the reference values used in the article.
As practical calculations shown, the value of the radius 0 ρ of the region containing the edge, where the solution is represented in the form of series (15), can be chosen in a first approximation so that the ratio of the number of functions w N to the length of the boundary of the region is equal to the average density of grid points per period, in accordance with the expression: where x N is the number of grid points along the x axis over the grating period p. If the high accuracy result is required, then for faster convergence it is convenient to use the formula: where the coefficients i с can be approximately chosen from the condition of the best convergence of the method to the result of intermediate accuracy already obtained. In the work we used this formula with two coefficients i с in the series. As examples, several triangular diffraction gratings Figs. 1b), the parameters of which are listed in Table 2, were used. Gratings with the same or similar characteristics were considered earlier in [7,8].  Table 2. Parameters of the triangular diffraction gratings Fig. 1b) used as examples. A TM-polarized electromagnetic wave is incident at an angle θ . Fig. 8a) shows the relative error of the zero-order reflected diffraction efficiency 0 R of a grating No. 1 of Table 2 versus the number of grid points x N along the x axis over the grating period p, calculated by the SEM with various numbers of functions w N in the subdomains containing edges. Previously this example was considered in [7], where the absence of the C-method (or Chandezon method [43]) convergence was observed due to the purely imaginary value of 1 τ . Note that in the case of a triangular grating Fig. 1b), the value of 1 τ for the edge at the top vertex of the profile is the same as at the bottom; see the remark following the Eq. (41) in Appendix 1. It can be seen that the situation is the same as in the case of the analytical modal method: when the number of functions w N is constant, the convergence of the spectral method is algebraic, and is exponentially fast or close to it if w N is increasing proportionally with the number of nodes x N . In Fig. 8b) the spatial distribution of the squared magnitude of the magnetic field  It is seen that near the grating boundary, there are regions in which the magnetic field is much larger than the amplitude of the incident plane wave Eq. (18), which is equal to unity. Fig. 9a) shows the relative error of the reflected minus first-order diffraction efficiency 1 R − of a triangular grating No. 2 of Table 2 versus the number of grid points x N , calculated by the SEM with various numbers of functions N w in the subdomains containing edges. In this example, as in the previous one 1 Re( ) 0 τ = + , therefore, the conventional numerical methods also do not work here, which was pointed out in [7].  1b). So that the angles of the subdomains are as close to 90 0 as possible, this reduces the condition number of the linear system matrix when we solve the Dirichlet problem inside these subdomains. The magnitude of the magnetic field reaches its maximum near the lower angle of the grating, in the region containing the edge that is shown beside.
In Fig. 10a) the convergence of several methods was compared using the example of a triangular grating No. 3 of Table 2. We examined: a) the conventional spectral element method without increasing the accuracy of the solution in the edge regions; b) the parametric C-method with adaptive transformation of coordinates near edges (or Adaptive Spatial Resolution (ASR) as it is called in the original paper [13]) and a transformation parameter   Table 2 versus the number of grid points x N over the grating period p for the spectral element method (SEM) and versus the number of modes M for the parametric C-method with Adaptive Spatial Resolution (ASR); b) the results of calculation by the spectral element method with a fixed number of functions N w used for representing the solution near the edges.
The spectral element method employing expansion (15) demonstrates the best convergence. In Fig.  10b) the convergence of this method is shown for a different number of functions N w in the subdomains containing the edges.  how it was done in the last example of the previous section, we get a physically incorrect result: 0 0 R < . The polynomial extrapolation of 0 R using the same points gives a relative error of about 30%. To improve the result, it is necessary to shift the extrapolation points closer to Im( ) 0 s ε = , which requires an exponential increase of the computation time in the case of conventional methods. Fig. 11b) demonstrates how the convergence rate of the parametric C-method diminishes, with a decrease in the imaginary part of s ε and a corresponding decrease in the real part of 1 τ .
For 1 Re( ) 0 τ = + the parametric C-method, like the conventional SEM, does not converge to a certain value. The reason for the divergence, consisting in the behavior of the field near the edges, has already been discussed in the papers [7,8]. Having the ability to calculate the field near the edges, we shall illustrate this using the example of the triangular grating No. 3 of table 2. In Fig. 12 Table 2 for two values of the dielectric constant of the grating material. The thin white lines indicate the boundaries of the subdomains, and the thick one marks the grating profile.
In the presented resolution the change in the field in both pictures seems quite smooth, and although in the second case its maximum value is several times larger, the behavior of the 2 | ( , ) | H x y distribution cannot say why the conventional numerical methods work in one case and not in the other. However, if we sequentially increase the scale of the region near the upper edge of Fig. 12b), we obtain a picture of the distribution of the squared magnitude of the magnetic field shown in Fig. 13.  Fig. 13 is the same as in the Fig.   12. In each subsequent image of Fig. 13 the radius decreases by about 1/2 d =23.3 times, the scale on the right is the same for all images. Despite of the reduction in scale, the field continues to change without decreasing the amplitude, without reaching a constant value. And if we wanted to represent such a distribution in the form of a complex Fourier series as in the C-method, or by using an orthogonal polynomials expansion as in the spectral element method, we would need more and more functions to take into account field changes.
To approximately represent the magnetic field near the top of the profile for small values of ρ , we restrict ourselves to the first terms of series (15) containing only 0 τ and 1 τ : where we consider that 0 0 τ = and 0 ( ) 1 ϕ Φ = .
where ( ) a ϕ and ( ) b ϕ do not depend on ρ . The squared magnitude of the magnetic field in this case will oscillate with an increasing frequency when ρ goes to zero. And for the radii 1 2 1 / exp( 2 / ) d ρ ρ π β = − ≡ the picture will begin to repeat, as we saw in Fig. 13. The components of the electric field will oscillate with increasing amplitude. Fig. 15 shows the squared magnitude of the magnetic field depending on the distance to the apex of the triangular grating (see Fig. 12) when moving along the profile from the left side to the top. The distance is plotted on a logarithmic scale. The dashed line shows the asymptotic approximation calculated by Eq. (34). In the case of 1.5 1.5 s i ε + = − , see Fig. 15a), when approaching the peak ( 0 ρ → ) the oscillations decay quickly, and the magnetic field converges to a constant. In this case the field can be represented with any predetermined accuracy by using a finite number of basis functions, and the usual calculation methods work well here. However, the smaller the value of 1 1 Re( ) α τ ≡ in Eq. (34), the slower the oscillations damping will be, and the more functions will be needed to represent the solution at a given accuracy. This behavior is seen in Fig. 11b). For 1.5 0 s i ε = − + the oscillations, in the case when the quantity ρ decreases, persist with a constant amplitude and with an increasing frequency, therefore, it was not possible to obtain the result by conventional methods [7,8]. Nevertheless, the explicit calculation of the field near the edges, as shown in the present paper, copes with this problem.

Conclusion
In this work we studied an approach which improves the accuracy and accelerates the convergence of numerical methods for solving the Maxwell's equations by means of taking into consideration the singularities of the electromagnetic field that arise near the geometric edges of scattering objects. Several algorithms to incorporate the singularities processing into methods for solving the Maxwell's equations in two-dimensional structures were analyzed. As implementation examples we used the analytical modal method and the spectral element method. In the first case, the resulting system of equations for the coefficients of the eigenmodes underwent changes. In the second case, for each region containing a singularity a Dirichlet-to-Neumann map operator was constructed and subsequently these regions were processed, like all the others.
In test examples, for which we used the calculation of the lamellar and the triangular diffraction gratings characteristics, a significant increase in accuracy and an acceleration of convergence was achieved.
In the examples considered, the possibility to increase the type of convergence of spectral methods from algebraic to exponential or close to exponential was demonstrated. Diffraction efficiencies of the gratings, for which the conventional methods fail to converge due to the special values of permittivities, were calculated. In this work we confined ourselves to the two-dimensional case only; however, the method presented in the article can be extended to the three-dimensional case as well.
In conclusion we note, that the approach presented in the fourth paragraph of the paper, which employs the Dirichlet-to-Neumann map for joining subdomains, is universal and allows to easily em-bed the processing of singularities in various 2D multi-domain methods, regardless of what the numerical algorithm for solving Maxwell's equations is used inside individual domains. In addition to the spectral methods analyzed in the article, it can be the boundary integral method, the finite-difference method, the finite-element method, and etc.

some simple cases
The examples considered here have already been investigated previously in [9,10,24]. For completeness, we present the main results for simple cases, writing out the whole series of roots of the eigenvalues , Re( ) 0 k k τ τ ≥ , taking into consideration the degree of degeneracy, which is important for the correct writing of system (12). Note that when we write the final formulae for k τ in this appendix, we do not sort k τ by the value of the real part. Let us at first solve the problem Eq. (6) in the case of a homogeneous infinite wedge with an angle ϑ made of a material with a permittivity 1 ε in a medium with a permittivity 2 ε , Fig. 16a). We position the x axis on the axis of symmetry, and we shall seek for the function ( ) k ϕ Φ in the upper half-plane 0 ϕ π ≤ ≤ , in the lower half it will be either antisymmetric or symmetric. Accordingly, at points 0, ϕ π = either the function itself or its derivative must vanish. In the first case, the solution can be written as: and in the second: Similarly for the symmetric case: 1 2 cos / 2 sin ( / 2 ) cos ( / 2 ) sin / 2 0 It is convenient to rewrite equations (39), (40) in the form: where the "+" sign in front of δ corresponds to the antisymmetric case and "-" sign to the symmetric one. Note that replacing 1 2 ε ε ↔ is identical to changing the sign " ± " before δ in Eq. (41) Such a notation explicitly takes into account that the eigenvalues corresponding to 2, 4,6,...

τ =
are doubly degenerate. Note that the middle terms in the brackets of expression (44) are not sorted by the value of the real part, as well as in the other formulae in this appendix, and to determine the 1 τ characterizing the singularity it is necessary to choose one of them whose real part is minimal. In numerical calculation only the first w N terms of the series { } k τ sorted in ascending order of the real parts of k τ are used. In general case of arbitrary ϑ when k τ is calculated numerically, for example, by method represented in Ref. [31], the symmetry of the eigenvectors corresponding to nondegenerate k τ can be readily identified by substituting k where we have divided the equations for derivatives by the factor i τ ⋅ , since we consider the constant solution at 0 τ = already known. In order that the linear system (50) has a nonzero solution, it is necessary that its determinant must be equal to zero, from this condition the values of k τ can be found.
Or, after simple conversions: Zero solutions 0 x = have no physical meaning, two solutions with 1 x = correspond to 3,4 2 τ = , and solutions corresponding to the last bracket vanishing: are:

Appendix 2. Calculation of radial functions R 2.1 Iterative calculation of radial functions R
We have to solve the system of differential equations (12) in the domain 0 1 r ≤ ≤ , to obtain a set of functions , ( ) k j R r satisfying condition (10). In the numerical solution we restrict ourselves to the first functions , ( ) − corresponding to the roots k τ , sorted in ascending order of the real part. We shall solve the system of equations (12) iteratively, finding the new l+1 approximation [ 1] , l k j R + by substituting of the previous [ ] , l k j R into the right side of the system: 1 The property of the continuous dependence of the roots k τ on the dielectric constant is used in some numerical implementations of the analytical modal method to calculate the values of k τ , see, for example, [30].

( )
As a zeroth approximation we take the static solution: only in the initial approximation (59). Solution of Eq. (58) is the sum of the general solution of the corresponding homogeneous system and the particular solution of the inhomogeneous one. As a solution of the homogeneous system that satisfies the finiteness condition at zero and Eq. (10), we take the static solution (59) at each iteration. At the first iteration the solution of equation (58) is In Eq. (60) we obtained the sum of the zeroth approximation with a term containing the powers and logarithms of r multiplied by 2 γ . If we continue this process the result will be a series in even powers of γ . Besides, where , , ( ) At the next iterations we have to solve equations containing the right-hand side of the form: where the degree of the logarithm β is an integer, A particular solution to the equation with such a right-hand side can be calculated using formal integral representation: where in the third line for 2 k τ α ≠ + the lower limit of the first integral must be assigned as 0 σ = if Re( ) Re( ) 2 k τ α < + and as σ = +∞ otherwise. Integrals with a logarithm are calculated using the recurrence formula: So, Eq. (63) can be rewritten in expanded form: Thus, the solution will look like the series: where , ( ) option is realized when the angles of the wedges are the product of π and rational numbers. As practical calculations show, for better accuracy the number N Φ should be taken much greater than w N , and on finishing the iterative process , ( ) k j R r can be truncated to a square matrix if necessary. The required number of functions N Φ can be determined by taking the value several times larger than w N as the initial approximation, and then double N Φ until the changes in the final calculation results become less than the specified error.

Convergence of the series ( ) k, j R r
In this appendix we shall show that the series (66) converges absolutely for any finite values of γ , when the domain under consideration consists of wedges with angles multiples of / m π , where m is a positive integer. We also briefly discuss the case of angles that are multiples of an irrational factor of π .
If the iterative process (61) is briefly written in the form: , then the series (66) will look like: To prove the convergence of this series, below we will construct a majorizing series and show its convergence.
To write down analytical formulas for the terms of series (67), let us trace their evolution. During the first 1 where we have introduced the notation ( ) 2 The first upper index in [ ,0] , l k j A corresponds to the iteration number and the second to power of the logarithm by which this coefficient will be multiplied. Note that ( ) 2  The finite length of the period l ∆ is a fundamental feature of the case of angles that are multiples of the rational factor of π . This allows us to determine finite j d Eq. (72) and prove the convergence of the series (66) simply enough. Also in this case we can choose a continuous function ( ) j D l that grows with increasing l no faster than logarithmic function, so that: where 0 j D > is some constant. This statement is easily proved by dividing (73) into the inner and outer sums as in Eq. (70) and replacing the inner sum over n by a majorizing integral with appropriately chosen limits.
Note that if the condition ( ) for all other equations of the infinite system (58), the solution remains in the form (68) with the coefficients: We specially added a factor j d to the denominator of the expression for the coefficient The solution to Eq. (77) for The sum in the expression (79) for the coefficient 1 [ 1,0] , and the product in the case of coefficient 1 [ 1,1] , where for brevity we introduced the notation for some 2 k is satisfied at the next iterations, the logarithm of the first degree will be present on the right side of each line of Eq. (58), and for iteration 1 2 L l L < < the solution will be written as: For all other lines of the solution everything will be as before, in accordance with Eq. (82 If we further continue this process, then in each iteration L β , when ( ) is satisfied for some k β , in the solution for the k β -th row of system (58) (or in more rows, if there is a degeneracy), the logarithm to power β is added and the factor of this logarithm decreases with increasing β : The factorial in the denominator of Eq. (87) arises due to the fact that each logarithm, when it appears in Eq.
where the coefficients [ , ] , L q k j A β are defined recurrently as nested series: ,   , l q k j A . If with an increase in the iteration number l they change according to the power law: [ , ] , l q l k j A α ∼ , where α is a constant, then the parameter ( ) 0 2 / γ π ρ λ = characterizing the radius of convergence 0 ρ is bounded above by some value where we denote The absolute value of [ , ] , For definiteness, let us analyze the case, when in each iteration l the equality ( )

L q k j
A β β will be [ , ] , Passing to the next ( 1) L β + -th iteration, let us estimate the sums over s in Eq. (89), considering Eq. (94)-(96): [ , ] [ , ] 2 0, , , where η is the order of degeneracy for Note that the term with the factor ( )( In the (  2) L β + -th iteration the sums over s in case of q β ≤ can be estimated as:   [44]. And we assume that the order of degeneracy η is the same for