A Nonlocal Operator Method for Partial Differential Equations with Application to Electromagnetic Waveguide Problem

A novel nonlocal operator theory based on the variational principle is proposed for the solution of partial differential equations. Common differential operators as well as the variational forms are defined within the context of nonlocal operators. The present nonlocal formulation allows the assembling of the tangent stiffness matrix with ease and simplicity, which is necessary for the eigenvalue analysis such as the waveguide problem. The present formulation is applied to solve the differential electromagnetic vector wave equations based on electric fields. The governing equations are converted into nonlocal integral form. An hourglass energy functional is introduced for the elimination of zeroenergy modes. Finally, the proposed method is validated by testing three classical benchmark problems.


Introduction
For the analysis of complex structures in engineering, mesh generation remains a laborious and time-consuming process, which often requires human interaction with meshing program and corrections for local mesh. To alleviate the mesh entanglement, socalled meshless or meshfree methods have been proposed during the 1990s [Viana and Mesquita (1999) ;Xuan, Zeng, Shanker et al. (2004); Razmjoo, Movahhedi and Hakimi (2011);Nicomedes, Bathe, Moreira et al. (2017)]. Meshless methods have been afterwards developed, enriched and applied to analyze a wide variety of problems in engineering including electromagnetic problems. For example, Viana and Mesquita applied the meshless Moving Least Square Reproducing kernel Method to solve twodimensional static electromagnetic problems. Liu et al. [Liu, Yang, Chen et al. (2004)] proposed an improved formulation of the element-free Galerkin method for electromagnetic field computations and studied the selection of the weight function, the treatment of imposing boundary conditions and interface conditions. The present original nonlocal operator method (NOM) proposed by the authors can be used to solve the partial differential equations by replacing the local operators in PDEs with newly defined nonlocal operators. The nonlocal operator is defined in the integral form based on the nonlocal interaction but converges to the local operator in the continuous limit. The nonlocal operator method is consistent with the variational principle and the weighted residual method, based on which the tangent stiffness matrix can be obtained with ease. The nonlocal operator method can be viewed as a generalization of dual-horizon peridynamics [Ren, Zhuang, Cai et al. (2016); Ren, Zhuang and Rabczuk (2017)] or peridynamics [Silling (2000); Silling, Epton, Weckner et al. (2007)]. In this paper, we develop and apply the nonlocal operator method to electromagnetic problem. Electromagnetic analysis has been an indispensable part of many engineering and scientific study since Maxwell established a unified electromagnetic field theory-the Maxwell equations-in the 19th century. The Maxwell equations describing electromagnetic waves have numerous applications including radar, remote sensing, bioelectromagnetics, wireless communication and optics, just to name a few. Several computational methods have been developed for the solution of the Maxwell equations including the method of moments [Gibson (2007)], finite element method [Jin (2015)], time domain finite difference method [Taflove and Hagness (2005)], ray theory [Deschamps (1972)], meshless/meshfree methods [Ho, Yang, Machado et al. (2001)], asymptotic-expansion methods [Bouche, Molinet and Mittra (2012)] and eigen expansion method [Chew, Jin, Lu et al. (1997)]. The finite element method and time domain finite difference method can capture complex shapes and inhomogeneous materials while the moment method is well known for its high precision and efficiency for mainly linear problems and simple geometries [Jin (2015)]. The Finite-Difference Time-Domain (FDTD) method [Yee (1966)] is a method for directly solving the Maxwell equation in the time domain. The FDTD method has merits such as the explicit time integration, high efficiency in storage and natural parallelization. The method of moments (MoM) or boundary element method (BEM) is a numerical computational method of solving linear partial differential equations which have been formulated as integral equations (i.e., in boundary integral form) [Harrington (1993)]. This method possesses high accuracy but requires artificial intervention to handle the integral equations. In addition, this method is only applicable for regular shapes and homogenous materials. The multi-layer fast multipole technology based on the moment method is often employed for the purpose of computational efficiency. The accuracy of the finite element and FDTD is lower than the moment method since both finite element and time domain finite difference have numerical dispersion errors, which does not occur for the moment method. The purpose of this paper is to develop a framework of nonlocal operator method exploiting variational principles and to reformulate the electromagnetic governing equations. Therefore, the local differential equation is converted into nonlocal integral form. The content of the paper is outlined as follows: The governing equations for electromagnetic fields in the time domain and frequency domain are described in Section 2. The concept of nonlocal operator method including definitions of the nonlocal curl and gradient operators are introduced in Section 3. Furthermore, the nonlocal formulation based on the variation of nonlocal operators in discrete form are presented in details to finally obtain the consistent tangent stiffness. In Section 4, we convert the differential equation into the nonlocal form. In Section 5, the hourglass energy functional to remove zero-energy modes is proposed, which are inherited from particle-based formulations based on nodal integration. The residual and tangent stiffness matrix of the hourglass functional is also derived. The nonlocal integral form of the electromagnetic equations in the time-domain is derived in Section 6. Three benchmark problems are solved in Section 7 to verify the method. Finally, we conclude our manuscript in Section 8.
The divergence-free requirement in Eq. (1d) can be imposed for example with the penalty method [Rahman and Davies (1984)], vector finite elements [Whitney (2012);Nédélec (1980); Hano (1984)] or specially designed shape functions as presented in [Evans and Hughes (2013)]. The constitutive relations can be written as , where the constitutive parameters  , µ and σ denote, respectively, the permittivity (farady/meter), permeability (henrys/meter), and conductivity (siemens/meter) of the medium. For the time-harmonic fields with a single frequency, the time dependent parts of Maxwell's equations can be written in simplified form as , j j j ω ω ωρ where j is the imaginary unit and ω is the angular frequency. The vector wave equations for E can be obtained by eliminating H from Eq. (3b) and considering the constitutive relations Eq.
The boundary conditions for equations based on E are Consider a domain as shown in Fig. 1(a), Let x be the spatial coordinates in the domain Ω ; : ′ = − r x x is the spatial vector, the relative distance vector between x and ′ x ; : ( , ) t = F F x and : ( , ) t ′ ′ = F F x are the electric field vectors for x and ′ x , respectively; : ′ = − r F F F is the relative electric vector for vector r .
The vector wave equations can also be formulated by using only H . In this paper, we will employ the vector wave equations based on electric fields.

