Mie Scattering with 3D Angular Spectrum Method

Mie theory is a powerful method to model electromagnetic scattering from a multilayered sphere. Usually, the incident beam is expanded to its vector spherical harmonic representation defined by beam shape coefficients, and the multilayer sphere scattering is obtained by the T-matrix method. However, obtaining the beam shape coefficients for arbitrarily shaped incident beams has limitations on source locations and requires different methods when the incident beam is defined inside or outside the computational domain or at the scatterer surface. We propose a 3D angular spectrum method for defining beam shape coefficients from arbitrary source field distributions. This method enables the placement of the sources freely within the computational domain without singularities, allowing flexibility in beam design. We demonstrate incident field synthesis and spherical scattering by comparing morphology-dependent resonances to known values, achieving excellent matching and high accuracy. Additionally, we present mathematical proof to support our proposal. The proposed method has significant benefits for optical systems and inverse beam design. It allows for the analysis of electromagnetic forward/backward propagation between optical elements and spherical targets using a single method. It is also valuable for optical force beam design and analysis.


Introduction
Electromagnetic scattering from a homogeneous sphere illuminated by an arbitrary incident beam can be computed with conventional methods, including full-wave simulations [1], geometrical optics [2], or physical optics [3,4].These techniques are well-studied and accurate, given the model fidelity and a suitable wavelength range.However, without considerable computational effort, they cannot assess the internal and scattered electric fields from multilayered spherical objects.Especially when the sphere's radius is of the order of wavelengths, complete classical electromagnetic wave theory is needed [5].
Mie theory and the generalized Lorentz-Mie theory are accurate methods to evaluate the internal and scattered fields from the multilayered spheres [6][7][8], where the incident beam is presented in vector spherical harmonics (VSH) expansion defined by beam shape coefficients (BSCs) [9].However, obtaining the BSCs for arbitrarily shaped incident beams can be difficult; often needing a combination of several complex methods [6,10].The BSCs can be computed from a known function or field distribution with certain constraints, such as polarization, source shape, and location limitations [11,12].In these cases, BSCs for the incident beam are obtained by the Bromwich method [13] or multipole expansion [14][15][16].These methods define the sources inside a closed volume, and the VSH expanded fields can be computed only outside of this area, limiting the source location.
BSCs can also be computed from known electric and magnetic field distributions on the surface of the spherical scatterer using closed surface orthogonality (CSO).CSO also has limitations as the fields must be defined across the entire spherical surface [17].For example, CSO cannot be applied to fields defined only on spherical subregion thus excluding it from many inverse beam synthesis tasks, where the beam only illuminates a small area of the sphere [18].On the other hand, BSCs can be computed from the standard 2D ASM, which does not have these location restrictions due to its eigenfunction property.However, 2D ASM is limited to the planar surface distributions [19,20]; see the comparison of the methods in Table (1).Bromwich [13] Multipole exp.[14][15][16] CSO [17] Sanford et al. [4] 2D ASM [19,20] Stepwise ASM [21,22] 3D ASM This article aims to create a VSH-expanded electromagnetic beam from a known electric field distribution on arbitrarily shaped and positioned surface.Then the incident beam could be modified in terms of propagating field distribution, polarization, and local phase variation.Also, the source location can be positioned inside the simulation area without restrictions.This goal can be achieved by expanding the 2D ASM to a 3D AS approach to accommodate arbitrary fields defined on arbitrarily shaped surfaces and construct modified BSCs for an incident beam VSH presentation.
2D ASMs present the incident field in the angular spectrum domain as the sum of differently oriented plane waves.When this angular summation of plane waves on a source plane is presented in a direction cosine coordinate system, these angles can be used directly to compute the BSCs for the VSH expansion.Additionally, the internal and scattered fields from a multilayered dielectric sphere can be mapped with the extended boundary condition method (EBCM), which is referred as the T-matrix method in this article [20].
In [21], the planar AS method was expanded to approximate diffractive fields from 2D curved surfaces.The expansion was obtained by dividing the curved 2D surface into the step-wise subregions windowed by the Gaussian function.This approach was later expanded to the 3D surfaces in [22], where the authors used the same planar step-wise subsections with Gaussian distribution.Both methods approximate well the diffraction, reflections, and transmission from the curved single-layer boundaries when the radius of curvature (RoC) is much larger than the wavelength [23].
The 3D ASM derivation begins with the same approach of dividing an arbitrary surface into planar subregions, presented as equal amplitude areas without Gaussian distribution or windowing (e.g.Stepwise AS).Then the areas of these subregions are shrunk into infinitesimally small points, which approach a Dirac delta function in the limit.This approach yields exact results (in the limit of infinitesimally small partitions) of the diffractive field from an arbitrary surface considering surface structures comparable to or larger than the illumination wavelength, i.e.,  ≥ 2, where  is the wavenumber and  is the radius of the sphere.
Moreover, this result is analogous to the angular spectrum presentation of the Stratton-Chu method with magnetic dipole sources, where the electromagnetic radiation from a source point is calculated as a curl of a magnetic dipole multiplied by scalar Green's function [24].This presentation satisfies Maxwell's equations and approximates well curved surfaces when the RoC is larger than a wavelength [25].Also, the advantages of the presented method are non-singularity at the source points, allowing the synthesis of a continuous beam through the source surface from any known electric field distribution.Finally, the incident beam can be expanded into a VSH presentation to compute scattered fields from multilayered spheres.

