Файл: Doicu A., Wriedt T., Eremin Y.A. Light scattering by systems of particles (OS 124, Springer, 2006.pdf

+G (ks, r, r ) · [E (r ) − E (r)] dV (r )

Dr

+G (ks, r, r ) − Gs (ks, r, r ) · E (r) dV (r )

Dr

The first term in (3.4) is a nonsingular integral that can be computed numerically. When the integration volume Di,m is spherical and su ciently small, we can set Di,m = Dr and the first term is zero. For the third term, we use the coordinate-free representations

where R = r − r , R = |R|, eR = R/R and g(R) = exp(jksR)/(4πR), and obtain

G (ks, r, r ) − Gs (ks, r, r ) · E (r) dV (r ) = M E (r)

Dr

with a = (1/ks2) (1/(4π|r −r|)), yields

194 3 Simulation Results

The usual tradition is to neglect the second term assuming that the field is constant or claiming that as the integration volume becomes small the integrand vanishes. As pointed out by Draine and Goodman [56], vanishing of the second term is of the same order as with the third term and cannot be omitted. Using the second-order Taylor expansion for each Cartesian component Ek of E,

Ek (r ) = Ek (r) + (r −r) · Ek (r) + 12 (r −r) · [(r −r) · ] Ek (r) ,

and the vector Helmholtz equation

E + m2r,mks2E = 0 , in Di,m

we obtain [185]

G (ks, r, r ) · [E (r ) − E (r)] dV (r ) = M1E (r)

Dr

with

Within each volume element, the total field can be approximated by

where W k are referred to as the shape or basis functions. Frequently, nodebased shape functions are used as basis functions, but edge-based shape functions, using an expansion of a vector field in terms of values associated with element faces, are better suited for simulating electromagnetic fields at corners and edges. A detailed mathematical analysis of edge-based shape functions for various twoand three-dimensional elements has been given by Graglia et al. [86]. Linear as well as higher order curvilinear elements have been discussed by Mur and deHoop [176], Lee et al. [139], Wang and Ida [249], Crowley et al. [42] and Sertel and Volakis [209]. In general, the edge-based shape functions can be categorized as either interpolatory or hierarchical. An element is referred to as hierarchical if the vector basis functions forming the nth order element are a subset of the vector basis functions forming the (n + 1)th order element. This property allows to use lowest order elements in regions where the field is expected to vary slowly and higher order elements in regions where rapid field variation is anticipated. For a tetrahedral element, a hierarchical edge-based element up to and including second order has been introduced by Webb and Forghani [259]. Additional to nodeand edge-based shape functions, low order vector spherical wave functions have been employed by Peltoniemi [185]. In our analysis we assume that the electromagnetic field is constant over each element and consider the simplest type of basis functions, the constant pulse functions: W k (r) = ek for r Di,m, and W k (r) = 0 for r / Di,m, where k = x, y, z and (ex, ey , ez ) are the Cartesian unit vectors. This approach is called collocation method, since a certain value of the electromagnetic field is assigned to the centroid rm of each element, i.e., E(r) ≈ E(rm) for r Di,m. Substituting this expansion into (3.5) and

−ks2 m2r,n − 1 (1 − δmn) VnG (ks, rm, rn) ,

and Vn is the volume of the n element.

In computing the self-term Em we extracted from the original integrand a singular function which approximates the integrand in the neighborhood of the singularity. However, there are other e cient schemes to accurately compute

196 3 Simulation Results

singular integrals [3]. For example, in the framework of the volume integral equation method, Sertel and Volakis [209] transferred one of the derivative of the scalar Green function to the testing function using the divergence theorem. The resulting integrals are evaluated by the annihilation method, which converts first-order singular integrals into nonsingular integrals by using an appropriate parametric transformation (change of variables).

Whereas the volume integral equation method deals with the total field, the discrete dipole approximation exploits the concept of the exciting field, each volume element being explicitly modeled as an electric dipole moment. In this regard, we consider the field that excites the volume element Di,m

It is clear from (3.7) that the exciting field consists of the incident field and the field contributions of all other volume elements Di,n with n = m. The heart of the discrete dipole approximation is (3.7), while (3.8) is used to express the total field in terms of the exciting field. In dyadic notations, (3.8) read as

Inserting (3.9) into (3.7), we obtain the expression of the field that excites the volume element Di,m in terms of the exciting fields of all other volume elements Di,n

× G (ks, rm, rn) · A−nn1 · Eexc (rn) .

Defining the polarizability dyad by

tn = εs m2r,n − 1 VnA−nn1

and taking into account that ks = k0√εs, yields

n=1

where the dyadic kernel is

Bmn = Iδmn − k02 (1 − δmn) G (ks, rm, rn) · tn .

In terms of the dipole moment pm

pm = εs m2r,m − 1 VmE (rm) ,

(3.9) can be written as

pm = tm · Eexc (rm)

and the e ect of the exciting field can be interpreted as inducing a dipole moment pm in each volume element Di,m. This is the origin of the term “discrete dipole approximation.”

Parenthetically we note that the discrete dipole approximation can be used to derive the size-dependent Maxwell–Garnett formula. For a collection of identical spheres of volume V and relative refractive index mr, the polarizability dyad is

The size-dependent Maxwell–Garnett formula can now be obtained by combining (3.11) and (2.168) [129].

The maximum linear cross-sectional extent 2rm of Di,m is such that

rmmr,m < 0.1 ,

λ

that is, the dimension rm is no more than a tenth of the wavelength in the volume element Di,m. This condition ensures that the spatial variations of the electromagnetic field inside each volume element are small enough, so that each volume element can be thought of as a dipole scatterer. Consequently, for large scale simulations, the systems of equations (3.6) and (3.10) are excessively large and iterative solvers as the conjugate or bi-conjugate gradient method are employed. The coupling submatrices in (3.6) and (3.10) depend on the absolute distance among elements, and the overall matrices are block T¨oplitz. Thus, (3.6) and (3.10) can be cast in circulant form, and the fast Fourier transform can be employed to carry out the matrix–vector products [68, 80].