This article studies numerically the electromagnetic scattering properties of three dimensional (3D), arbitrary shaped dielectric resonator antennas which are composed of single and multi - layered (composite) dielectric materials . Using the equivalence principle and the integral equation techniques, we first derive a surface integral equation (SIE) formulation which produces well - conditioned matrix equation. We then develop an algorithm to speed up the matrix -vector multiplication s by employing the well -known method of moments (MoM) and the multi level fast multipole method (MLFMM) on personal computer (PC) clusters. To solve the obtained integral equations, we apply a Galerkin scheme and choose the basis and testing functions as Rao -Wilton - Glisson (RWG) defined on planar patches. Finally, we present s ome 3D numerical examples to demonstrate the validity and accuracy of the proposed approach. INTRODUCTION The dielectric resonator antennas (DRAs) have received significant attention d uring the last two decades because of their advantages over conventional metallic antennas which suffer from problems like power loss, low bandwidth, low radiation efficiency and fabrication difficulties [1]. Figure 1 shows the configuration of a multi -layer dielectric antenna which comprises several la yers of varying dielectric constant. Figure 1: multi -layer dielectric resonator antenna There is a clear need for computational and numerical analysis and design tools for calculating the performance and optimizing the parameters of dielectric antennas prior to costly prototype development. These computational methods are the building blocks of the commercial simulation software for modeling of three -dimensional me tallic and dielectric antennas. One of the most general tactics to the analysis of dielectric antennas is the finite -difference time -domain (FDTD) method. The FDTD has been used by a number of authors to analyze dielectric antennas [2-5]. In [6] it is fo und that the FDTD is a suitable choice for axially symmetric antennas. However, in practice, the practical structures are often asymmetric such as slotted and aperture dielectric antennas. Although the FDTD is capable of modeling and analyzing an extensive range of arbitra ry-shaped dielectric antennas, the complexity cost is very environment dependent and suffers from stair -casting problem [7]. The s urface integral equation method is one of the most popular numerical methods in electromagnetic analysis of homogenous dielectric and composite metallic and dielectric objects. In [8] SIE based on the method of moments (MoM) has been used to simulate the interaction of light and plasmonic nanostructures. SIEs are also popular in microwave regime, and have been extensively used for the analysis of periodic lossy [9] or metallic [10] structures. In this paper, we first derive a surface integral equation formulation to analyze electromagnetic scattering and radiation problems of all -dielectric resonator antennas. By discretizing the obtained integral equations, we acquire dense matrix equations. We then apply the well -known method of moments (MoM) along with the multilevel fast multipole method (MLFMM) to speed up vector - matrix multiplications. Proposed SIE can be extremely useful by enabling de tailed analysis of novel antenna designs before their actual realizations and preventing waste of sources and time required for building and testing prototypes. INTEGRAL EQUATION FORMULATION Consider the case in figure 2 involving the problem of electromagnetic scattering by an arbitrary shaped, homogenous dielectric object B characterized by electrical properties ( , )ii which is surrounded by an infinite medium having ( , )oo . Figure 2: A homogenous region and boundary surfaces Using the equivalence principle, the inner and outer normal (N) and tangential (T) electric and magnetic fields can be considered as: N-EFIE(l)  ˆ ˆ ˆ J M Einc l l l l l l ln n n −  =−  (1) T-EFIE(l)  tanˆ ˆ ˆ ˆ ˆ ˆ M J ( E )inc l l l l l l l l l ln n n n n n    −   =   (2) N-MFIE(l)  1ˆ ˆ ˆ J M Hinc l l l l l l ln n n − −  =−  (3) T-MFIE(l)  1 tanˆ ˆ ˆ ˆ ˆ ˆ J M ( H )inc l l l l l l l l l ln n n n n n −−   −   =   (4) In (1) -(4), ˆlndenotes the inner unit normal for l=i and outer unit normal for l=o, l l l  = is the impedance of the host medium with   ( ), ( ) l i inner o outer= , (E , H )inc inc ll denote the incident electric and magnetic fields respectively, and Jl and Ml are equivalent electric and magnetic currents on the surface of the object which are relat ed to the total electric ( E ) and magnetic ( H ) fields as ˆ J ( ) H( )llr n r= (5) ˆ M ( ) E( )llr r n= (6) Figure 3: Boundary conditions for two dielectrics For figure the boundary conditions can be written as 1 1 1 2 2 2ˆˆ J ( ) H ( ) J ( ) H ( )r n r r n r   =  =− =−  , (7) 1 1 2 2 2ˆˆ M ( ) E( ) M ( ) E ( )r r n r n r   =  =− =  (8) The operators and  operators a re defined as 2X [ X( ) .X( ) ]G ( , )l l lssjk r ds k ds r r r−       = +   (9) X X( ) G ( , )lsr r r ds     =  (10) Where l l lk = is the wavenumber and we define G(r,r ) ( )4ljk ReR r rR− = = − (11) Which is the homogenous Green’s function . Among the various possibilities of combination of the normal and tangential electric and magnetic fields, the following mixed formulation is called the electric and magnetic current combined field integral equation (JMCFIE ) [11, 12], i.e.      22 12 11 22 2 11T- EFIE(l) ( 1) N- MFIE(l) T- MFIE(l) ( 1) N- EFIE(l)l l ll l l ll − == ==+−  −−  (12) Using boundary conditions, it can be written as: 1 11 12 tan 21 22 tan[ ] [ ] Jˆ (E ) H [ ] [ ] Mˆ (H ) Einc inc N N N N ll inc inc N N N N llZZ n ZZ n −   −  =  +    (13) For a composite dielectric object involving B -1 regions B 1, B2,…, B n with constitutive parameters ( , )nn located in a host medium 00( , ) , the surface integral equation can be derived as 1 1 1 1 11, 12, tan 1, 0 0 0 1 1 1 21, 22, tan 2, 0 0 0ˆ [ ] [ ] ( (E ) H )J Mˆ [ ] [ ] ( (H ) E )B B B inc inc n N N n N N l l n B B B B B B inc inc n N N n N N l l n B B BZ Z n Z Z n − − − −  = = = − − −  = = =   −       =    +            (14) DISCRETIZATION AND ITERATIVE SOLUTION In order to discretize the surface integral equation in ( 14), one popular choice is the Rao -Wilton - Glisson (RWG) basis function [13], which discretizes the surface of the dielectric object into small elements each of which is triangle. Figure 3 depicts a n RWG f unction on a pair of triangles. Figure 4: An RWG function An RWG function associated with the thi edge can be defined as ,2 f ( ) ,2 0,ii i i ii ii ilrT lr r T otherwise + + + − − −  = −     (15) Where il is the length of the main edge, i is area of triangle iT associated with the edge, and i+ and i+ are respectively position vectors with respect to the free vertex of iT+ and iT− . Also note that in figure 3, t he total charg e is zero for each RWG function, ensuring that all edges of iT+ and iT− are free of line charges. Discretization of ( 14) using RWG functions leads to the following:  n n n n ,n ,n n n n n1TT + NK NT TK E J H MTK - NT TT + NKn n  −  =  (16) In which discretized operators can be simplified to the following equations after solving the interaction between the basis and testing ( mt) functions as , , , ,1[ , , , ] ,, , , , ,G ( , ) G ( , )4TT p r q s p r q sp q r s n p q n p q p r n q s nS S S Sn p r q s p r q sjk l l jk l ldr r r dr dr r r dr− +−    =  −        (17) , , , ,1[ , , , ] , , , , , , ,ˆˆ G ( , ) G ( , )42NT p r q s p r q sp q r s n p q n p q p r n q s p r nS S S Sn p r q s p r q sjk l l jk l ln dr r r dr n dr r r dr  − + − +    =   −          (18) ,,[ , , , ] ,, ,,ˆ G ( , )4NK p r q sp q r s pq p r q s nSSn p r q slln dr r r dr++ =     (19) ,,[ , , , ] ,, ,,G ( , )4TK p r q sp q r s pq p r q s nSSn p r q slldr r r dr++ =    (20) And ,n m m1ˆ E t E ( ) t H ( )inc inc nnSS nr dr n r dr=  −   (21) ,n m mˆ H t H ( ) t E ( )inc inc n n nSSr dr n r dr=  +   (22) In the far -field scattering electric field can be written as: 00 1 0 0 E J( ) M( )4jk r jk r sca SSjr e dr r e dr−  −  −     =− (23) The discretization gives: 00 0 0J J 1 M M 0 0 11 2 J 0 , 11 2 M 0 , 11E ( , ) ( ) ( )4 42 42NN jk r jk r sca n n n nSSnn N jk r n n n qSnq nq N jk r n n n qSnq nqjc v r e dr c v r e dr jlc e dr jk lc e dr    −  −  − == −+ == −+ ==    =−  =   −   (24) In which the electric (J) and magnetic (M) currents obtain from ( 16). EXAMPLES In order to show the accuracy and applicability of the proposed surface integral equation (SIE), we have considered in this section two examples of all -dielectric resonator antenna s found in the literature. The first example is the simulation of a P-shape wideband single layer dielectric antenna which is made of a dielectric material with the relative permittivity of 30 [14]. To illuminate the case we used plane wave excitation in the ˆ ik direction with the electric field polarized in the ˆe direction 00ˆ ˆ E (r) .inc i eE exp(ik k r)= (25) Where 0E is the amplitude of the plane wave. 01ˆ H (r) E (r)inc inck= (26) (a) Electric field distribution at 5.4GHz (b) Magnetic field distribution at 5.4GHz (c) Tangential electric field distribution at 5.4GHz (d) Tangential magnetic field distribution at 5.4GHz (e) Electric field distribution at 7GHz (f) Magnetic field distribution at 7GHz (g) Tangential electric field distribution at 7GHz (h) Tangential magnetic field distribution at 7GHz Figure 5: Simulation results of the P -shape dielectric antenna Figure 5 presents the total electric field intensity (V/m) and magnetic field intensity (A/m) on the surface of the P-shape dielectric antenna when it is illuminated by a unit plane wave (1 V/m) in the +z direction. As a second example we have considered the analysis of a concentric half -split dielectric resonator antenna (CDRA) which is developed using permittivity variation in radial direction. In this paper we have considered two CDRAs which are proposed in [15]. Figure 6: Fabricated two and three layer concentric dielectric antenna s [15] The fabricated two -layer and three -layer concentric composite dielectric resonator antennas are shown in figure. The dielectric constant and radius of the two -layer antenna are 1 2 1 210.2, 2.32, 11 , 20 and 11.4rr r mm r mm h mm = = = = = . The optimal design parameters of the three -layer antenna are: 1 2 3 1 210.2, 6.15, 2.32, 5 , 10 ,r r r r mm r mm   = = = = = 320 and 11.4r mm h mm== . ɛr1 ɛr2 ɛr2 ɛr1 ɛr3 (a) Electric field distribution for two -layer DA at 6.92GHz (b) Magnetic field distribution for two -layer DA at 6.92GHz (c) Tangential electric field distribution for two -layer DA at 6.92GHz (d) Tangential magnetic field distribution for two -layer DA at 6.92GHz (e) Electric field distribution for three -layer DA at 5.97GHz (f) Magnetic field distribution for three -layer DA at 5.97GHz (g) Tangential e lectric field distribution for three -layer DA at 5.97GHz (h) Tangential magnetic field distribution for three - layer DA at 5.97GHz Figure 7: Simulation results of the two and three -layer concentric dielectric antenna The electric and magnetic field distribution on the surface of the proposed two and three -layer CDRA have been shown in Figure 7. The antenna is placed in the x -y direction and illuminated by a plane wave (1 V/m) in the +z direction. Computational results in figures 5 and 7 agree well with the experimental discoveries repo rted in the literature. CONCLUSION A surface integral equation formulation obtained from the Maxwell’s equations has been applied to the electromagnetic scattering and radiation problems of three -dimensional (3D) , arbitrary shaped single and multi -layer dielectric antennas. The obtained inte gral equations were discretized with the well -known method of moments (MoM) , where the multilevel fast multipole method (MLFMM) has been employed to speed up matrix -vector multiplication and reduce memory requirements . Some 3D examples were presented to co nfirm the validity and versatility of this approach on dealing with dielectric antennas .