Theory
This chapter presents a method to compute the forward/backscatter of multilayered spheres under arbitrary beam illumination.The chapter is divided into three sections introducing the modified 3D ASM, global coordinate mapping system, and modified BSCs coefficients for the VSH expansion.
Section 2.1 begins by introducing a general parametrization of an arbitrary surface.The local base vectors are derived from parametrization to determine a Cartesian coordinate system for surface sources.Further, the electric field polarization is modified by base vector rotation.Also, a transformation matrix is presented for mapping the local source coordinates in a global coordinate system and vice versa.After coordinate mapping, an arbitrary surface parametrization is discretized into infinitesimally small surface elements.The key concept is that Riemann's surface integral combines differential surface elements, which can be approximated as locally planar elements.The ASM models the radiated field from each element, and the total electromagnetic field is the sum of the modeled fields from each source point by the superposition principle [19].This approach relies on the assumption that when the areas of the elements are small enough, the synthesized field from the differential elements approaches the field from the original surface.
Section 2.2 presents a system to map local coordinates into a global origin-centered coordinate system with each source point's orientation and position information, which is later used to construct the BSCs.
Section 2.3 presents the VSH expansion from an arbitrary surface with modified BSCs.These BSCs are modified to include the mapped source's orientation and position information while preserving the needed spherical symmetry for VSH expansion.