Support x
 is the domain where the nonlocal operator is defined, and any spatial point ′ x in support forms spatial vector r . Support x  is usually presented by a spherical domain with a radius of h x . A point interacts with other points which fall inside the support of that point through nonlocal interactions. In order to define the nonlocal operator, the shape tensor is defined as which is symmetric. The prerequisites of shape tensor are that it shall be invertible, which can be satisfied usually when enough particles fall inside the support. Numerical example shows that the minimal number of neighbors in support is 2 and 3 for two-dimensional and three-dimensional problems, respectively.
Dual-support is defined as the union of points whose supports include x , denoted by On the other hand, ′ r is the spatial vector formed in ′ x  . One example to illustrate the support and dual-support is shown in Fig. 1(b).

Nonlocal operators and definitions based on the support
In nonlocal operator method, key operators include the nonlocal operators for divergence, curl and gradient since they can be used to replace the local operators in the partial differential equations. We use ∇  to denote the nonlocal operator, while the local operator is ∇ . The nonlocal gradient of field F for point x in support x  is defined as The nonlocal curl of field F for point x is defined as The nonlocal divergence of field F for point x is defined as The field value near a point ′ x can be approximated by Taylor series expansion by neglecting higher order terms as , or .
Inserting Eq. (11) into the RHS of Eqs. (8), (9), (10) and integrating in support x  , it can be shown that the nonlocal operators are identical to the local operators. For example, When ′ x is close enough to x or when support x is small enough, the nonlocal operator can be considered as the linearization of the nonlinear field. The nonlocal operator converges to the local operator in the continuous limit. On the other hand, the nonlocal operator defined by integral form, still holds in the case where strong discontinuity exists and the local operator cannot be defined. The local operator can be viewed as a special case of the nonlocal operator.

Variation of nonlocal operators
The nonlocal operators defined above are in vector or tensor form. The variation of the nonlocal operators leads to a higher-order tensor form, which is not convenient for implementation. We need to express the high-order tensor into to vector or matrix form. Before we derive the variation of nonlocal operator, some notations to denote the variation and how the variations are related to the first-and second-order derivatives is to be discussed. Assuming a functional ( , ) x are unknown functions in unknown vector [ , ] u v , the first and second variation can be expressed as : It can be seen that the second variation 2 ( , ) u v δ  is the double inner product of the Hessian matrix and the tensor formed by the variation of the unknowns, while the first variation ( , ) u v δ  is inner product of the gradient vector and the variation of the unknowns. The gradient vector and the Hessian matrix represent the residual vector and tangent stiffness matrix of the functional, respectively, with unknown functions , u v being the independent variables, [ , ] 2 [ , ] ( , ) [ , ] ( , ) .
The inner product or double inner product indicates that the location of an element in the residual or the tangent stiffness matrix corresponds to the location of the variation of the unknowns.
In this paper, we use a special variation δ The special first-order and second-order variation of a functional lead to the residual and tangent stiffness matrix directly. The traditional variation can be recovered by the inner product of the special variation and the variation of the unknown vector.
The variation of ∇ ⋅ x F  is given by The number of dimensions of δ ∇ ⋅ x F  is infinite, and discretization is required.
After discretization of the domain by particles, the whole domain is represented by on nodal level is sometimes called the nodal assembly. In the following section, we mainly discuss the special variation of the nonlocal operator and functional, while the actual variation can be recovered with ease.

