Solution of Volume Integral and Hydrodynamic Equations to Analyze Electromagnetic Scattering from Composite Nanostructures

A coupled system of volume integral and hydrodynamic equations is solved to analyze electromagnetic scattering from nanostructures consisting of metallic and dielectric parts. In the metallic part, the hydrodynamic equation relates the free electron polarization current to the electric flux and effectively"updates"the constitutive relation to enable the modeling of nonlocality. In the metallic and the dielectric parts, the volume integral equation relates the electric flux and the free electron polarization current to the scattered electric field. Unknown electric flux and free electron polarization current are expanded using Schaubert-Wilton-Glisson basis functions. Inserting these expansions into the coupled system of the volume integral and hydrodynamic equations and using Galerkin testing yield a matrix system in unknown expansion coefficients. An efficient two-level iterative solver is proposed to solve this matrix system. This approach"inverts"the discretized hydrodynamic equation for the coefficients of the free electron polarization current and substitutes the result in the discretized volume integral equation. Outer iterations solve this reduced matrix system while the inner iterations invert the discretized hydrodynamic equation at every iteration of the outer iterations. Numerical experiments are carried out to demonstrate the accuracy, the efficiency, and the applicability of the proposed method.


I. INTRODUCTION
I N RECENT years, with the dramatic advances in fabrication technologies, the use of plasmonic nanostructures to manipulate high-frequency electromagnetic fields has become more prevalent than ever [1], [2]. Often, metals are used as the building blocks of these nanostructures since they support surface plasmon modes at optical frequencies. These modes localize the electromagnetic fields in the proximity of the nanostructure and significantly enhance those scattered from it in the far-field region. This has enabled the use of metallic nanostructures as nanoantennas [3], resonators [4], waveguides [5], couplers [6], and sensors [7].
Depending on the frequency, interaction of electromagnetic fields with metals can be accounted for using various models and equations under certain assumptions and approximations. At microwave frequencies, free electrons in a metal have high mobility, which leads to large conductivity and small skin depth (compared with the size of the structure) [8]. Therefore, an electric current, which is confined to the surface of the metal, is used to represent the electromagnetic field interactions on the metal. At optical frequencies, the free electron mobility decreases. As a result, there is a time delay in the response of the electrons to the electromagnetic excitation [9]. To account for this frequency dispersion effect, the classical Drude model [10] is used to represent the permittivity of the metal. Furthermore, at the optical frequencies, the skin depth is usually comparable to the size of a typical nanostructure and the metals are modeled as "penetrable" materials and volume electric currents are used to represent the electromagnetic field interactions on or inside them.
When the frequency is further increased into the ultraviolet part of the spectrum, spatial dispersion appears in the response of the free electrons to the electromagnetic excitation. Effectively, the permittivity becomes nonlocal, i.e., it depends not only on the observation point but also on the source point in space [11]. This spatial dispersion effect is due to the fact that at this frequency regime, a free electron exhibits quantum behavior [12], and the interactions of the electromagnetic fields with the electrons should ideally be modeled using full quantum mechanics simulation methods (e.g., density functional theory [13]). However, the computational cost of This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ these methods is very high, and, therefore, they can only be used when the structure is very small, i.e., only a few nanometers in size [14]. This makes them unsuitable for full-scale simulations of plasmonic nanostructures in real-life scenarios.
This bottleneck can be addressed using a hydrodynamic equation to model the mechanical motion of the free electrons [11]. This equation assumes that the electrons can collectively be accounted for as a moving fluid of charges driven by the electromagnetic fields in the medium. Naturally, these moving charges (termed as the free electron polarization current in the rest of the text) generate electromagnetic fields. This interaction between the free electrons and the electromagnetic fields can be described by a coupled system of the Maxwell and hydrodynamic equations [15], [16], [17], [18], [19]. This system of coupled equations can account for the nonlocality/spatial dispersion and its numerical solution is not as costly as that of the full quantum mechanics simulation methods.
Several methods have been developed to numerically solve the coupled system of the Maxwell and hydrodynamic equations [15], [16], [17], [18], [19]. A majority of these methods are differential equation solvers, e.g., finite element method [16], finite-difference time-domain method [17], and discontinuous Galerkin method [18], [19]. These solvers, just like their traditional versions, which are developed to solve only the Maxwell equations, suffer from several well-known shortcomings that might limit their accuracy and efficiency (see, for example, [20] for details).
Surface integral equation solvers do not suffer from these shortcomings (see, for example, again [20] for details), and indeed, they have been extended to analyze scattering from metallic objects in which the motion of electrons is described by the hydrodynamic equation [15]. However, this method requires the derivation of a new Green function for every type of boundary condition and their combination enforced by the hydrodynamic equation (for example, boundary conditions for the free electron polarization current on metal-metal or metal-dielectric interfaces are different). In addition, this solver is applicable only when the scatterer has homogeneous or piecewise homogeneous material properties.
In this work, these shortcomings are avoided by switching to a volume integral equation formulation. The proposed scheme represents the scattered electromagnetic field in the form of a (volumetric) convolution between the background medium's Green function and the electric flux (induced in the metallic and dielectric parts of the scatterer) and the free electron polarization current (induced in the metallic part). In both the metallic and dielectric parts, the volume integral equation, which relates the electric flux and the free electron polarization current to the scattered electric field, is enforced. In the metallic part, the hydrodynamic equation, which relates the free electron polarization current to the electric flux, is enforced. To numerically solve this coupled system of equations, first the scatterer is discretized into a mesh of tetrahedrons. The electric flux and the free electron polarization current are expanded using a combination of "full" and "half" Schaubert-Wilton-Glisson (SWG) basis functions [21] defined on these tetrahedrons. Inserting these expansions into the coupled system of the volume integral and the hydrodynamic equations and applying Galerkin testing yield a matrix system in unknown expansion coefficients. The boundary condition for the normal component of the free electron polarization current (which should be set to zero on a metal-dielectric interface) is enforced by excluding the half SWG functions, which are defined on the tetrahedrons that have a face on the metal-dielectric interface, from the expansion of the free electron polarization current.
An efficient two-level iterative solver is developed to solve the matrix system resulting from this discretization. This solver "inverts" the discretized hydrodynamic equation for the coefficients of the free electron polarization current and substitutes the result in the discretized volume integral equation. Outer iterations solve this reduced matrix system, while the inner iterations invert the discretized hydrodynamic equation at every iteration of the outer iterations. Note that a preliminary version of the method proposed in this work has been described in [22] as a conference contribution.
The remainder of this article is organized as follows. Sections II-A-II-C describe the formulation of the coupled system of the volume integral and the hydrodynamic equations, its discretization, and the two-level iterative solution, respectively. Section II-D explains how the mesh element size is selected and provides several comments on possible extensions of the proposed method and its applications. Section III presents the numerical results that demonstrate the accuracy, the efficiency, and the applicability of the proposed scheme. Finally, Section IV summarizes this work and outlines future research directions.

