An Optimized Schwarz Method for the Optical Response Model Discretized by HDG Method

An optimized Schwarz domain decomposition method (DDM) for solving the local optical response model (LORM) is proposed in this paper. We introduce a hybridizable discontinuous Galerkin (HDG) scheme for the discretization of such a model problem based on a triangular mesh of the computational domain. The discretized linear system of the HDG method on each subdomain is solved by a sparse direct solver. The solution of the interface linear system in the domain decomposition framework is accelerated by a Krylov subspace method. We study the spectral radius of the iteration matrix of the Schwarz method for the LORM problems, and thus propose an optimized parameter for the transmission condition, which is different from that for the classical electromagnetic problems. The numerical results show that the proposed method is effective.


Introduction
Nanophotonics or nano-optics [1] is the study of the behavior of light on the nanometer scale, and of the interaction of nanometer-scale objects with light. Moreover, such a discipline is reshaping our worldview in many ways with fascinating (potential) applications such as novel biological detection and new storage media [2]. These applications require fine control of the propagation of light waves. So, it is important to use appropriate mathematical models to describe the behavior of the light-matter interactions [3]. Light waves are regarded as electromagnetic (EM) waves and modeled by Maxwell's equations. Classical or semi-classical models can be employed to model the light-matter interactions, such as the Drude model, the hydrodynamic Drude model [4], and the nonlocal optical response model [5]. In this paper, we consider the classical local Drude model which is a fairly simple yet efficient oscillator model for free electrons in metals that performs well when the size of the considered nanostructure increases beyond ≥50 nm [6][7][8][9].
For a frequency-domain simulation with the Drude model, we are indeed required to solve the time-harmonic Maxwell's equations whose closed-form solutions are not available. Thus, various numerical methods, such as the finite element (FE) method, discontinuous Galerkin (DG) method [10][11][12], and hybridizable discontinuous Galerkin (HDG) method, have been developed to solve Maxwell's equations. Discretization by either an FE method [13] or an HDG method [14][15][16] can yield a large sparse discretized linear system. It is still difficult to solve the resulting system of linear algebraic equations by either a direct solver or a standard preconditioned iterative method. On the other hand, the domain decomposition method (DDM) of Schwarz-type is considered to be one of the most efficient solving strategies for Helmholtz-type problems and then has been extended for the time-harmonic Maxwell's equations in [17][18][19]. Moreover, DDMs should be very suitable for implementing the high-performance parallel computations, because they can decompose the large-scale and complex boundary value problems (BVPs) into a series of small-scale and simple BVPs that can be solved separately. In short, DDMs are usually employed to deal with such large-scale problems [20]. In [17], where the permittivity is a real number, an optimized Schwarz method combined with an HDG method discretization was used for EM problems, and the coupling between the Schwarz method and the HDG method was shown to be natural.
In this paper, the permittivity is a complex number in the Drude model, which adds the complexity to the optimized Schwarz method. We derive an optimized transmission condition and a formulation of the spectral radius of the iteration matrix of the Schwarz method. Furthermore, the parameters for an optimized transmission condition are discussed and tested. The subdomain problems are discretized by an HDG method [16] and the resulting linear systems can be solved by a sparse direct solver. The popular Krylov subspace method, namely GMRES [21], is considered to accelerate the solution of the interface linear system.
There are five sections in the rest of this paper. First of all, the Drude model and some notations are briefly introduced in Section 2. The discretization of Drude model by an HDG formulation are described in Section 3. In Section 4, we present the formulations of a Schwarz algorithm and study the parameters of an optimized transmission condition in a two subdomains setting. Numerical tests are presented in Section 5 to show the effectiveness of the proposed method. Finally, we draw some concluding remarks in Section 6.

Problem and Notations
In this section, we will introduce several concepts and notations which are essential for our present study.

Maxwell's Equations with Drude Model
We consider the 2D Maxwell's equations in the frequency domain with a first-order Silver-Müller absorbing boundary condition (i.e., an artificial absorbing boundary condition) [22]      where i = √ −1 stands for the imaginary unit, ω refers to the angular frequency of the light wave, ε, ε 0 , and µ 0 represent the relative permittivity, the permittivity of free space, and the permeability of free space, respectively, E = E x , E y and H = H z denote the electric and magnetic fields, the superscript "· inc " means the incident field, and n is the outward unit normal vector. The differential operators in this 2D setting are CurlH = (∂ y H, −∂ x H) and curlE = ∂ x E y − ∂ y E x . The computational domain is denoted by Ω, and the artificial absorbing boundary is denoted by Γ a [22].
Note that ε = 1 in the free space. According to the Drude model, ε accounts for the interactions between the time-varying electric field and the electron gas [2]. It varies with the angular frequency of the incoming light, i.e., where ω p denotes the bulk plasma frequency of the material and γ is a damping constant.
In [23], the authors consider the time-domain Maxwell-Drude model with a DG timedomain method.