The variation of
where k is the index of particle k i j N ∈ . The minus sign denotes the reaction from the dual-support, which guarantees the regularity of the stiffness matrix in the absence of external constraints. The nodal assembly for the variation of the vector cross product can be finally obtained by The indices of R correspond to the locations in × F .

Similarly, the variation of
where k is the index of particle k i j N ∈ . The sub-index of δ ∇ x F  can be obtained by the way similar to Eq. (27).

Waveguide
A waveguide is a structure that guides waves, such as electromagnetic waves or sound, with minimal loss of energy by restricting expansions to one or two dimensions. The study of waveguide requires the eigenmode of the wave propagation which in turn requires the tangent stiffness matrix of the field. In this section, we derive the tangent stiffness matrix in the framework of variational principle using nonlocal operator method. The partial differential equation with boundary conditions for the waveguide problem can be expressed as The sum of Eq. (31) and Eq. (32) is Applying second vector Green's theorem in Eq. (34) The divergence-free condition is enforced by the penalty method, so the functional becomes where p is the penalty parameter which is set to 1 in our examples as suggested in Rahman et al. [Rahman and Davies (1984)]. Finally, the eigenvalue problem of the waveguide problem reads 1 2 ( ) 0, 0 on , 0 on .
Eq. (38b) is the electric wall and Eq. (38c) is the magnetic wall, which is enforced for the sake of better accuracy and the elimination of some spurious solutions. For rectangular waveguide, the normal direction is parallel to a certain axis, for example (1, 0, 0), can be applied the same as Dirichlet boundary conditions.
and its first variation is written as Consider the first variation of all particles, and let The relation ( ) ( ) ( ) ⋅ × = ⋅ × = ⋅ × a b c c a b b c a is used in the third step. In the third and fourth steps of above derivation, the dual-support has been employed, i.e., in the third step, the term δ ′ x F is the vector from x 's support, but is added to particle ′ x . In the fourth step, all the terms with δ x F are collected from other particles whose supports contain x and therefore form the dual-support of x . For any δ x F , the first order variation ( ) 0 δ = E  leads to the nonlocal form of the governing equation of the waveguide problem: When the particle's volume 0 V ′ ∆ → x , the continuous form is Eq. (43) is the nonlocal governing equation of the waveguide on the electric field. For the eigenvalue problem, the stiffness matrix is required.
The special second variation of ( ) x E  leads to the tangent stiffness matrix, leading to the generalized eigenvalue problem 2 0 Note that the nodal integration of the above integrals results in hourglass modes which can be removed by introducing so-called hourglass energy, which will be addressed in the next section.