A. Coupled System of the Volume Integral and the Hydrodynamic Equations
Let V D represent a composite scatterer that consists of dielectric and metallic parts (Fig. 1). V diel and V H represent these parts, respectively. Both V diel and V H are non-magnetic but their relative permittivity (for bound electron polarization) as denoted by ε b (r) can be inhomogeneous. The boundary surface of V H is represented by S H . The scatterer resides in an unbounded homogeneous background medium with permittivity ε 0 and permeability μ 0 . The electric field of the excitation is denoted by E inc (r) and its frequency is denoted by ω.
Upon this excitation, equivalent volumetric electric current J(r) is induced in V D and this current generates the scattered electric field E sca (r). The incident electric field E inc (r), the scattered electric field E sca (r), and the total electric field E(r) satisfy the fundamental field relationship. (1) Using the volume equivalence principle [20], J(r) is expressed in terms of E(r) and the electric flux D(r) as follows: and E sca (r) is expressed as follows: Here, L V [X](r) is the volume integral operator given by the following equation: is the Green function, k 0 = ω √ μ 0 ε 0 is the wavenumber in the background medium, and R = r − r denotes the distance between the observation point r and the source point r .
The motion of the free electrons in V H is accounted for using a nonlocal hydrodynamic equation. Let J H (r) represent the polarization current associated with these free electrons. Then, J H (r) is expressed in terms of E(r) and D(r) as [15] In addition, J H (r) is "driven" by E(r). This relationship is described by the nonlocal hydrodynamic Drude equation as [11] Here, ω p is the plasma frequency, β 2 = 0.6v 2 F , v F is the Fermi velocity, and γ is a damping constant. The hydrodynamic (5) has to be complemented by a boundary condition on the boundary surface of V H [as denoted by S H (Fig. 1)] [11] n(r) · J H (r) = 0, r ∈ S H .
Here,n(r) is the outward pointing unit normal vector on S H . This boundary condition ensures that the normal component of J H (r) vanishes on S H (i.e., free electrons do not flow from the metallic part into the dielectric part or the background medium). Note that the hydrodynamic (5) together with the boundary condition (6) introduces an electromagnetic wave solution with a longitudinal electric field in V H in addition to the one with a transverse electric field [11]. The proposed scheme defines D(r) and J H (r) as unknowns. To eliminate E(r), E(r) from (4) is inserted into (2). This yields an expression for J(r) in terms of D(r) and J H (r) as follows: where κ(r) = 1 − 1/ε b (r). Substituting (7) into (3) and inserting the resulting equation and E(r) from (4) into (1) yield the volume integral equation in unknowns D(r) and J H (r) as follows: Similarly, the hydrodynamic (5) should be expressed in only D(r) and J H (r). Inserting E(r) from (4) into (5) yields Equations (8) and (9) are the final form of the coupled system of the volume integral and the hydrodynamic equations in unknowns D(r) and J H (r). This system is discretized using the scheme described in Section II-B. Note that this discretization scheme ensures that the boundary condition (6) is enforced correctly.