Notations
We write here a triangulation T h of Ω with K denoting an element of discrete mesh. F h , F I h , and F B h represent the set of all edges of T h , the set of all the edges of T h associated with the nanostructure, and the union of all the boundary edges of T h , respectively. For an element K 1 ∈ T h and its adjacent element K 2 , F = K 1 K 2 is the common edge of K 1 and K 2 . Let (v 1 , v 1 ) be the traces of (v, v) on F from the interior of K 1 and (v 2 , v 2 ) be the traces of (v, v) on F from the interior of K 2 . n 1,2 and t 1,2 stand for the outward unit normal vectors to K 1,2 and the unit tangent vectors to the boundaries ∂K 1,2 , respectively, so we have t 1 × n 1 = 1 and t 2 × n 2 = 1. On the face F, {·} and · can be defined as For each K ∈ T h (F ∈ F h ), p ≥ 0 refers to the local interpolation order, and P p (K) (P p (F)) refers to the space of polynomial functions of degree at most p. We define the discontinuous FE spaces V Note that L 2 (Ω) represents the space of a squared integrable functions over Ω, where Γ m satisfies Γ m ∪ Γ a = ∂Ω, Γ m ∩ Γ a = ∅. Note that ∂Ω denotes the boundary. M p h consists of the functions which are not continuous at its ends, but continuous on an edge. For a domain D in R 2 , u, v in (L 2 (D)) 2 and u, v in L 2 (D), (u, v) D refers to D u ·vdx where · denotes the complex conjugation and (u, v) D stands for the inner product D uvdx. On an interface F, u, v F stands for the inner product F uv ds. So on the whole domain Ω, we have

HDG Formulations
We consider an approximate solution Use the Green's formula for the above equations and replace the boundary terms with the numerical traces E h , H h . One can have A proper choice of numerical trace E h , H h affects the correctness and the convergence of the discrete problem (3). According to the ideas in [16], we choose a hybrid variable λ h ∈ M p h , and set E h and H h as follows where τ > 0 is the local stabilization parameter. Considering the contributions of Equation (3) over all elements, the artificial absorbing boundary condition in the formulation of this conservativity condition and enforcing the continuity of the tangential component of E h , we have Using the variable λ h to express E h and H h , then we can obtain a global problem with only since all the interior faces satisfy the conservativity condition, we have and inserting (7), we can obtain where the superscript and the subscript 1 and 2 denote the values from the two elements coupled by edge. Therefore, the numerical traces can be expressed as and using a similar method we also can obtain

An Optimized Schwarz Method
To introduce the Schwarz method, we divide the domain Ω into Ω 1 and Ω 2 , and note that Ω 1 and Ω 2 are two non-overlapping subdomains. One can easily obtain the case that the domain Ω is divided into many subdomains, because the transmission condition only involves the adjacent subdomains. For the given initial guesses (E l,0 , H l,0 ), l = 1, 2, on the interface between the subdomains, we can compute (E l,n+1 , H l,n+1 ) from (E l,n , H l,n ) with the following Schwarz method [17].
in Ω 1 , in Ω 2 , The computational domain Ω is displayed in Figure 1, where Γ 1,2 denotes the interface between the two adjacent subdomains. Γ 1a and Γ 2a denote the artificial absorbing boundary in each subdomain. The transmission condition is defined as B n l (E, H) = S l H + n × E, In the following, we will show that the coupling between the Schwarz method and the HDG method is natural. We set K 1 ∈ Ω 1 , K 2 ∈ Ω 2 to be two elements sharing a common face F between two adjacent subdomains. We denote n 1,2 as the outward unit normal vectors to K 1,2 and impose Then one can set S 1 = τ 2 , S 2 = τ 1 , and we have the transmission condition B n (E, H) = S l H + n × E. , l = 1, 2, the transmission condition will be the Silver-Müller condition, where Re(·) takes the real part of a complex number. We call it the classical transmission condition [17].