Modified angular spectrum method
Let the global Cartesian base-vectors be (e  , e  , e  ) in space R 3 , where position vector is r = e  + e  + e  and global coordinates are marked as  glob = (, , ).Consider a generalized, compact surface Ω in R 3 , which has a continuously differentiable parametrization with parameters  and  as Let us place the origin of the local coordinate system at point o( , ) on the surface Ω, where two base vectors f 1 and f 2 define tangential plane at the point, and the third base vector f 3 is a local surface normal vector, see Figure (1).When normalized, these vectors form an orthonormal base (f 1 , f 2 , f 3 ) and are defined as The final local base (e 1 , e 2 , e 3 ) is obtained from the base (f 1 , f 2 , f 3 ), when the local polarization of the electric field is chosen.Let us set where the sign is chosen according to the desired electromagnetic propagation direction.In a simple case, when polarization is known on the surface Ω, the local electric field vector e 1 on Ω can be presented with the base (f 1 , f 2 ), and the local magnetic field vector e 2 is obtained as a cross-product with local propagation direction e 3 as e 2 = e 3 × e 1 .
Polarization on the surface Ω can also have externally sourced initial conditions, for example, from a polarization of an incident beam, that creates the electric field distribution on the surface Ω.Let p = p( , ) be a vector that is perpendicular to the plane of that incident beam polarization.Then, local direction e 1 is an intersection of that plane, and the tangential plane of surface Ω spanned by (f 1 , f 2 ).
The vectors (e 1 , e 2 , e 3 form a right-handed orthonormal base, which together with the origin o defines a local coordinate system.Let a vector r be presented in this base as r = xe 1 + ȳe 2 + ze 3 .Its coordinate triple is marked briefly as  loc = ( x, ȳ, z).When the vector r is expressed with both (e  , e  , e  ) and (e 1 , e 2 , e 3 ) bases, as it is known, the dependence of the coordinates on each other is determined by an orthogonal transformation matrix Θ( , ) = e 1,glob e 2,glob e 3,glob , where the columns are the coordinate triples e Let s denote briefly r = r − o.Thus, we have presentations of the location r in both local and global coordinate systems.The dependence of the corresponding coordinates can be written using the transformation matrix Eq. ( 5) as Let us divide the segment of the surface Ω containing the electric field distribution  0 ( , ) into a finite partition of small separate sets ∪  Ω  = Ω.Select the point o( , ) on differential element Ω  , the vectors associated with it are e 1 and e 2 .Area Ω  is projected to the plane spanned by the vectors e 1 and e 2 .The resulting area at the plane is marked as Ω  .Likewise, function  0 is projected into that plane, limited to set Ω  and zero elsewhere.This geometry is presented in Figure 2. The resulting projection is marked as   0 .In the local coordinate system it holds z = 0 on the piece Ω  , because this is located on the x ȳ−plane.  0 is valid when x ≈ 0 and ȳ ≈ 0. The angular spectrum theory is locally applied on the local coordinate system to function   0 valid on piece Ω  .Let's observe the field it creates at point r, where the coordinates are as in Eq. ( 6).At a fixed point r holds [19] where F   0 is the Fourier transform of the local electric field   0 and  = ( 2  +  2  +  2  ) 1/2 is the wavenumber.The total field from the arbitrary surface is obtained by summing the fields from the differential sources by the superposition principle as where local coordinates ( x, ȳ, z) as in Eq. ( 6) and the total field is polarized along the global e 1 unit-vector.The obtained field Eq. ( 9) also satisfies the Helmholtz equation as unit vector e 1 is locally constant.To be precise, if the limit with respect to the areas |Ω  | of the sets Ω  exists, we define As pointed out at the beginning of the Theory section, this limit also works physically when considering smooth surfaces and is even more precise when the wavelength decreases compared to the local radius of the surface.Finally, a further examination is made to obtain a closed form for clause Eq. ( 10); we present a proof in Appendix A.
Heuristically, when |Ω  | shrinks, the function normed in volume, approaches a Dirac delta  (0,0) in the sense of the distribution theory; since  1 ( , ) =   1 ( 0, 0) when o is fixed.Then Eq. ( 9) and Eq. ( 11) yield first the form of a Riemann sum and then an integral as a limit The final form can be written as where and x, ȳ and z are as in Eq. ( 6) and notation (r; , ) defines electric field on observation point r related to the location (, ) on the parametrizied surface.Propagating waves are obtained by integrating Eq. ( 14) over the real disk  2 x +  2 ȳ ≤  2 , and evanescence waves can be obtained by expanding the integration over the real domain, i.e.,  2 x +  2 ȳ >  2 .There are two key observations regarding this derivation.First, Eq. ( 14) is similar to E  = −∇ × (M) after replacing the scalar Green's function () by its angular spectrum given by Weyl identity, where M = −2n × e 2 is a small magnetic dipole [26].Thus, Eq. ( 13) is a particular case of the Stratton-Chu equation for the electric field in which only the magnetic sources M = −n × E   are presented and the electric ones J = n × H   have been removed [27].Secondly, if the surface Ω is planar, Eq. ( 13) and Eq. ( 14) return to the original two-dimensional angular spectrum method.

Polarization and global coordinate mapping
An electromagnetic beam in the plane wave spectrum representation also has polarization vector components along the propagation direction.To account for polarization let us expand Eq. ( 14) along plane (e 1 , e 3 ) using the relation E  • k = 0 [28].Let k =  x +  ȳ +  z be the wave vector, where  = k is the wave number and Hence the summarized exponent in Eq. ( 14) can be written as where e 0 ( x ,  ȳ ,  z ) defines the polarization.The local Cartesian coordinate system presents differential source Eq. ( 15) in the base (e 1 , e 2 , e 3 ).When expanded to the VSH presentation, the influence of each differential source must be presented by means of position vector r to preserve the spherical symmetry of VSH.Next, the global mapping of local coordinates is presented.
We have a simple connection (r − o) glob = Θ rloc as in Eq. ( 6).Let us define s | | = s, when z > 0, but, when z < 0, the sign of z−coordinates in s is exchanged.Consequently, we can write where written briefly t = s| |, .As a result, we get Based on Eq. ( 17), the term   (s•r || ) , which is in the local coordinate system, can be presented in the global coordinate system as where  − (t•o loc ) =  − (s || •o) presents the position and phase shift of the source beam's focus (waist) from global origin to the source surface.Each local beam created from the differential source is presented in a global coordinate system.

VSH expansion with modified beam shape coefficients
Modified BSCs coefficients for each source point are derived for presenting any parameterized field distribution in a global coordinate system, including the locations and orientations of the sources.These coefficients map complicated coordinate geometries into one presentation.This allows synthesizing a total electromagnetic beam as a superposition of VSH expanded fields.The integration wave numbers  x ,  ȳ in Eq. ( 15) are analogous to the BSC angles [20].Thus Eq. ( 15) is expressed in the direction cosine coordinate system as where  and  are direction cosines related to Cartesian coordinates as  cos  =  x and  cos  =  ȳ .The polarization e 0 in terms of direction cosines is given as and the vector s is Then, Eq. ( 19) is evaluated by numerical integration over a disk of radius  using the trapezoidal rule method, which is sufficiently accurate for periodic functions [29], with uniform-width  [20] where indices  and  refer to individual plane waves with propagation constants   x and   ȳ respectively to x and ȳ directions, and  =  +1 −   =  +1 −   is the step size of the numerical integration and t   = (s | |,  ) loc .The integration radius is limited to  ≤ 2/ for obtaining only propagating waves, and the evanescence fields can be obtained by expanding the integration domain.
The transformation matrix from local coordinates ( x, ȳ, z) to global coordinates (, , ) is Θ( , ) as in Eq. ( 5), and the radiated electric field can be written at an arbitrary point r by VSHs as [20] where M 1 emn , M 1 omn , N 1 emn and N 1 omn are the VSH of the first kind and   is a normalization factor [20].For r we have to substitute its spherical coordinates (, , ) and corresponding spherical base vectors (e  , e  ), taken with respect to the local base (e 1 , e 2 , e 3 ).They can be computed from the triple  loc .Modified BSCs   emn ,   omn ,   emn and   omn are given as The BSCs for each  − plane wave   emn   ,   omn   ,   emn   and   omn   are defined in Appendix B. Now the magnetic field is obtained by rearranging the modified BSCs and by multiplying by constant −/ as The total incident electric field Eq. ( 13) and total incident magnetic field from the source distribution, accounting for polarization as in Eq. ( 15), is written as the integration of vector spherical harmonics expansions Eq. ( 23) and Eq. ( 25) over the source surface as Similar scattered field presentations mapped with the T-matrix method are presented as where E  and H  are introduced in Appendix B.
In summary, we have expanded planar 2D ASM to 3D ASM by dividing the surface into infinitesimally small sub-regions used as local source points.We mapped source points' global location and orientation information into the modified BSCs and obtained the total fields as a superposition of VSH-expanded source points.Also, surface polarization can be modified by parameterization.

Results
The results are divided into three sections; first, we present the radiation patterns and nonsingularity of the differential sources.In the second section, we demonstrate the practicality and accuracy of the presented theory by comparing the Gaussian beam s morphology-dependent resonances from a sphere with the traditional 2D ASM and the presented 3D ASM, with a high agreement.In the last section, we synthesize an electromagnetic beam from an ellipsoidal surface, verify the incident field with physical optics simulations, and compute scattered fields from a 100-dielectric layer sphere.

Differential source
The total field from arbitrary surface Ω is obtained as a superposition of differential source fields E  = Ω E  .As shown in derivation, when the surface area of the differential source plane approaches zero, the result Eq. ( 14) is analagous to the angular spectrum presentation of a curl of the magnetic dipole multiplied with scalar Green's function E  = −∇ × (−2n × e 2 ).As is well known, the fields created in this way are exact from an infinitely large plane, and the method accurately approximates field synthesis from curved surfaces when the RoC is above a wavelength.
The radiation pattern of a differential source E  is presented, where the magnetic dipole is positioned to the origin along the y-axis.Radiated fields are computed at transverse xy-plane and yz-plane; see Figure (3).This arrangement's main polarization is along the x-axis, and the field propagates to the positive z-direction.Simulation shows the free placement of source points inside the computation domain due to the absence of singularities.The maximum amplitude of the electric fields is one without normalization.

Morphology-dependent resonance test
Mie scattering region is known for frequency-dependent backscatter intensity.This phenomenon is due to the scattering resonances, which trap energy temporally inside the sphere.Backscatter intensity is presented as a function of size parameter , where = ˛0  is a refractive index dependent wavenumber inside the sphere and  is the radius of the sphere.These backscatter resonances are called morphology-dependent resonances (MDRs) [30,31].MDRs are highly responsive to simulation errors and are a good test of the accuracy of the derived 3D ASM combined with VSH expansion [32].
High accuracy manifests with simulations, where a Gaussian beam is created from its beam waist and its field components on a plane outside the beam waist.In the latter case, each source point has varying propagation direction and complex amplitude, making field synthesis on this computation inefficient with traditional methods.MDRs are computed with both source plane locations when the beam waist is positioned at the origin of the homogeneous sphere and on the side of the coated sphere, defined precisely later.Likewise, MDRs from an individual source point are computed in both focus scenarios, see Figure (5).
Simulations are executed by illuminating homogeneous and coated spheres by Gaussian beam in size parameter range  ∈ [32,36] with beam waist radius  0 = 1.5.Homogeneous spheres have a refractive index of   = 1.36, and coated spheres have a   = 1.36 core with a   = 1.5 shell with radius of 0.7 .These values were selected from article [28] for comparing results.In both scenarios, the Gaussian beam is first synthesized on the transverse Plane 1 at the beam waist from the planar Gaussian electric field distribution with the traditional 2D AS and the proposed 3D AS method with an exact match.Then the incident field is computed on Plane 2 at −5  distance from the origin.The Cartesian electric and magnetic field components of the incident field at Plane 2 are presented in Figure (7).These components are synthesized from Plane 1. Gaussian beam MDRs synthesized from both source plane scenarios are computed at -500  distance from the origin [20] and compared in Figure (8).The resonances are labeled as   , or   , , where  indicates the mode number, and  the number of radial peaks in the angle-averaged internal energy density distribution [20].First, the homogeneous sphere is illuminated with a Gaussian beam focused at the origin.Second, a coated sphere is illuminated with a Gaussian beam focused on the side of the coated sphere  =  to verify even more specific resonance behavior.Also, the MDRs from a source point on Plane 2 in both scenarios are presented to underline the resonance difference with varying beam shapes.Electromagnetic backscatter intensities show excellent MDRs matching from both source plane synthesis scenarios.The resonance shape, location, and peak spacing in Figure 8 match with values from [28].In Figure (8 a), the peaks correspond to modes   39,1 ,   40,1 , and   41,1 respectively, and the spacing between the peaks satisfies the Δ = /2  ≈ 1.15 condition.Also, in Figure (8 b), the MDRs peaks from coated sphere correspond to the modes   39,1 ,   40,1 ,   41,1 and   42,1 respectively.
The MDRs of an individual differential source do not match the Gaussian beam resonances.However, a total beam created as a superposition of the differential sources from the surface Ω synthesizes the original beam with highly matching MDRs.This leads to the conclusion that the presented 3D AS method combined with the VSH expansion synthesizes electromagnetic scattering from multilayered spheres with high accuracy.

Elliptical source distribution
An elliptical surface was used as an example of a source distribution to synthesize a focused linearly polarized incident field using /6 sampling, see Figure (9).The incident field is compared to the Physical optics simulations at the yz-plane with an excellent agreement to less than a -39 dB difference in amplitude and less than the 0.1-degree difference in phase.Also, the scattered fields are computed from a  = 7.8 mm radius sphere in xz-plane with 100 dielectric layers with linearly varying permittivity   ∈ [1 − 0.001, 3 − 0.01] from the surface to the center.Base vectors for elliptical surfaces are derived from the Eq.(1 -4) with the parametrization with axial scaling factors   = 0.6,   = 1.5,   = 0.8 and an unscaled radius  = 3.The origin of the ellipsoid is shift ( 0 ,  0 ,  0 ) = (3, 0, −0.5) compared to the sphere.

Conclusions
We propose a method to synthesize electromagnetic beams by arbitrary surface electric field distributions.The goal is to compute scattering from multilayered spheres illuminated by incident fields that can be modified freely by adjusting the source distribution's shape, amplitude, phase, and location.The incident field is presented in VSH expansion defined by the BSCs, and the scattered field is obtained by mapping the incident field BSCs with the T-matrix method.
The goal is achieved by computing BSCs with the proposed 3D ASM derived from the nominal 2D ASM in the theory chapter.The key concept is to divide the arbitrary surfaces into differential elements, which are used as local AS sources.Location and orientation information of AS sources is transformed in a global coordinate system and mapped into the BSCs to preserve spherical symmetry for VSH expansion.As a result, the total incident field is obtained as a superposition of VSH-expanded differential sources, in which orientation, locations, and complex amplitude can be adjusted.
The proposed method approximates incident field synthesis from an arbitrary surface, considering the wavelength compared to the surface details.Then the scattered fields are rigorously computed from approximated incident fields with the Mie theory.In addition, two essential points of the method are (1) 3D ASM is analogous to the magnetic dipole presentation of the Stratton-Chu method, which satisfies Maxwell's equations and well approximates scattering/radiation from surfaces with RoC larger than /2, and (2) when the 3D ASM is applied to planar surfaces, it returns the original 2D AS representation.
Simulations verify the proposed method's incident field synthesis and scattering accuracy.First, we demonstrate the source points magnetic dipole behavior with the radiation patterns without singularities.Then we verify the scattering accuracy by comparing the MDRs resonances from spheres illuminated by the Gaussian beam created from the beam waist by nominal 2D ASM and by more complex electric field distribution with the proposed 3D ASM.The MDRs resonances of the Gaussian beam synthesized by the 3D ASM were in excellent agreement with the reference values, validating the high accuracy on the scattered fields.At last, we demonstrate the method's practicality by synthesizing an incident field from an elliptical surface and computing scattered fields from a 100-layer sphere.Additionally, we present mathematical proof to support the theory.
The novelty of this approach lies in a straightforward simulation algorithm, where the incident field can be defined at any parametrized surface.Furthermore, the parameterized surface can be located inside the computational domain, on the surface, or inside the spherical scatterer without restrictions and source singularities.Unrestricted placement of the source distribution has clear advantages, especially in the inverse beam design, where the desired beam is defined on the sub-region on the spherical scatterer.Additionally, the synthesized beam can be reradiated from any evaluation surface to another, enabling beam simulation between optical elements with the ability to consistently compute the scattered fields from multilayered spheres.

Appendix A
Theorem 5.1.Let the function  0 be continuous in a closed segment of surface Ω and let r ∈ R 3 be such that (r − o) • ( o   × o  ) ≠ 0 for all o ∈ Ω.The field  created by function  0 is presented by Eq.( 13) and Eq.( 14) at point r.
Proof.Let us define  > 0. Let us show first that when the partition Ω  is small enough, in other words, the areas |Ω  | are small enough regardless of t, by replacing the function   0 supported by each piece Ω  with a constant   0 ( 0, 0) =  0 ( , ), the error made in Eq. ( 8) is smaller than /4.Function  0 is uniformly continuous on the compact segment of surface Ω, and based on this, with sufficiently small |Ω  | it holds for all x, ȳ ∈ Ω  and , where the constant  will be defined later on.In the Fourier transform on the plane holds for all ȳ is integrable on ( x ,  ȳ ) plane, because outside of the circle  2 x +  2 ȳ =  2 , the exponential becomes real and negative.Let  be, at first, a continuous elementary function Due the monotone of integral,  (| z( , )|) ≤  on all  regardless of partitions.The upper limit for the error due  0 ( , ) on the Eq. ( 8) is obtained as where Eq. ( 29) and Eq. ( 32) have been used.Because of ||e 1 || = 1 and the triangle inequality the error in Eq. ( 9) is In the above approximation, a constant  0 ( , ) is used on piece Ω; as a function, it is zero on the plane outside of the piece.Let us denote this function simply as  0 ( , ).Next, we consider the Fourier transform of the function  0 ( , ).Because function ȳ is integrable, there can be found an origin centered closed disk B such that for all  because | z| ≥ z .For all ( x ,  ȳ ) ∈ R 2 it holds Let us now consider the error ( 2 ) made in integral Eq. ( 8) when F { 0 ( , )} is replaced by the constant  0 (   ,   )|Ω  |.Under the estimates Eq. ( 35) and Eq. ( 36) it follows For all ( x ,  ȳ ) ∈ B and ( x, ȳ) ∈ Ω  , we can estimate by continuity of the function Replacing by F { 0 ( , )} with the constant  0 ( , )|Ω  | in Eq. ( 8) we make an error where the estimates Eq. (37) and Eq. ( 39) have been used.The error for sum Eq. ( 9) is obtained as Thus, when the partition Ω  is dense enough and r ∈ R (42) Finally, by continuity of the integrand, the integral in Eq. ( 13) exists, and, because of the estimate Eq. ( 42), for the points r ∈ R 3 like in Theorem there also exists The proof above also holds for the second Eq. ( 15) component; only a few adaptions are needed.Moreover, Theorem 5.1 holds with an expectation r ∉ Ω for both Eq. ( 15) components, and the solutions are visible in the simulations.The required extended proof is reasonably complicated and will be presented in a separate article.Its difficulty lies primarily in that integrals do not exist in the usual sense when z = 0.

Appendix B
The VSH beam shape coefficients for each  − plane wave of the incident field is defined as [ where    is the associated Legendre function of the first kind of degree  and order .Additionally ,  are spherical coordinates and e  , e  spherical base vectors of r in the local base (e 1 , e 2 , e 3 ) as in Eq. (23) The scattered field from a differential element in VSH presentation is obtained as E  (r) =  2  ∑︁  ∑︁     emn M (3)  emn (r) +  omn M (3)  omn (r) +  emn N (3)  emn (r) +  omn N (3)  omn (r) , where the superscripts (3) present the vector spherical harmonics (outgoing wave) with spherical Hankel function of the first kind ℎ (1)  n ().The  emn ,  omn ,  emn and  omn are vector spherical harmonic coefficient for the scattered field calculated from the T-matrix method as where the T-matrix elements for the multilayered sphere are obtained by the algorithm defined in [34].

Fig. 1 .
Fig. 1.Vector relations between source point base-vectors on an arbitrary surface and a global origin.

Fig. 3 .
Fig. 3. Magnetic dipole in E  along the y-axis.The magnetic dipole is marked as a green arrow on the evaluation plane presented as a red area.On a) and b) are the simulation arrangements on xy-and yz-planes, respectively.