B. Discretization
To numerically solve the coupled system of (8) and (9), first, V D is discretized into a mesh of tetrahedrons. Then, the unknowns D(r) and J H (r) are expanded as follows: Here, {I D } n and {I H } n are the unknown coeffcients, and f D n (r) and f H n (r) are the basis functions constructed using the SWG functions defined on the triangles of the tetrahedral mesh [21]. The SWG basis function associated with triangle n is defined as follows: Here, V + n and V − n are the tetrahedrons "touching" triangle n on its two sides, r ± n are the corners of V ± n that are not on S n (i.e., free nodes), |S n | is the area of S n , and |V ± n | are the volumes of V ± n . The basis set f D n (r) includes "full" SWG basis functions as defined by (11)  The basis set f H n (r) consists of only the full basis functions as defined by (11) on every pair of tetrahedrons in V H and does not include the half SWG basis functions defined in single tetrahedrons that have their S n on S H (metal-dielectric interface). Note that J H (r) does not flow from the metallic part into the dielectric part or the background medium and its normal component on S H is zero as described by the boundary condition in (6). Exclusion of the half SWG basis functions from f H n (r) ensures that this boundary condition is correctly enforced. N H in (10) is the total number of full SWG basis functions used in the expansion of J H (r) in V H .
Inserting the expansions (10) into (8) and (9) and Galerkin testing the resulting equations using In (12), the entries of the block matrices Z DD of dimension and the tested incident field vector V inc of dimension N D are given by the following equations: respectively. Here, the inner product between vector functions a(r) and b(r) is defined as follows: where V a is the support of a(r).
are "global" operators, one can see from (13) and (14) that Z DD and Z DH are dense matrix blocks. On the other hand, since f D n (r) and f H n (r) have "local" supports (two tetrahedrons for a full SWG function and one tetrahedron for a half SWG function), one can see from (15) and (16) that Z HD and Z HH are sparse matrix blocks (the maximum number of nonzero entries in one row of these blocks is seven). Note that the detailed expressions for the matrix entries in (13)- (16) and the vector entries in (17) are provided in the Appendix.