Optimized Parameters for Optimized Schwarz Method
In the following, we try to give an analysis of the theoretical spectral radius of the iteration matrix of the Schwarz iteration, which is similar to that in [24]. Suppose that Ω 1 is the left half plane (−∞, 0] × (−∞, +∞) and Ω 2 is the right half plane (0, +∞) × (−∞, +∞). Taking a Fourier transform in the y direction with Equation (1), we obtain where k is the Fourier coefficient. According to the first equation of Equation (13), replacing E x with H, we obtain Because of the radiation condition, the solution of Equation (14) in Ω l (l = 1, 2) is given by where (14), v 1 , v 2 are their corresponding eigenvectors, and the coefficients a 1 , a 2 are uniquely determined by the transmission conditions. Set (15) and (16) into the last equation of Equation (11), we have S 1 a n+1 1 iωε 0 ε λ e λx + a n+1 1 e λx = −S 1 a n 2 iωε 0 ε λ e −λx + a n 2 e −λx , −S 2 a n+1 2 iωε 0 ε λ e −λx − a n+1 2 e −λx = S 2 a n 1 iωε 0 ε λ e λx − a n 1 e λx . At the n-th step of the Schwarz algorithm with x = 0, the coefficients a 1 , a 2 satisfy the system a n+1 then we have the spectral radius ρ in the form In order to derive an optimized transmission condition, one can set the second-order approximation of the operator S l = α l + β l ∂ 2 τ , l = 1, 2, see Reference [25], where ∂ 2 τ denotes the second-order derivative along with the interface and α l , β l are the parameters to be determined [17]. With the zeroth order approximation of the operator S l , i.e., S 1 = α 1 , S 2 = α 2 , where α l (l = 1, 2) are two complex numbers, then Therefore one can consider the optimization problem [24,25] as follows to determine the optimized parameters α l (l = 1, 2) α * l = arg min α 1 ,α 2 (max(|ρ|)), l = 1, 2.
Unfortunately, it is difficult to solve this optimization problem explicitly because this problem is an open problem. For classic Maxwell's equations with real permittivity, α l = (iω) (−1) (p l + ip l ) is often used [17,19,24], where where k 1 and k 2 are the highest and lowest possible frequency allowed [24]. However, this choice does not work well for the Drude model. We present four possible guesses of α l in Table 1 which lead to different optimized transmission conditions. The term "classical" in Table 1 represents classical transmission condition mentioned in Remark 1. Case3 is the above common choice. Table 1. Spectral radius ρ and parameter α.
As seen from Table 1, we find that the Schwarz methods with the last two cases do not converge at all. Case1 is theoretically the best choice which can be seen from the spectral radius ρ. In the next sections, numerical results will confirm this theoretical observation.

HDG Discretization
Let Γ i a = Γ a ∂Ω i and Γ i,j = ∂Ω i Ω j , i, j = 1, 2. Using the optimized transmission conditions on the interface Γ i,j , one can make discrete Equation (11) where the quantity 1 − ω 2 p ω(ω+iγ) is just equal to ε(ω) mentioned in Equation (2). For an element K e , we rewrite the local solution (E h , H h ) and hybrid variable λ h like the form in Equation (17) of [17], i.e., at this moment, the discretized system is transformed to solve the following problem where g represents the according degrees of freedom (DOFs) on Γ 1,2 and i indicates the according DOFs in Ω 1 or Ω 2 . Moreover, this resulting linear system (22) is large, sparse but complex non-Hermitian, so the sparse direct solvers are always very expensive and prohibitive [26,27]. Thus, in the next section we use the Krylov subspace methods [21], which only depend on the information of the coefficient matrix-vector products. Now we summarize the main steps of solving the Drude model with an optimized Schwarz method discretized by HDG method as follows: • The model problem is split into some sub-problems with the corresponding subdomains which are discretized using an HDG method; • Then we solve the resulting system of linear algebraic Equations (6) in each subdomain by a sparse direct solver; • Finally, for the interface system between the two subdomains, solving the resulting linear systems (22) in the domain is accelerated using a Krylov subspace method.