Figure ( 4 )
Figure (4) presents the radiated fields from a differential source on xy-and yz-planes.

Fig. 4 .
Fig. 4. Electromagnetic radiation from a source point, where a) electric and d) magnetic fields at the transverse xy-plane.On b) electric and e) magnetic fields on propagation yz-plane, and on c) electric and f) magnetic fields phases on propagation yz-plane.

Fig. 5 .
Fig. 5. Morphology-dependent resonance simulation of a sphere illuminated by Gaussian beam.On a) is a Gaussian beam synthesis on the source Plane 1 at the origin, and field components are computed on Plane 2. On b) is Simulation 2, where Plane 2 is used as a 3D ASM source, where each source point's local propagation direction is defined by the local Poynting vector (marked as blue), and source points are multiplied by complex amplitudes, obtained from Simulation 1.

Fig. 7 .
Fig. 7. Electric and magnetic field components of Gaussian beam one at the plane 2. On a)-c) are electric field components, and on d)-f) are magnetic field components.All the components are normalized by max |  |.

Fig. 8 .
Fig. 8. Morphology-dependent resonances of a Gaussian beam and a differential source illumination presented in Figure (4).On a) is a situation where the Gaussian beam waist is positioned at the origin of the homogeneous sphere with refractive index   = 1.36.On b) is a situation where the Gaussian beam waist is on the side of the coated sphere, with a core refractive index of   = 1.5 and shell as   = 1.36.The backscatter intensities of one source point are scaled to fit the same plot.

Figure ( 10 )
illustrate the amplitude and phase of the total field E  = E  + E  on xz-plane where ,  ∈ [−5, 5].Simulations are performed with 500 × 500 evaluation points with size parameter of  = 13 and VSH modes of  = 50.

Fig. 10 .
Fig. 10.Total fields from a 100-layer sphere illuminated with the incident beam defined in the ellipsoidal surface.a) total field magnitude, and b) total field phase.

Table 1 .
Characteristics of the selected methods in interest.
1,glob , e 2,glob and e 3,glob .Let the position  be presented in global coordinate system as r = e  + e  + e  .Position  presented in the local coordinate system with a global base is r − o = ( −   )e  + ( −   )e  + ( −   )e  .