C. Solution of the Matrix equation
To take advantage of the sparsity of Z HD and Z HH directly, the coupled matrix system (12) is solved iteratively for the unknown coefficient vectors I D and I H . Two approaches can be used for this purpose.
1) Single-Level Iterative Solver: The coupled system (12) is iteratively solved as a whole using a transpose-free quasiminimal residual (TFQMR) scheme [23]. The computational cost of matrix-vector multiplication ZĨ required at every iteration of TFQMR scales as O( The four terms in this expression represent the computational cost of multiplying matrix blocks Z DD , Z DH , Z HD , and Z HH with the relevant part ofĨ, respectively. Then, the overall computational cost of this single-level iterative solver scales as follows: (19) where N it is the number of iterations required for the relative residual error to converge to a user-defined value.
2) Two-Level Iterative Solver: In this approach, before an iterative solver is used, the coupled system (12) is first reduced into a smaller matrix system. This is done by inverting the second row of (12) for I H , i.e., I H = −Z −1 HH Z HD I D , and inserting this expression into the first row. This yields a smaller matrix system of dimension N D × N D in unknown I D as follows: Then, TFQMR is used to solve (20) for I D . Matrix-vector multiplication (Z DD − Z DH Z −1 HH Z HD )Ĩ D required at every iteration of TFQMR is carried out as described below.
Step 1: Compute the first term Z DDĨD .
Step 2: Compute the second term Z DH Z −1 HH Z HDĨD in three steps as Step 2.1: Compute y = Z HDĨD .
Step 2.2: Compute x = Z −1 HH y by solving y = Z HH x for x. This is done iteratively using TFQMR.
Step 2.3: Compute Z DH x.
it is the number of iterations required for the relative residual error to converge to a user defined value during the solution of the matrix equation y = Z HH x at Step 2.2 (i.e., inner iterations). Then, the overall computational cost of this two-level iterative solver scales as follows: where N out it is the number of iterations required for the relative residual error to converge to a user defined value during the solution of the matrix (20) (i.e., outer iterations).
Comparing (21) to (19), one can see that the single-level iterative solver would certainly be faster than the two-level iterative solver for N it ≤ N out it . However, numerical results presented in Section III-A show that N out it is much smaller than N it and N in it is small, and therefore the two-level iterative solver is significantly faster than the single-level iterative solver. The difference between N it and N out it can be explained by the fact that the dimension of (20) is almost half of that of (12) and inserting I H = −Z −1 HH Z HD I D into the first row of (12) [to obtain (20)] effectively takes care of the scaling difference between the volume integral and the hydrodynamic equations.
Note that for both the approaches, the computational cost of matrix-vector multiplications Z DDĨD and Z DHĨH can be reduced using the fast multipole method [24], [25], [26], [27] and its multilevel versions [28], [29], [30], [31] as well as other matrix compression schemes like those described and referred to in [32], [33], [34], and [35]. But this does not change the conclusions of the above comparison since the difference in the computational cost of the two iterative solvers is mainly due to the difference in the number of iterations.