Hourglass energy functional
In order to remove the hourglass or zero-energy modes, a penalty term is added to achieve the linear completeness of the electric field, in which the penalty is proportional to the difference between the current value of a point and the value predicted by the field gradient. The electric field in the neighborhood of a particle is required to be linear. Therefore, it has to be exactly described by the gradient of the electric field, and the hourglass modes are identified as that part of the electric field, which is not described by the electric gradient. The difference between the current vector r E and the predicted vector by the field gradient ( : The above definition of the hourglass energy is similar to the variance in probability theory and statistics. In the above derivation, we have used the relations, : : ( ), : , : tr( ) where the capital letter denotes the matrix and small letter the column vector. The purpose of m K is to make the energy functional independent of the support since the shape tensor K is involved in : It should be noted that the zero-energy functional is valid in any dimensions and there is no limitation on the shape of the support. Consider the variation of the zero-energy functional The residual of the hourglass mode is required in the explicit time integration method. In this paper, we only discuss the implicit analysis.
The electric flux of the hourglass mode for one vector r is given by ( Eq. (50) is an explicit formula for the hourglass flux. The term on δ E is the hourglass term from its support, while the terms on δ ′ E are the hourglass terms for the dual The second variation of the zero-energy functional is its stiffness matrix on one point. The global tangent stiffness matrix for hourglass energy functional can be assembled by The above equations indicate, when the electric field is consistent with the field gradient, then the hourglass energy residual is zero.
Once the hourglass mode is eliminated, the residual of the hourglass mode is zero. The stiffness matrix of the hourglass mode overcomes the rank deficiency of the matrix system when nodal integration is used. The generalized eigenvalue problem becomes 2 0 ( ) 0, where E is the eigenvector for all unknowns.

Nonlocal operator method for electromagnetic in the time domain
Consider a volume Ω bounded by the surface S . The electromagnetic field generated by an electric current density x J satisfies the Maxwell equations. Eliminating the magnetic field with the aid of the constitutive relations, the curl-curl equation for the electric field E is obtained by Jin [Jin (2015)]: .. . .
σ µ In order to obtain the equivalent nonlocal form of Eq. (54), let us consider Eq. (40) from the previous section. From the variational derivation of the waveguide problem, Based on Eq. (43), one can get the nonlocal form of The Dirichlet boundary conditions are 1 , , where x P is the specified electric wall on point x .
The Neumann boundary condition on 2 S can be written as Finally, the central difference scheme can be used for the time integration yielding ) The tangent stiffness matrix is obtained as 2 [ 10,10] ( ) is the nonlocal operator in 1D. The hourglass energy functional in 1D can be obtained with the procedure similar to that in Section 5.
We calculate the lowest eigenvalue and compare the numerical result with 0 0.5 λ = . The convergence plot of the error with different weight function and discretizations is shown in Fig. 2. It can be seen that the convergence rate is around 2. The weight function and inhomogeneous discretization have limited effect on the convergence. Good agreement is observed between the numerical result and the exact solution. weight function 1/r 3 uses an inhomogenous discretization in Fig. 3; the particle spacing in dual-form is selected as the minimal particle spacing in the discretization The discretization of the dual-form is given in Fig. 3. The first three wave functions are given in Fig. 4.
where φ denotes potential, ρ is the charge density of the domain, n is outward normal direction of the boundary, Ω is the solution domain, and its boundary q φ ∂Ω = Γ ∪ Γ , φ is the specified potential value at the boundary φ Γ , q is the potential derivative value given on the boundary q Γ .
In the simulation, the boundary conditions are applied with the penalty method. The equivalent energy functional is where 6 1 10 α = × is the penalty coefficient. The stiffness matrix can be obtained the similar way in Section 5. In order to validate the accuracy of nonlocal operator formulation on the electrostatic field problem, we calculate the electrostatic field of a rectangular plate, which is a benchmark problem with the analytical solution given in Eq. (67). In this example, the potential on the upper and lower sides is 0, and the right side is The support radius is selected as The hourglass penalty of 1 hg p = is used in all simulations. The plate is discretized with different mesh densities to test the convergence. The contour plot of the electric potential from the numerical solution and analytical solution are shown in Fig. 6 and Fig. 7, respectively. A satisfactory agreement is observed between Fig. 6 and Fig. 7. The L2 norm of the electric potential decreases with the refinement of the mesh with the convergence rate of 0.8709 r = , as shown in Fig. 8.

Rectangular Waveguide problem
A hollow waveguide is a transmission line that looks like an empty metallic pipe. It supports the propagation of transverse electric (TE) and transverse magnetic (TM) modes, but not transverse electromagnetic (TEM) modes. There is an infinite number of modes that can propagate as long as the operating frequency is above the cutoff frequency of the mode. The notation TEmn and TMmn are commonly used to denote the type of the wave and its mode, where m and n are the mode number in the horizontal and vertical directions, respectively. The mode with the lowest cutoff frequency is called the fundamental mode or dominant mode. For a hollow rectangular waveguide, the dominant mode is TE10. The analytical solution for the E-field in the TE mode is expressed as The electromagnetic analysis of a rectangular waveguide is well known [Pozar (2009)]. Let us focus on the results used to verify our formulation, i.e., 10.16 b = mm and 40 l = mm; the boundaries in blue illustrate the electric walls and the red boundary is the magnetic wall, see Fig. 9. The domain of waveguide is discretized with two different particle spacings, as shown in Fig.  10. The support is selected as 2.2 h x = ∆ and weight function 2 ( ) 1/ w r r = . Figure 9: Section of a rectangular waveguide, where a=22.86 mm, b=10.16 mm and l=40 mm. Blue boundary denotes electric wall and red boundary is magnetic wall The calculated frequencies are given in Tab. 1. The error in the frequency for Case 2 is less than 4%. Good agreements are obtained between the numerical results and theoretical results. The modes of the E-Field for two cases are shown in Fig. 11 and Fig. 12. (c) TE01 Figure 12: The TE modes for case 2 8 Conclusion In this paper, we proposed a nonlocal operator formulation for electromagnetic problems employing variational principles. The formulation is implicit and provides the tangent stiffness matrix, which is needed for the solution of the eigenvalue problem. We presented a scheme for assembling the global stiffness matrix based on nonlocal operators. The nonlocal form of the electromagnetic vector wave equations based on the electric field is obtained by means of the variational principles. Three numerical examples, including the Schrödinger equation in 1D, electrostatic field problem in 2D and waveguide in 3D are tested and show good agreement to available analytical solutions. In the future, we intend to solve also the transient problem and study problems involving strong discontinuities which are one of the key strength of nonlocal operator method.