Numerical Tests
In this section, we present two numerical results to show that the optimized Schwarz method is effective. All the numerical simulations are implemented in MATLAB R2012a and performed on a desktop with an AMD A6-6310 APU with AMD Radeon R4 Graphics CPU of 1.80 GHZ and 4.0 GB memory. We only employ the zeroth order approximation of S l in our tests. We use Gmsh (see https://gmsh.info/ (accessed on 15 September 2022)) to decompose the domains. In the following, "HDG-P 1 " denotes the HDG discretization method relying on a nodal Lagrange basis interpolation of order p = 1. For the numerical solution, we set the stopping criterion of the iteration process as 10 −6 ; that is, when the relative residual r k 2 r 0 2 < 10 −6 , then we stop the iteration of the Krylov subspace method, namely GMRES (DD-Gmres).

Cylindrical Nanowire Problem
In this test, we set the radius of the cylinder to 20 nm. The interband transitions are ignored here. The computational domain Ω = [−L, L] 2 is a square with L = 200 nm. We impose the artificial absorbing boundary condition on the boundary of Ω. In our test, we set ω p = 8.65 × 10 15 , γ = 0.01ω p [28]. A typical subdomain decomposition is shown in Figure 2a. Meshes for the cylindrical nanowire problem are shown in Figure 2b. In Table 2, we divide the domain into 3288 elements with 1645 nodes. In Table 3, we divide the domain into 13,024 elements with 6513 nodes. We chose ω = 1.4 ω p and ω = 0.4 ω p since they are closest to the resonance frequency of the material. In fact, the Maxwell's equations are often considered to be more difficult to solve around the reso-nance frequency than other frequencies [5]. The results of different optimized parameters with different subdomains at ω = 1.4 ω p and ω = 0.4 ω p are presented in Tables 2 and 3. The results show that the parameters with case3 do not work well for the Drude model problem. With the same number of subdomains, case1 and case2 converge much faster than the classical Schwarz method, and case3 which is consistent with Tables 2 and 3. Furthermore, case1 outperforms case2 under the same number of subdomains. For each case, the number of iterations increases with the increasing of number of subdomains.  The domain in the following tests contains 13,024 elements with 6513 nodes. We show how the number of iterations required by the GMRES method varies with the number of subdomains at ω = 1.4ω p and ω = 0.4ω p in Figure 3. As seen from Figure 3, we can find that the number of iterations for GMRES increases with the number of subdomains. Additionally, case1 increases more slowly than that of the classical Schwarz method. For two subdomains, how the interpolation order p in the HDG formulation affects the convergence at the frequency ω = 1.4ω p is shown in Table 4. We observe that the number of iterations increases with the increasing of p for each case. Furthermore, case1 outperforms the classical Schwarz method under the same p. Field distributions at the frequency ω = 0.4ω p are presented in Figure 4. Note that in the above tests we set the value of the local stabilization parameter τ = 1. The local stabilization parameter τ varies how it affects the number of iterations at ω = 1.4ω p ; this is shown in Table 5.   Table 5. The influence of the local stabilization parameter τ in the HDG method on the number of DD-gmres iterations.

Dimer of Cylindrical Nanowires
The computational domain is a rectangle with length and width of 300 nm and 200 nm, respectively. We divide the computational domain into 2226 elements with 1114 nodes. Then we consider the plasmonic dimer structures with small gaps [3,29]. We use the parameters in [3]: ω p = 1.34 × 10 16 , γ = 1.14 × 10 14 . The radius of the cylinder is 30 nm. We present a typical subdomain decomposition in Figure 5.  Table 6 shows how the number of iterations varies with the number of subdomains at ω = 0.4 ω p . As we can see from Table 6, case1 converges much faster than the classical Schwarz method with the same number of subdomains. For each case, the number of iterations increases as we increase the number of subdomains; how the number of iterations varies with the number of subdomains is shown in Figure 6. For two subdomains, how the interpolation order p in the HDG formulation affects the convergence is shown in Table 7. As seen from Table 7, we can find that the number of iterations increases as we increase the interpolation order. The field distributions at this particular frequency are displayed in Figure 7.

Conclusions
In the previous study [17], where the permittivity is a real number, the authors have solved the Maxwell's equations with an optimized Schwarz method discretized by an HDG method, which performs well. In the current paper, the permittivity is a complex number in the Drude model, which adds the complexity to the optimized Schwarz method. We employ an optimized Schwarz method combined with an HDG discretization to solve the local optical response model. The domain is arbitrarily divided into several subdomains. New transmission parameters are proposed and tested. Numerical tests show that the optimized Schwarz method with a proposed parameter works quite well for Drude model problems.
For future work, it is noted that the coefficient matrix of the resulting linear system (22) is complex symmetric, which means that some particular Krylov subspace solvers [30,31] with suitable preconditioners for such linear systems can be employed to reduce the computational cost. In addition, it will be meaningful to extend the proposed method for solving three-dimensional model problems.  Data Availability Statement: Data will be made available on request.

Acknowledgments:
The authors would like to thank the referees and academic editor for their insightful suggestions and comments that lead to a significant improvement of this article.

Conflicts of Interest:
The authors declare no conflict of interest.