D. Comments
A metallic medium that is described by the hydrodynamic (5) supports the propagation of electromagnetic waves with electric fields along the transverse and longitudinal directions [11]. Let the (complex) wavenumbers associated with these electromagnetic waves be denoted by k T (r) and k L (r), respectively. The expressions of k T (r) and k L (r) are given by [11] For example, for gold, ε b (r) = 1, ω p = 1.20 × 10 16 s −1 , β = 1.07 × 10 6 m/s, and γ = 1.36 × 10 14 s −1 [36], [37]. Fig. 2 shows the plots of the values of k T (r) and k L (r) computed for gold in the frequency range ω ∈ [0.5ω p , 1.5ω p ].
Note that during the computation of square roots in (22), the positive imaginary part of the result is selected to avoid a nonphysical growing wave. The figure clearly shows that in this frequency range, both the real and imaginary parts of k L (r) are significantly larger than those of k T (r), respectively. This means that to accurately capture the behavior of the electromagnetic fields in V H (inside the metallic part), the mesh of tetrahedrons must resolve the wavelength associated with k L (r). Note that within the frequency range considered here, since k 0 is significantly smaller than both k L (r) and k T (r), the mesh in V diel = V D − V H (inside the dielectric part) can ideally be coarser than the one in V H . But since these two volumes share a surface and a conformal discretization is used, the mesh in V diel is denser than what it would ideally be. This unnecessary computational overhead can be alleviated by switching to a nonconformal discretization scheme, such as the one described in [38] and [39]. Development of such a scheme is underway. Several comments about the possible extensions of the proposed method are in order. 1) The formulation in Section II-A assumes that the scatterer includes only one type of homogeneous metalic structure (described by a set of hydrodynamic equation parameters γ , β, and ω p ). The formulation can easily account for additional metal types described by assigning different values to γ , β, and ω p in relevant matrix entries as long as the structures made of different metals do not touch each other. This limitation stems from the fact that the boundary condition (6) is not valid on metal-metal interfaces and therefore the discretization scheme described in Section II-B is not applicable anymore. Development of a formulation that removes this limitation and allows for modeling of metal-metal interfaces is currently underway. 2) In [40], the nonlocal hydrodynamic Drude model [as described mathematically in (5)] has been extended to account for the classical kinetic effects of the charge carrier diffusion in the nonlocal response. This so-called generalized nonlocal optical response (GNOR) model expresses the relationship between J H (r) and E(r) as follows: Here, β, γ , and ω p are the same as those in (5) and the additional parameter D is the charge carrier diffusion constant. Comparing (5) and (23), one can easily see that the only difference between the two equations is the additional term D(γ + j ω) in (23). Furthermore, both the nonlocal models use the same boundary condition given in (6). Therefore, the implementation of the GNOR model within the numerical scheme proposed here is rather trivial and can simply be done by replacing β 2 by β 2 + D(γ + j ω).
3) In [41] and [42], an analytical method that relies on the Mie series expansion of the fields has been used to study the nonlocal response of nanospheres in three different scenarios: electromagnetic scattering, electron energy-loss spectroscopy, and atomic spontaneous emission.
Even though examples presented in Section III involve only electromagnetic scattering problems under plane-wave excitation, the proposed numerical scheme is applicable to other problems with different types of excitation, including electron energy-loss and atomic spontaneous emission.
III. NUMERICAL RESULTS In this Section, several numerical examples are presented to demonstrate the accuracy, the efficiency, and the applicability of the proposed solver. All the scatterers considered in these examples reside in free space with permittivity ε 0 and permeability μ 0 . For all the examples, the excitation is a plane wave with electric field E inc (r) =pE 0 e − jk 0k inc ·r .
(24) Here, the unit vectorp represents the direction of the polarization, E 0 is the electric field amplitude, and the unit vector k inc represents the direction of propagation. Unless otherwise stated,p =x, E 0 = 1 V/m, andk inc =ẑ for all the examples considered here. All the TFQMR iterations (outer and inner iterations for the two-level iterative solver and the iterations for the single-level iterative solver) are terminated when the relative residual error (RRE) reaches a desired level, i.e., when the condition ||b − AI n ||/||b|| ≤ χ RRE is satisfied. Here, I n is the solution at iteration step n, A is the matrix, b is the righthand side vector, and χ RRE is the convergence threshold.
Two sets of simulations are carried out. For the first set of simulations, frequency is set to ω = 0.5 ω p and three levels of mesh are used. These meshes use N D = {5494, 9216, 17 546} and N H = {5030, 8612, 16 694} basis functions to discretize D(r) and J H (r) induced inside the sphere, respectively. The single-level and two-level iterative schemes are used to solve the matrix system (12), which is preconditioned from left using a diagonal preconditioner, and the matrix system (20), respectively. For the iterations of the single-level scheme and the outer iterations of the two-level scheme, χ RRE = 10 −4 . For the inner iterations of the two-level scheme, χ RRE = 10 −8 . Note that the accuracy of the inner iterations has to be high to ensure that the outer iterations converge. This does not increase the computational cost significantly because the converge rate of the inner iterations is already very fast. Fig. 3(b) shows plots of the radar cross section (RCS) computed on the x y-plane from the solutions obtained by the single-level and two-level iterative schemes for the mesh with N D = 17 546 and N H = 16 694 and the analytical Mie series with nonlocal response material model [43], [44]. The results agree very well. Table I shows comparison of the performance of the single-level and two-level iterative solvers. The two-level iterative solver is significantly faster due to fact that N out it is much smaller than N it and N in it is small (see the computational complexity comparison in Section II-C). Note that the numbers for N in it presented in Table I are the range of the inner iterations.

TABLE I PERFORMANCE OF THE SINGLE-LEVEL AND TWO-LEVEL ITERATIVE SOLVERS IN ANALYZING SCATTERING FROM A GOLD NANOSPHERE
For the second set of simulations, a total of 60 simulations are carried out using the two-level iterative solver at equally spaced points in the frequency range ω ∈ [0.5ω p , 1.5ω p ]. In these simulations, D(r) and J H (r) induced inside the sphere are discretized using N D = 143 254 and N H = 140 022 basis functions, respectively. For the outer and inner iterations of the two-level iterative scheme, χ RRE = 10 −4 and χ RRE = 10 −8 , respectively. Fig. 3(c) shows plots of the extinction cross section (ECS) computed from the solutions obtained using the two-level iterative scheme and the analytical Mie series with nonlocal and local response material models versus ω/ω p [43], [44]. The figure clearly shows that the result obtained using the proposed method matches well with the result obtained using the Mie series with the nonlocal response material model. The figure also shows that the error in the result obtained by the proposed method increases with ω. This error can be reduced using a denser mesh that can capture the behavior of the fields with large k L (r) more accurately (see Section II-D and Fig. 2).
The first peak observed in all three ECS curves is caused by the transverse field resonance. Also, a "blue shift" phenomenon is shown in this figure, i.e., the transverse field resonance peak shifts toward higher frequencies when the nonlocal material response is taken into account [40]. Furthermore, three other peaks are identified at ω = 1.13 ω p , ω = 1.29 ω p , and ω = 1.50 ω p in ECS computed by the proposed solver and the Mie series with the nonlocal response material model. These are caused by the longitudinal field resonance. It also can be concluded from Fig. 3(c) that the transverse field response is more dominant at ω < ω p while the longitudinal field response is more dominant at ω > ω p .

B. Metallic Dimer
In this example, electromagnetic scattering from a nanodimer [ Fig. 4(a)] is analyzed using the proposed solver. The radius of the spheres is 1 nm, and the shortest distance between them is 0.2 nm. The hydrodynamic equation parameters for the material making up the spheres are ω p = 1.20 × 10 16 s −1 , γ = 1.36 × 10 14 s −1 , v F = 1.39 × 10 6 m/s, and ε b (r) = 1. Note that for this problem, V D = V H and V diel = ∅.
A total of 60 simulations are carried out using the twolevel iterative scheme at equally spaced points in the frequency range ω ∈ [0.5ω p , 1.  Fig. 4(b) shows plots of the scattering cross section (SCS) computed from the solutions obtained using the two-level iterative scheme versus ω/ω p [45]. The transverse field resonance peaks are observed at ω = 0.6 ω p and ω = 0.75 ω p , while the longitudinal field resonance peaks are observed at ω = 1.13 ω p , ω = 1.29 ω p , and ω = 1.50 ω p . Furthermore, the figure shows that the nanodimer supports bonding and antibonding modes (generated by the transverse field resonances) at ω = 0.72 ω p and ω = 0.75 ω p , respectively [15].

C. Composite Sphere
In this example, electromagnetic scattering from a silicacoated gold nanosphere [ Fig. 5(a)] is analyzed using the proposed solver. The radius of the gold sphere is 1 nm, and the thickness of the silica coating is 1 nm. The hydrodynamic equation parameters for gold are the same as those in Section III-A, and the relative permittivity of silica is ε b (r) = 2.25. Note that for this problem, A total of 60 simulations are carried out using the two-level iterative scheme at equally spaced points in the frequency range ω ∈ [0.4ω p , 1.4ω p ]. In these simulations, D(r) induced inside the coating and the sphere and J H (r) induced inside the sphere are discretized using N D = 97 642 and N H = 80 486 basis functions, respectively. For the outer and inner iterations of the two-level iterative scheme, χ RRE = 10 −4 and χ RRE = 10 −8 , respectively. Fig. 5(b) shows plots of ECS computed from the solutions obtained by the two-level iterative scheme and the Mie series with nonlocal response material model for the sphere and local response material model for the coating versus ω/ω p . The results agree well; however, as expected, the error in the result obtained by the proposed method increases with ω. This error can be reduced using a denser mesh that can capture the behavior of the fields with large k L (r) more accurately (see Section II-D and Fig. 2).
The transverse field resonance peak is observed at ω = 0.55 ω p and the longitudinal field resonance peaks are observed at ω = 1.13 ω p and ω = 1.31 ω p . Comparing Figs. 3(c) and 5(b), one can see that the transverse field resonance peak shifts toward lower frequencies due to the presence of silica coating.

D. Metallic Cylinder
In this example, electromagnetic scattering from a gold nanocylinder [ Fig. 6(a)] is analyzed using the proposed method. Note that the ends of the cylinder are rounded here since sharp edges cannot be often fabricated at nanoscales. Four nanocylinders with length H ∈ {2.5 nm, 3 nm, 4 nm, and5 nm} are considered. All the four cylinders have a radius of 1 nm. The hydrodynamic equation parameters for gold are the same as those in Section III-A. Note that for this problem, V D = V H and V diel = ∅.
A total of 60 simulations are carried out using the two-level iterative scheme at equally spaced points in the frequency range ω ∈ [0.5ω p , 1.25ω p ]. Four levels of mesh are used for the four nanocylinders and these meshes use N D = {58 312, 60 948, 61 010, 62 186} and N H = {56 384, 58 900, 58 962, 60 006} basis functions to discretize D(r) and J H (r) induced inside the cylinder, respectively. For the outer and inner iterations of the two-level iterative scheme, χ RRE = 10 −4 and χ RRE = 10 −8 , respectively. Fig. 6(b) show compares of ECS computed from the solutions obtained using the two-level iterative scheme for all the four nanocylinders versus ω/ω p . As seen from the figure, the transverse field resonance peak shifts from ω = 0.67 ω p to ω = 0.69 ω p when the length of the cylinder is increased from 2.5 to 5 nm. The longitudinal field resonance peaks are observed at ω = 1.1 ω p and ω = 1.25 ω p . In addition, this figure shows an extra resonance peak that shifts from ω = 0.81 ω p to ω = 0.73 ω p with increasing length. Note that this resonance (associated with the length of the cylinder along theŷ-direction) is induced even though the incident electric field E inc (r) is polarized in thex-direction.

E. Metallic Cylinder on Top of a Dielectric Slab
For the last example, electromagnetic scattering from a gold nanocylinder located on top of a silica substrate [ Fig. 7(a)] is analyzed using the proposed method. The length and the radius of the cylinder are 4 and 1 nm, respectively. The width, the length, and the height of the slab are 6, 6, and 1 nm, respectively. The shortest distance between the cylinder and the slab is 2 nm. The hydrodynamic equation parameters for gold are the same as those in Section III-A, and the relative electric permittivity of silica is ε b (r) = 2.  N H = 0). The simulations of these three scenarios are carried out using the two-level iterative scheme [reduces to a "traditional" volume integral equation solver [46], [47] with single-level iterations for scenario (iii)] at 60 equally spaced points in the frequency range ω ∈ [0.5ω p , 1.25ω p ]. Note that for scenarios (ii) and (iii), only 284 basis functions are used to discretize D(r) induced inside the slab (see the discussion in Section II-D). For the outer and inner iterations of the two-level iterative scheme, χ RRE = 10 −4 and χ RRE = 10 −8 , respectively.  Fig. 7(b) shows plots of ECS computed in the three scenarios described above versus ω/ω p . As expected, ECS of the silica slab [scenario (iii)] does not involve any resonance peaks since silica does not have any plasmonic properties. Second, the slab is located far away from the gold nanocylinder, so the coupling between them is not expected to be strong, and ECS in scenario (i) and ECS in scenario (ii) are close to each other especially at the lower end of the frequency range. The difference increases at the higher end, which might be explained by the fact that ECS of the silica slab is larger at higher frequencies [scenario (iii)] and therefore contributes more to the total ECS in scenario (ii).

IV. CONCLUSION
Electromagnetic scattering from nanostructures consisting of the metallic and dielectric parts is analyzed by solving a coupled system of the volume integral and hydrodynamic equations. The hydrodynamic equation, which is enforced only in the metallic part, relates the free electron polarization current to the electric flux. This equation effectively updates the constitutive relationship and permits modeling of the nonlocality. The volume integral equation, which is enforced in both the metallic and dielectric parts, relates the electric flux to the scattered electric field. Unknown electric flux and free electron polarization current are expanded using a combination of full and half SWG basis functions. The boundary condition associated with the free electron polarization current on the metal-dielectric interface is enforced by excluding the half SWG basis functions, which are located on this interface, from the expansion of the free electron polarization current. Inserting these expansions into the coupled system of the volume integral and hydrodynamic equations and using Galerkin testing yield a matrix system.
An efficient two-level iterative solver is developed to solve this matrix system. This approach inverts the discretized hydrodynamic equation (bottom rows of the matrix system) for the coefficients of the free electron polarization current and substitutes the result in the discretized volume integral equation (top rows of the matrix system). Outer iterations solve this reduced matrix system, while the inner iterations invert the discretized hydrodynamic equation at every iteration of the outer iterations.
Numerical experiments, which involve the computation of RCS, ECS, and SCS for metallic and composite nanostructures, are carried out to demonstrate the accuracy, the efficiency, and the applicability of the proposed method.
Future research directions include incorporation of different nonlocal hydrodynamic equations, implementation of boundary conditions on metal-metal interfaces, acceleration of the matrix solution using matrix compression schemes, and application of the proposed solver to different problems with different types of excitations.

APPENDIX ENTRIES OF THE MATRIX AND THE RIGHT-HAND SIDE VECTOR IN (12)
While computing the entries of Z DD , Z DH , Z HD , and Z HH , it is assumed that ε b (r) and κ(r) are constant in a given tetrahedron, and the values of these constants are obtained by sampling ε b (r) and κ(r) at the center of that tetrahedron. Therefore, one can define local functions as follows:

A. Entries of Z DD
If both f D m (r) and f D n (r) are full SWG functions, then