Discrete-dipole approximation for periodic targets: theory and tests

Discrete-dipole approximation for periodic targets: theory and tests Bruce T. Draine 1,* and Piotr J. Flatau 2 1Princeton University Observatory, Princeton, New Jersey 08544-1001, USA 2Scripps Institution of Oceanography, University of California San Diego, La Jolla, California 92093, USA *Corresponding author: draine@astro.princeton.edu Received June 24, 2008; accepted August 19, 2008; posted September 3, 2008 (Doc. ID 97834); published October 14, 2008 The discrete-dipole approximation (DDA) is a powerful method for calculating absorption and scattering by targets that have sizes smaller than or comparable to the wavelength of the incident radiation. The DDA can be extended to targets that are singly or doubly periodic. We generalize the scattering amplitude matrix and the 4 /H110034 Mueller matrix to describe scattering by singly and doubly periodic targets and show how these ma- trices can be calculated using the DDA. The accuracy of DDA calculations using the open-source code DDSCAT is demonstrated by comparison with exact results for infinite cylinders and infinite slabs. A method for using the DDA solution to obtain fields within and near the target is presented, with results shown for infinite slabs. © 2008 Optical Society of America OCIS codes: 050.1755, 050.5298, 260.0260, 290.5825 . 1. INTRODUCTION Electromagnetic scattering is used to study isolated par- ticles, but increasingly to characterize extended targets ranging from nanostructure arrays in laboratories to planetary and asteroidal regoliths. To model the absorp- tion and scattering, Maxwell’s equations must be solved for the target geometry . For scattering by isolated particles with complex geom- etry , a number of different theoretical approaches have been used, including the discrete-dipole approximation (DDA) [1–4], also known as the coupled-dipole approxima- tion or coupled-dipole method. The DDA can treat inho- mogeneous targets and anisotropic materials, and has been extended to treat targets near substrates [ 5,6]. Other techniques have also been employed, including the finite-difference time-domain (FDTD) method [ 7,8]. For illumination by monochromatic plane waves, the DDA can be extended to targets that are spatially periodic and (formally) infinite in extent. This could apply , for ex- ample, to a periodic array of nanostructures in a labora- tory setting, or it might be used to approximate a regolith by a periodic array of “target unit cells” with complex structure within each unit cell. Generalization of the DDA (or coupled-dipole method) to periodic structures was first presented by Markel [ 9] for a 1-dimensional chain of dipoles, and more generally by Chaumet et al. [10], who calculated the electric field near a 2-dimensional array of parallelepipeds illuminated by a plane wave. Chaumet and Sentenac [ 11] further ex- tended the DDA to treat periodic structures with a finite number of defects. From a computational standpoint, solving Maxwell’s equations for periodic structures is only slightly more dif- ficult than calculating the scattering properties of a single “unit cell” from the structure. Since it is now feasible to treat targets with N /H1140710 6 dipoles (target volume /H11407200/H92613, where /H9261is the wavelength), it becomes possible to treat extended objects with complex substructure. The objective of the present paper is to present the theory of the DDA applied to scattering and absorption by structures that are periodic in one or two spatial dimen- sions. We also generalize the standard formalism for de- scribing the far-field scattering properties of finite targets (the 2 /H110032 scattering amplitude matrix and 4 /H110034 Mueller matrix) to describe scattering by periodic targets. We show how to calculate the Mueller matrix to describe scat- tering of arbitrarily polarized radiation. The theoretical treatment developed here has been implemented in the open-source code DDSCAT 7 (see Appendix A). The theory of the DDA for periodic targets is reviewed in Section 2, and the formalism for describing the far-field scattering properties of periodic targets is presented in Sections 3–5. Transmission and reflection coefficients for targets that are periodic in two dimensions are obtained in Section 6. The applicability and accuracy of the DDA method are discussed in Sections 7 and 8, where we show scattering properties calculated using DDSCAT 7 for two geometries for which exact solutions are available for comparison: (1) an infinite cylinder and (2) an infinite slab of finite thick- ness. The numerical comparisons demonstrate that, for given /H9261, the DDA converges to the exact solution as the interdipole spacing d → 0. 2. DDA FOR PERIODIC TARGETS The DDA is a general technique for calculating scattering and absorption of electromagnetic radiation by targets with arbitrary geometry . The basic theory of the DDA has been presented elsewhere [ 3]. Conceptually , the DDA con- sists of approximating the target of interest by an array of polarizable points with specified polarizabilities. Once the B. T. Draine and P . J. Flatau Vol. 25, No. 11 / November 2008 / J. Opt. Soc. Am. A 2693 1084-7529/08/112693-11/$15.00 © 2008 Optical Society of America

polarizabilities are specified, Maxwell’s equations can be solved accurately for the dipole array . When applied to fi- nite targets, the DDA is limited by the number of dipoles N for which computations are feasible—the limitations may arise from large memory requirements or the large amount of computing that may be required to find a solu- tion when N is large. In practice, modern desktop comput- ers are capable of solving the DDA equations, as imple- mented in DDSCAT [3,12], for N as large as /H11011106. Developed originally to study scattering from isolated, finite structures such as dust grains [ 1], the DDA can be extended to treat singly or doubly periodic structures. Consider a collection of N polarizable points defining a target unit cell (TUC). Now consider a target consisting of a 1-dimensional or 2-dimensional periodic array of identi- cal TUCs, as illustrated in Figs. 1 and 2; we will refer to these as 1-D or 2-D targets, although the constituent TUC may have an arbitrary 3-dimensional shape. For a mono- chromatic incident plane wave E inc/H20849r,t/H20850= E0 exp/H20849ik0 · r − i/H9275t/H20850, /H208491/H20850 the polarizations of the dipoles in the target will oscillate coherently . Each dipole will be affected by the incident wave plus the electric field generated by all of the other point dipoles. Let index j =1,..., N run over the dipoles in a single TUC, and let indices m, n run over replicas of the TUC. The /H20849m , n/H20850replica of dipole j is located at r jmn = rj00 + mLu + nLv , /H208492/H20850 where Lu and Lv are the lattice vectors for the array . For 1-D targets we let m vary , but set n = 0. For 2-D targets, the area per TUC is ATUC /H11013/H20841Lu /H11003Lv/H20841= LuLv sin /H9258uv , /H208493/H20850 where /H9258uv is the angle between Lu and Lv. The replica dipole polarization Pjmn/H20849t/H20850is phase-shifted relative to Pj00/H20849t/H20850: Pjmn/H20849t/H20850= Pj00/H20849t/H20850exp/H20851i/H20849mk0 · Lu + nk0 · Lv/H20850/H20852. /H208494/H20850 Define a matrix A such that − Aj,kmnPkmn gives the elec- tric field E at rj00 produced by an oscillating dipole Pkmn located at rkmn. Expressions for the 3 /H110033 tensor elements of A have been presented elsewhere (e.g., [ 3]); A depends on the target geometry and wavelength of the incident ra- diation, but not on the target composition nor on the di- rection or polarization state of the incident wave. Using Eq. (2) we may construct a matrix A ˜ such that, for j /HS11005k,− A˜ j,kPk00 gives the electric field at rj00 produced by a dipole Pk00 and all of its replica dipoles Pkmn, and for j = k it gives the electric field at rj00 produced only by the replica dipoles: A˜ j,k = /H20858 m=−/H11009 /H11009 /H20858 n=−nmax nmax /H208491− /H9254jk/H9254m0/H9254n0/H20850Aj,kmn exp/H20851i/H20849mk0 · Lu + nk0 · Lv/H20850/H20852, /H208495/H20850 where nmax = 0 for 1-D targets and nmax = /H11009for 2-D targets, and /H9254ij is the Kronecker delta. For /H20841m/H20841, /H20841n/H20841→ /H11009, location j00 is in the radiation zone of dipole kmn, and the electric field falls off in magnitude only as 1 / r. The sums in Eq. (5) would be divergent were it not for the oscillating phases of the terms, which ensure convergence. Evalua- tion of these sums can be computationally demanding when k 0Ly or k0Lz are small. Chaumet et al. [10] have dis- cussed methods for efficient evaluation of these sums. We evaluate Eq. (5) numerically by introducing a factor exp/H20851−/H20849/H9253k0r/H208504/H20852to smoothly suppress the contributions from large r, and truncating the sums: Fig. 1. (a) Target consisting of a 1-D array of TUCs and (b) showing how an infinite cylinder can be constructed from disk- like TUCs (lower). The M = 0 scattering cone with /H9251s = /H92510 is illustrated. Fig. 2. (a) Target consisting of a 2-D array of TUCs and (b) showing how an infinite slab is created from TUCs consisting of a single “line” of dipoles. 2694 J. Opt. Soc. Am. A / Vol. 25, No. 11 / November 2008 B. T. Draine and P . J. Flatau

A˜ j,k /H11015/H20858 m,n /H11032Aj,kmn exp/H20851i/H20849mk0 · Lu + nk0 · Lv/H20850− /H20849/H9253k0rj,kmn/H208504/H20852, /H208496/H20850 where rj,kmn /H11013/H20841rkmn − rj00/H20841and the summation is over /H20849m , n/H20850with rj,kmn /H333552/ /H9253k0, i.e., out to distances where the suppression factor exp /H20851−/H20849/H9253k0r/H208504/H20852/H11015e−16. For given k0, Lu, Lv, the A˜ j,k depend only on rj00 − rk00, and therefore only O/H208498N/H20850distinct A˜ j,k require evaluation. Ideally , one would use a very small value for the inter- action cutoff parameter /H9253, but the number of terms [ /H11008/H9253−1 for 1-D or /H11008/H9253−2 for 2-D] in Eq. (6) diverges as /H9253→ 0. We show that setting /H9253/H110150.001 ensures accurate results for the cases studied here. The polarizations Pj00 of the dipoles in the TUC must satisfy the system of equations Pj00 = /H9251j/H20875Einc/H20849rj/H20850− /H20858 k/HS11005j A˜ j,kPk00/H20876. /H208497/H20850 If there are N dipoles in one TUC, then Eq. (7) is a system of 3 N linear equations where the polarizability tensors /H9251j are obtained from lattice-dispersion-relation theory [13,14]. After A˜ has been calculated, Eq. (7) can be solved for Pj00 using iterative techniques when N /H112711. 3. IN THE RADIATION ZONE In the radiation zone kr /H112711, the electric field due to dipole jmn is Ejmn = k0 2 exp/H20849ik0/H20841r − rjmn/H20841/H20850 /H20841r − rjmn/H20841/H208751− /H20849r − rjmn/H20850/H20849r − rjmn/H20850 /H20841r − rjmn/H208412 /H20876Pjmn , /H208498/H20850 /H20841r − rjmn/H20841= /H20851r2 −2 r · rjmn + rjmn 2 /H208521/2 /H11015r/H208771− r · rjmn r2 + 1 2r2 /H20875rjmn 2 − /H20873 r · rjmn r /H20874 2 /H20876+ ¯ /H20878. /H208499/H20850 Define the unit vector kˆ s /H11013ks / k0. We seek to sum the con- tribution of all the dipoles to the electric field propagating in direction kˆ s. At location r = rkˆ s, the dominant contribu- tion will be from dipoles located within the Fresnel zone (see, e.g., [ 15]), which will have a transverse radius R F /H11015/H20849r / k0/H208501/2. For dipoles within the Fresnel zone, 1 /H20841r − rjmn/H20841/H208751− /H20849r − rjmn/H20850/H20849r − rjmn/H20850 /H20841r − rjmn/H208412 /H20876Pjmn /H11015 1 r /H208511− kˆ skˆ s/H20852Pjmn . /H2084910/H20850 Thus, in the radiation zone, E/H20849r/H20850= k0 3 k0r exp/H20849ik0r/H20850/H208511− kˆ skˆ s/H20852/H20858 j Pj00/H20858 m,n exp/H20849i/H9274jmn/H20850, /H2084911/H20850 with /H9274jmn /H11013mk0 · Lu + nk0 · Lv − ks · rjmn + k0 2r /H20851rjmn 2 − /H20849kˆ s · rjmn/H208502/H20852/H20849 12/H20850 /H11015− ks · rj00 + m/H20849k0 − ks/H20850· Lu + n/H20849k0 − ks/H20850· Lv + 1 2k0r /H20851m2/H20849k0 2 − ksu 0 /H20850Lu 2 + n2/H20849k0 2 − ksv 2 /H20850Lv 2 +2 mn/H20849k0 2Lu · Lv − ksuksvLuLv/H20850/H20852 + O/H20873 mL r /H20874, /H2084913/H20850 where ksu /H11013ks · Lu / Lu, ksv /H11013ks · Lv / Lv, and terms of order /H20849mL / r/H20850may be neglected because mL / r /H11011RF / r /H11008r−1/2 as r → /H11009. Thus, for r → /H11009, the electric field produced by the os- cillating dipoles is Es/H20877 k0 2 r exp/H20849ik0r/H20850/H208511− kˆ skˆ s/H20852/H20858 j Pj00 exp/H20849− iks · rj00/H20850/H20878G/H20849r,ks/H20850, /H2084914/H20850 G/H20849r,ks/H20850/H11013/H20858 m,n exp/H20849i/H9021mn/H20850, /H2084915/H20850 /H9021mn /H11013m/H20849k0 − ks/H20850· Lu + n/H20849k0 − ks/H20850· Lv + 1 2k0r /H20851m2/H20849k0 2 − ksu 2 /H20850Lu 2 + n2/H20849k0 2 − ksv 2 /H20850Lv 2 +2 mn/H20849k0 2Lu · Lv − ksuksvLuLv/H20850/H20852. /H2084916/H20850 It is convenient to define FTUC/H20849kˆ s/H20850/H11013k0 3/H208511− kˆ skˆ s/H20852/H20858 j=1 N Pj00 exp/H20849i/H9275t − iks · rj00/H20850, /H2084917/H20850 so that the electric field produced by the dipoles is Es = exp/H20849iks · r − i/H9275t/H20850 k0r FTUC/H20849kˆ s/H20850G/H20849r,ks/H20850. /H2084918/H20850 FTUC depends on the scattering direction kˆ s and also on the direction of incidence kˆ 0 and polarization E0 of the in- cident wave. A. Isolated Finite Target: /H9263=0 We will refer to finite isolated targets—consisting only of the dipoles in a single TUC—as targets that are periodic in /H9263= 0 dimensions. For this case, we simply set G =1 i n Eq. (14); the scattered electric field in the radiation zone is Es = exp/H20849ik0r − i/H9275t/H20850 k0r FTUC/H20849kˆ s/H20850. /H2084919/H20850 The time-averaged scattered intensity is B. T. Draine and P . J. Flatau Vol. 25, No. 11 / November 2008 / J. Opt. Soc. Am. A 2695

Is = c/H20841Es/H208412 8/H9266= c 8/H9266k0 2r2 /H20841FTUC/H208412 , /H2084920/H20850 and the differential scattering cross section is dCsca d/H9024= 1 k0 2 /H20841FTUC/H208412 /H20841E0/H208412 . /H2084921/H20850 B. Target Periodic in One Dimension: /H9263=1 Without loss of generality , we may assume that targets with 1-D periodicity repeat in the yˆ direction. It is easy to see from Eqs. (15) and (16) that G = 0 except for scattering directions satisfying ksy = k0y + M 2/H9266 Ly M =0 ,±1 ,±2 ,. . ., /H2084922/H20850 where energy conservation /H20849ks 2 = k0 2/H20850limits the allowed val- ues of the integer M: /H20849− k0y − k0/H20850 Ly 2/H9266/H33355M /H33355/H20849− k0y + k0/H20850 Ly 2/H9266. /H2084923/H20850 If /H20849k0 + /H20841k0y/H20841/H20850Ly /H110212/H9266, then only M = 0 scattering is allowed. Define polar angles /H92510 and /H9251s for the incident and scat- tered radiation, so that k0y = k0 cos /H92510 , /H2084924/H20850 ksy = k0 cos /H9251s . /H2084925/H20850 For each allowed value of ksy, the scattering directions de- fine a cone: ks = ksyyˆ + /H20849k0 2 − ksy 2 /H208501/2 sin /H9251s sin /H92510 /H20851/H20849kˆ 0 − yˆ cos /H92510/H20850cos /H9256 + yˆ /H11003kˆ 0 sin /H9256/H20852, /H2084926/H20850 where /H9256is an azimuthal angle measured around the tar- get axis yˆ . The sum for G/H20849r , ks/H20850is [since M must be an in- teger; see Eq. (22)] G = /H20858 m=−/H11009 /H11009 exp/H20849i/H9021m0/H20850 = /H20858 m=−/H11009 /H11009 exp/H20875−2 /H9266iMm + i m2 2k0r /H20849k0 2 − ksy 2 /H20850Ly 2 /H20876/H2084927/H20850 → lim /H9280→0+ /H20885 −/H11009 /H11009 dm exp/H20875 i/H208491+ i/H9280/H20850 2k0r m2/H20849k0 2 − ksy 2 /H20850Ly 2 /H20876/H2084928/H20850 = /H208492/H9266ik0r/H208501/2 /H20849k0 2 − ksy 2 /H208501/2Ly = /H208492/H9266ik0r/H208501/2 k0Ly sin /H9251s , /H2084929/H20850 and the scattered electric field Es = /H20873 2/H9266i k0r /H20874 1/2 exp/H20849ik0r − i/H9275t/H20850 k0Ly sin /H9251s FTUC/H20849kˆ s/H20850/H2084930/H20850 shows the expected r−1/2 behavior far from the scatterer (the distance from the cylinder axis is R = r sin /H9251s). For each allowed value of kˆ s, the total time-averaged scattered power P¯ sca, per unit length of the target, per unit azimuthal angle /H9256, may be written d2P¯ sca dLd/H9256= /H20841E0/H208412 8/H9266c d2Csca dLd/H9256, /H2084931/H20850 where the differential scattering cross section is d2Csca dLd/H9256= 8/H9266 /H20841E0/H208412 1 c d2P¯ sca dLd/H9256= 8/H9266 /H20841E0/H208412c /H20841E/H208412c 8/H9266R sin /H9251s /H2084932/H20850 = 2/H9266 k0 3Ly 2 /H20841FTUC/H208412 /H20841E0/H208412 . /H2084933/H20850 C. Target Periodic in Two Dimensions: /H9263=2 For targets that are periodic in two dimensions, it is ap- parent from Eqs. (15) and (16) that G = 0 unless /H20849ks − k0/H20850· Lu =2 /H9266M, M =0 ,±1 ,±2 ,. . . , /H2084934/H20850 /H20849ks − k0/H20850· Lv =2 /H9266N, N =0 ,±1 ,±2 ,. . . . /H2084935/H20850 The 2-D target constitutes a diffraction grating, with scat- tering allowed only in directions given by Eqs. (34) and (35). It is convenient to define the reciprocal lattice vec- tors u /H11013 2/H9266xˆ /H11003Lv xˆ · /H20849Lu /H11003Lv/H20850 , v /H11013 2/H9266xˆ /H11003Lu xˆ · /H20849Lv /H11003Lu/H20850 . /H2084936/H20850 The wave vector transverse to the surface normal is ks/H11036/H11013k0/H11036+ Mu + Nv. /H2084937/H20850 Energy conservation requires that ksx 2 = k0 2 − /H20841k0/H11036+ Mu + Nv/H208412 /H110220. /H2084938/H20850 For any /H20849M , N/H20850allowed by Eq. (38), there are two allowed values of ksx differing by a sign; one (with ksxk0x /H110220) cor- responds to the /H20849M , N/H20850component of the transmitted wave, and the other (with ksxk0x /H110210) to the /H20849M , N/H20850compo- nent of the reflected wave. Define sin /H92510 /H11013 /H20841k0x/H20841 k0 , /H2084939/H20850 sin /H9251s /H11013 /H20841ksx/H20841 k0 . /H2084940/H20850 Note that /H92510 = /H9266/ 2 for normal incidence, and /H92510 → 0 for grazing incidence. For ks satisfying Eqs. (34) and (35),w e have 2696 J. Opt. Soc. Am. A / Vol. 25, No. 11 / November 2008 B. T. Draine and P . J. Flatau

G = /H20858 m,n exp/H20849i/H9021mn/H20850= /H20858 m,n exp/H20877 i 2k0r /H20851/H20849k0 2 − ksu 2 /H20850Lu 2m2 + /H20849k0 2 − ksv 2 /H20850Lv 2n2 +2 /H20849k0 2Lu:Lv − ksuksvLuLv/H20850mn/H20852/H20878 → lim /H9280→0+ /H20885 −/H11009 /H11009 dm/H20885 −/H11009 /H11009 dn exp/H20877 i/H208491+ i/H9280/H20850 2k0r /H20851/H20849k0 2 − ksu 2 /H20850Lu 2m2 + /H20849k0 2 − ksv 2 /H20850Lv 2n2 +2 /H20849k0 2Lu · Lv − ksuksvLuLv/H20850mn/H20852/H20878 = lim /H9280→0+ 1 ATUC /H20885 −/H11009 /H11009 dy/H20885 −/H11009 /H11009 dz exp/H20877 i/H208491+ i/H9280/H20850 2k0r /H20851k0 2/H20849y2 + z2/H20850− /H20849ksyy + kszz/H208502/H20852/H20878= 2/H9266ir k0ATUC sin /H9251s . /H2084941/H20850 The scattered electric field is Es = 2/H9266i exp/H20849iks · r − i/H9275t/H20850 k0 2ATUC sin /H9251s FTUC/H20849kˆ s/H20850. /H2084942/H20850 Note that /H20841Es/H20841is independent of distance x = r sin /H9251s from the target, as expected for a target that is infinite in two directions. For pure forward scattering, k s = k0, we must sum the incident wave Einc and the radiated wave Es: E = exp/H20849ik0 · r − i/H9275t/H20850/H20875E0 + 2/H9266iFTUC/H20849kˆ s = kˆ 0/H20850 k0 2ATUC sin /H9251s /H20876. /H2084943/H20850 The cross section per unit target area A for scattering into direction /H20849M , N/H20850is (for ks /HS11005k0) dCsca/H20849M,N/H20850 dA = /H20841E/H208412 sin /H9251s /H20841E0/H208412 sin /H92510 /H2084944/H20850 = 4/H92662 k0 4ATUC 2 sin /H92510 sin /H9251s /H20841FTUC/H20849kˆ s/H20850/H208412 /H20841E0/H208412 , /H2084945/H20850 where Csca can be evaluated for either transmitted or re- flected waves. For the special case ks = k0, the transmis- sion coefficient T/H20849M , N/H20850is obtained from the total forward-propagating wave (43): T/H208490,0 /H20850= 1 /H20841E0/H208412 /H20879E0 + 2/H9266iFTUC/H20849kˆ s = kˆ 0/H20850 k0 2ATUC sin /H92510 /H20879 2 . /H2084946/H20850 4. SCATTERING AMPLITUDE MATRICES Si „/H9263d… A. Isolated Finite Targets: /H9263=0 In the radiation zone, the scattered electric field is related to the incident electric field via the 2 /H110032 scattering ampli- tude matrix [ 16], defined so that /H20875 Es · eˆ s/H20648 Es · eˆ s/H11036 /H20876= i exp/H20849iks · r − i/H9275t/H20850 k0r /H20875 S2 /H208490d/H20850S3 /H208490d/H20850 S4 /H208490d/H20850S1 /H208490d/H20850/H20876/H20875 E0 · eˆ i/H20648 E0 · eˆ i/H11036 /H20876, /H2084947/H20850 where eˆ i/H11036= eˆ s/H11036/H11013 kˆ s /H11003kˆ 0 /H20841kˆ s /H11003kˆ 0/H20841 = kˆ s /H11003kˆ 0 1− /H20849kˆ s · kˆ 0/H208502 =− /H9278ˆ s , /H2084948/H20850 eˆ i/H20648/H11013kˆ 0 /H11003eˆ i/H11036= kˆ s − /H20849kˆ s · kˆ 0/H20850kˆ 0 1− /H20849kˆ s · kˆ 0/H208502 , /H2084949/H20850 eˆ s/H20648/H11013kˆ s /H11003eˆ s/H11036= − kˆ 0 + /H20849kˆ s · kˆ 0/H20850kˆ s 1− /H20849kˆ s · kˆ 0/H208502 = /H9258ˆ s /H2084950/H20850 are the usual conventions for the incident and scattered polarization vectors parallel and perpendicular to the scattering plane (see, e.g., Sec. 3.2 of [ 16]). B. Target Periodic in One Dimension: /H9263=1 For targets with 1-D periodicity , it is natural to generalize the scattering amplitude matrix, so that—for directions kˆ s for which scattering is allowed—the scattered electric field at a distance R = r sin /H9251s from the target is /H20875 Es · eˆ s/H20648 Es · eˆ s/H11036 /H20876= i exp/H20849iks · r − i/H9275t/H20850 /H20849k0R/H208501/2 /H20875 S2 /H208491d/H20850S3 /H208491d/H20850 S4 /H208491d/H20850S1 /H208491d/H20850/H20876/H20875 E0 · eˆ i/H20648 E0 · eˆ i/H11036 /H20876 /H2084951/H20850 for ks satisfying Eqs. (22)–(26). C. Target Periodic in Two Dimensions: /H9263=2 For targets with 2-D periodicity , it is natural to generalize the scattering amplitude matrix so that, for directions ks /HS11005kˆ 0 for which scattering is allowed, we write /H20875 Es · eˆ s/H20648 Es · eˆ s/H11036 /H20876= i exp/H20849iks · r − i/H9275t/H20850/H20875 S2 /H208492d/H20850S3 /H208492d/H20850 S4 /H208492d/H20850S1 /H208492d/H20850/H20876/H20875 E0 · eˆ i/H20648 E0 · eˆ i/H11036 /H20876 /H2084952/H20850 for ks satisfying Eqs. (37) and (38). For the special case of forward scattering ( M = N = 0 and ks = k0), where the scattering plane is not simply defined by k0 and ks, it is natural to use k0 and the target normal xˆ to define the scattering plane. Thus eˆ i/H11036= eˆ s/H11036/H11013 k0 /H11003ks/H11036 /H20841k0 /H11003ks/H11036/H20841, /H2084953/H20850 with eˆ i/H20648and eˆ s/H20648defined by Eqs. (49) and (50). For r → /H11009 /H20875 E · eˆ s/H20648 E · eˆ s/H11036 /H20876= i exp/H20849ik0 · r − i/H9275t/H20850/H20875 /H20849S2 /H208492d/H20850− i/H20850 0 0 /H20849S1 /H208492d/H20850− i/H20850/H20876 /H11003/H20875 E0 · eˆ i/H20648 E0 · eˆ i/H11036 /H20876. /H2084954/H20850 D. Far-Field Scattering Amplitude Matrices The scattering amplitude matrices Si /H20849/H9263d/H20850/H20849ks/H20850are directly related to the FTUC for the three cases /H9263=0,1,2: S1 /H20849/H9263d/H20850= C/H9263eˆ s/H11036· FTUC/H20849kˆ s,E0 = eˆ i/H11036/H20850, /H2084955/H20850 B. T. Draine and P . J. Flatau Vol. 25, No. 11 / November 2008 / J. Opt. Soc. Am. A 2697

S2 /H20849/H9263d/H20850= C/H9263eˆ s/H20648· FTUC/H20849kˆ s,E0 = eˆ i/H20648/H20850, /H2084956/H20850 S3 /H20849/H9263d/H20850= C/H9263eˆ s/H20648· FTUC/H20849kˆ s,E0 = eˆ i/H11036/H20850, /H2084957/H20850 S4 /H20849/H9263d/H20850= C/H9263eˆ s/H11036· FTUC/H20849kˆ s,E0 = eˆ i/H20648/H20850, /H2084958/H20850 C0 =− i, /H2084959/H20850 C1 =− /H20873 2/H9266i sin /H9251s /H20874 1/2 i k0Ly , /H2084960/H20850 C2 = 2/H9266 k0 2ATUC sin /H9251s , /H2084961/H20850 for finite targets /H20849C0/H20850and targets that are periodic in one or two dimensions ( C1 or C2). 5. FAR-FIELD SCATTERING MATRIX FOR STOKES VECTORS For a given scattering direction kˆ s, the 2 /H110032 complex am- plitude matrix Si /H20849/H9263d/H20850/H20849ks/H20850fully characterizes the far-field scattering properties of the target. The far-field scattering properties of an isolated finite target are characterized by the 4 /H110034 dimensionless Mueller matrix S /H9251/H9252 /H208490d/H20850, with the Stokes vector of radiation scattered into direction kˆ s at a distance r from the target given by Isca,/H9251/H11013 1 /H20849k0r/H208502 /H20858 /H9252=1 4 S/H9251/H9252 /H208490d/H20850Iinc,/H9252, /H2084962/H20850 where Iinc,/H9252= /H20849I , Q , U , V/H20850inc is the Stokes vector for the ra- diation incident on the target. For 1-D targets, we define the dimensionless scattering matrix S/H9251/H9252 /H208491d/H20850by Isca,/H9251/H11013 1 k0R /H20858 /H9252=1 4 S/H9251/H9252 /H208491d/H20850Iinc,/H9252, /H2084963/H20850 where R is the distance from the 1-D target. For 2-D tar- gets, we define S/H9251/H9252 /H208492d/H20850by Isca,/H9251/H11013/H20858 /H9252=1 4 S/H9251/H9252 /H208492d/H20850Iinc,/H9252. /H2084964/H20850 The 4 /H110034 scattering intensity matrix S/H9251/H9252 /H20849/H9263d/H20850is obtained from the scattering amplitude matrix elements Si /H20849/H9263d/H20850.E x - cept for the special case of forward scattering /H20849ks = k0/H20850for 2-D targets, the equations are the same as Eq. (3.16) of Bohren and Huffman [ 16]. For example, S 11 /H20849/H9263d/H20850= 1 2 /H20849/H20841S1 /H20849/H9263d/H20850/H208412 + /H20841S2 /H20849/H9263d/H20850/H208412 + /H20841S3 /H20849/H9263d/H20850/H208412 + /H20841S4 /H20849/H9263d/H20850/H208412/H20850, /H2084965/H20850 S21 /H20849/H9263d/H20850= 1 2 /H20849/H20841S2 /H20849/H9263d/H20850/H208412 − /H20841S1 /H20849/H9263d/H20850/H208412 − /H20841S4 /H20849/H9263d/H20850/H208412 + /H20841S3 /H20849/H9263d/H20850/H208412/H20850, /H2084966/H20850 S14 /H20849/H9263d/H20850=I m/H20849S2 /H20849/H9263d/H20850S3 /H20849/H9263d/H20850* − S1 /H20849/H9263d/H20850S4 /H20849/H9263d/H20850*/H20850. /H2084967/H20850 For the special case of forward scattering /H20849kˆ s = kˆ 0/H20850for 2-D targets, it is necessary to replace S1 /H208492d/H20850and S2 /H208492d/H20850with /H20849S1 /H208492d/H20850− i/H20850and /H20849S2 /H208492d/H20850− i/H20850[cf. Eq. (54)]. Thus, for example, S11 /H208492d/H20850/H20849ks = k0/H20850= 1 2 /H20849/H20841S1 /H208492d/H20850− i/H208412 + /H20841S2 /H208492d/H20850− i/H208412 + /H20841S3 /H208492d/H20850/H208412 + /H20841S4 /H208492d/H20850/H208412/H20850. /H2084968/H20850 6. TRANSMISSION AND REFLECTION COEFFICIENTS FOR 2-D TARGETS For targets with 2-D periodicity , it is natural to define generalized transmission and reflection coefficients for the Stokes vectors: for scattering order /H20849M , N/H20850, I sca,/H9251 = /H20858/H9252T/H9251/H9252/H20849M , N/H20850Iinc,/H9252is the Stokes vector component /H9251for radiation with ksxkinc,x /H110220, and R/H9251/H9252/H20849M , N/H20850is the fraction of the incident Stokes vector component /H9252that emerges in Stokes vector component /H9251with ksxkinc,x /H110210. These can be related to the S/H9251/H9252 /H208492d/H20850: R/H9251/H9252/H20849M,N/H20850= sin /H9251s sin /H92510 S/H9251/H9252 /H208492d/H20850 for ksxk0x /H110210, /H2084969/H20850 T/H9251/H9252/H20849M,N/H20850= sin /H9251s sin /H92510 S/H9251/H9252 /H208492d/H20850 for ksxk0x /H110220. /H2084970/H20850 The fraction of the incident power that is absorbed by the target is Pabs/Area /H20841E0/H208412c sin /H92510/8/H9266=1− /H20858 M,N /H20858 /H9252=1 4 /H20851R1/H9252/H20849M,N/H20850+ T1/H9252/H20849M,N/H20850/H20852 Iinc,/H9252 Iinc,1 , /H2084971/H20850 where Iinc,/H9252is the Stokes vector of the incident radiation. For unpolarized incident radiation, R11/H20849M , N/H20850is the fraction of the incident power that is reflected in diffrac- tion component /H20849M , N/H20850, T 11/H20849M , N/H20850is the fraction that is transmitted in component /H20849M , N/H20850, and 1 − /H20858M,N/H20851R11/H20849M , N/H20850+ T11/H20849M , N/H20850/H20852is the fraction of the inci- dent power that is absorbed. 7. EXAMPLE: INFINITE CYLINDER DDSCAT 7 has been used to calculate scattering and ab- sorption by an infinite cylinder consisting of a periodic ar- ray of disks of thickness d and period L y = d (where d is the interdipole spacing). Figure 3(a) shows S11 /H208491d/H20850for re- fractive index m = 1.33+ 0.01i, /H9266D / /H9261=5 0 (D is the cylinder diameter and /H9261the wavelength of the incident radiation), and incidence angle /H92510 = 60°. Because k0/H208491+ /H20841cos /H92510/H20841/H20850d /H110212/H9266, Eqs. (22) and (23) allow only M = 0 scattering, with /H9251s = /H92510. Also shown is the exact solution, calculated using a code written by Mackowski [ 17]. Light scattering by cyl- inders is generally described by scattering amplitudes Ti; in Appendix B we provide expressions relating these Ti to the Si used here. Figure 3(b) shows the fractional error in S11 /H208491d/H20850calculated using DDSCAT.A s d is decreased, the er- rors decrease. Excellent accuracy is obtained when the validity criterion [ 3] /H20841m/H20841kd /H113510.5 is satisfied: the fractional error in S 11 is typically less than a few%, except near deep minima in S11. Figure 3(c) shows S21 /H208491d/H20850, characterizing scattering of un- polarized light into the Stokes parameter Q (S21 /H110210 cor- responds to linear polarization perpendicular to the scat- tering plane). DDSCAT 7 and the exact solution are in very 2698 J. Opt. Soc. Am. A / Vol. 25, No. 11 / November 2008 B. T. Draine and P . J. Flatau

good agreement when /H20841m/H20841kd /H113510.5. Note that although the error in S21 /H208491d/H20850/H20849/H9258=0 /H20850/H110152 is large compared to S21/H208490/H20850=−6, this is small compared to S11 /H208491d/H20850/H208490/H20850/H110151500: the scattered ra- diation is only slightly polarized. The results in Fig. 3 were obtained using /H9253= 0.001 to truncate the integrations. To see how the results depend on /H9253, Fig. 4 shows S11 /H208491d/H20850computed for the problem of Fig. 3 but using different values of /H9253. For azimuthal angles /H9256 /H1102220°, the results for /H9253= 0.005 and 0.001 are nearly indis- tinguishable; the difference between the computed result and the exact solution is evidently due to the finite num- ber of dipoles used, rather than the choice of cutoff pa- rameter /H9253. However, the results for forward scattering are more sensitive to the choice of /H9253, as is seen in Figs. 4(c) and 4(d): it is necessary to reduce /H9253to 0.001 to attain high accuracy in the forward-scattering directions. Table 1 gives the CPU times to calculate A˜ , to then it- eratively solve the scattering problem to a fractional error /H1102110 −5 (using double-precision arithmetic), and finally to evaluate the scattering intensities for several of the cases shown in Figs. 3 and 4. For most cases the CPU time is dominated by the iterative solution using the conjugate gradient algorithm. While the time required to evaluate A ˜ might be reduced using the strategies suggested by [ 10], this step is generally a subdominant part of the computa- tion for targets with kd /H114070.1. Fig. 4. (Color online) Scattering by an infinite cylinder with diameter D and m = 1.33+ 0.01i, for radiation with /H9266D / /H9261= 50 and incidence angle /H92510 = 60°. (a) Exact solution (solid curve) and DDA results for D / d = 512 and various values of the interaction cutoff parameter /H9253.( b ) Fractional error in S11 /H208491d/H20850. (c), (d) Same as (a), (b), but expanding the region 0 /H11021/H9256/H1102120°. For this case, results computed with /H9253= 0.002 and 0.001 are nearly indistinguishable. Fig. 3. (Color online) Scattering by an infinite cylinder with diameter D and m = 1.33+ 0.01i, for radiation with x = /H9266D / /H9261= 50 and inci- dence angle /H92510 = 60°. (a) S11 /H208491d/H20850: solid curve, exact solution; broken curves, DDA results for D / d = 256, 360, and 512 ( N = 51,676, 102,036, 206,300 dipoles per TUC). (b) Fractional error in S11 /H208491d/H20850/H20849DDA/H20850. (c) S21 /H208491d/H20850. (d) Error in S21 /H208491d/H20850. B. T. Draine and P . J. Flatau Vol. 25, No. 11 / November 2008 / J. Opt. Soc. Am. A 2699

The above results have been for a weakly absorbing cyl- inder. To confirm that the DDA can be applied to strongly absorbing material, Fig. 5 shows scattering calculated for a cylinder with m =2+ i and x = /H9266D / /H9261= 25. Once again, the accuracy is very good, with small fractional errors pro- vided /H20841m/H20841kd /H113510.5. 8. EXAMPLE: PLANE-PARALLEL SLAB Consider a homogeneous plane-parallel slab with thick- ness h and refractive index m. Radiation incident on it at angle of incidence /H9258i will either be specularly reflected or transmitted. The reflection and transmission coefficients R and T can be calculated analytically , taking into ac- count multiple reflections within the slab [ 15]. With an exact solution in hand, we can evaluate the accuracy of the DDA applied to this problem. Figure 6 shows results for two cases: a dielectric slab with m = 1.50, and an ab- sorbing slab with m = 1.50+ 0.02i. DDSCAT 7 was used to calculate reflection, transmission, and absorption by an infinite slab, generated from a TUC consisting of a single line of dipoles extending in the x di- rection, with Ly = Lz = d. The selection rules (34), (35), and (38) allow only M = N = 0: transmission or specular reflec- tion. The reflection and transmission coefficients for ra- diation polarized parallel or perpendicular to the plane containing kˆ and the surface normal are R /H20648= S11 /H208491d/H20850/H20849ksx =− k0x/H20850+ S12 /H208491d/H20850/H20849ksx = k0x/H20850, /H2084972/H20850 R/H11036= S11 /H208491d/H20850/H20849ksx =− k0x/H20850− S12 /H208491d/H20850/H20849ksx = k0x/H20850, /H2084973/H20850 T/H20648= S11 /H208491d/H20850/H20849ksx = k0x/H20850+ S12 /H208491d/H20850/H20849ksx = k0x/H20850, /H2084974/H20850 T/H11036= S11 /H208491d/H20850/H20849ksx = k0x/H20850− S12 /H208491d/H20850/H20849ksx = k0x/H20850. /H2084975/H20850 The DDA results are in excellent agreement with the ex- act results when the validity condition /H20841m/H20841k0d /H110210.5 is sat- isfied, but results with moderate accuracy are obtained even when /H20841m/H20841k 0d /H110151. 9. NEAR-FIELD EVALUATION The polarizations Pj00 can be used to calculate the electric and magnetic fields at any point, including within or near the target, using the exact expression for E and B from a point dipole, modified by a function /H9278: E/H20849r,t/H20850= e−i/H9275t /H20858 j /H20858 m,n /H11032exp/H20849ik0Rjmn/H20850 /H20841Rjmn/H208413 /H9278/H20849Rjmn/H20850/H20877k0 2Rjmn /H11003/H20849Pjmn /H11003Rjmn/H20850+ /H208491− ik0Rjmn/H20850 Rjmn 2 /H208513Rjmn/H20849Rjmn · Pjmn/H20850 − Rjmn 2 Pjmn/H20852/H20878+ E0 exp/H20849ik0 · r − i/H9275t/H20850, /H2084976/H20850 B/H20849r,t/H20850= e−i/H9275t /H20858 j /H20858 m,n /H11032k2 exp/H20849ik0Rjmn/H20850 Rjmn 2 /H9278/H20849Rjmn/H20850/H20849Rjmn /H11003Pjmn/H20850 /H11003/H208731− 1 ik0Rjmn /H20874+ kˆ 0 /H11003E0 exp/H20849ik0 · r − i/H9275t/H20850, /H2084977/H20850 Rjmn /H11013r − rjmn , /H2084978/H20850 /H9278/H20849R/H20850/H11013exp/H20851− /H9253/H20849k0R/H208504/H20852/H11003/H20877 1 for R /H33356d /H20849R/d/H208504 for R /H11021d/H20878. /H2084979/H20850 The function /H9278/H20849R/H20850smoothly suppresses the (oscillating) contribution from distant dipoles in order to allow the summations to be truncated, just as in Eq. (6) for evalua- tion of A ˜ j,k.I f r is within the target or near the target sur- face, the summations over /H20849m , n/H20850are limited to /H20841Rjmn/H20841 /H333552/ /H9253k0. The /H20849R / d/H208504 factor suppresses the R−3 divergence of E as r approaches the locations of individual dipoles, and at the dipole locations results in E that is exactly equal to the field that is polarizing the dipoles in the DDA formulation. Evaluation of Eqs. (76) and (77) is computa- tionally intensive, because the summations /H20858 j/H20858m,n/H11032typi- cally have many terms. Table 1. CPU Time to Calculate Scattering by m = 1.33+ 0.01i Infinite Cylinders on Single-Core 2.4 GHz AMD Opteron Model 250 /H9266D / /H9261N /H9253 Calc. A˜ (min) Solution (min) Scat. (min) Total (min) 25 51,676 0.005 3.29 17.6 0.65 22.2 25 51,676 0.001 16.2 17.8 0.65 35.3 50 102,036 0.005 4.58 59.7 1.27 66.8 50 102,036 0.001 22.8 45.2 1.09 70.2 50 206,300 0.005 13.2 254 2.88 273 50 206,300 0.001 66.0 292 2.76 364 Fig. 5. (Color online) Light scattered by an infinite cylinder with m =2+ i for radiation with x =2 /H9266R / /H9261= 25 and incidence angle /H92510 = 60°. (a) S11 /H208491d/H20850: solid curve, exact solution; broken curves, DDA results for D / d = 128, 180, and 256 ( N = 12,972, 25,600, 51,676 dipoles per TUC). (b) Fractional error in S11 /H208491d/H20850. 2700 J. Opt. Soc. Am. A / Vol. 25, No. 11 / November 2008 B. T. Draine and P . J. Flatau

To illustrate the accuracy , we consider the infinite slab of Fig. 6(b), with refractive index m = 1.5+ 0.02i and radia- tion incident at an angle /H9258i = 40°. Figure 7 shows the time- averaged /H20841E/H208412 / /H20841E0/H208412 for slab thickness h = 0.2/H9261—near a minimum in transmission, and a maximum in reflection [see Fig. 6(b)]. The program DDFIELD (see Appendix A) was used to evaluate E along two lines normal to the slab: track 1 passes directly through dipole sites, and track 2 passes midway between the four nearest dipoles as it crosses each dipole layer. The E fields calculated along tracks 1 and 2 are very similar, although of course not identical. Within the slab, /H20841E/H20841along track 2 tends to be slightly smaller than along track 1, but for this example the difference is typically less than /H110111%. Figure 7 shows results for the slab represented by h / d = N x = 10 and 20 di- pole layers (with /H20841m/H20841kd = 0.19 and 0.094, respectively). Even for Nx = 10, the electric field at points more than a distance d from the edge is obtained to within /H110112% accu- racy at worst, which is perhaps not surprising because, as seen in Fig. 6, the calculated transmission and reflection coefficients are very accurate. The discontinuity in /H20841E/H20841 2 at the boundary is spread out over a distance /H11011d. The DDA Fig. 6. (Color online) Transmission and reflection coefficients for radiation of wavelength /H9261incident at angle /H9258i = /H20849/H9266/2− /H92510/H20850= 40° relative to the normal on a slab with thickness h, incident E /H20648and /H11036to the scattering plane, as a function of Re /H20849m/H20850h / /H9261. (a) Nonabsorbing slab with m = 1.5. (b) Absorbing slab with m = 1.5+ 0.02i. Solid curves exact solution; symbols, results calculated with the DDA using dipole spacing d = h / 10, h / 20, and h / 40. Fig. 7. (Color online) /H20841E2/H20841/ /H20841E0 2/H20841within and near the dielectric slab of Fig. 6(b) for slab thickness h =0 . 2/H9261, incidence angle /H9251i = 40°, and incident polarizations /H20648and /H11036to the scattering plane. Results were calculated using Eq. (76) with the slab represented by Nx = 10 and Nx = 20 dipole layers (i.e., dipole spacing d =0 . 1h and 0.05 h). The circles along track 1 are at points where dipoles are located. B. T. Draine and P . J. Flatau Vol. 25, No. 11 / November 2008 / J. Opt. Soc. Am. A 2701

obviously cannot reproduce field structure near the target surface on scales smaller than the dipole separation d, but fields on scales larger than d appear to be quite accu- rate. DDSCAT and DDFIELD should be useful tools for study- ing electromagnetic fields around arrays of nanostruc- tures, such as gold nanodisks [ 18,19]. 10. SUMMARY The principal results of this study are as follows: 1. The DDA is generalized to treat targets that are pe- riodic in one or two spatial dimensions. Scattering and ab- sorption of monochromatic plane waves can be calculated using algorithms that parallel those used for finite tar- gets. 2. A general formalism is presented for description of far-field scattering by targets that are periodic in one or two dimensions using scattering amplitude matrices and Mueller matrices that are similar in form to those for fi- nite targets. 3. The accuracy of the DDA for periodic targets is tested for two examples: infinite cylinders and infinite slabs. The DDA, as implemented in DDSCAT 7 , is accurate provided the validity criterion /H20841m/H20841kd /H113510.5 is satisfied. 4. We show how the DDA solution can be used to evalu- ate E and B within and near the target, with calculations for an infinite slab used to illustrate the accuracy of near- field calculations. APPENDIX A: DDSCAT 7 AND DDFIELD The theoretical developments reported here have been implemented in a new version of the open-source code DDSCAT [20]. DDSCAT 7 is written in Fortran 90 with dy- namic memory allocation and the option to use either single- or double-precision arithmetic. DDSCAT 7 includes options for various target geometries, including a number of periodic structures. A program DDFIELD for near-field calculations is also provided. DDSCAT 7 offers the option of using an implementation of BICGSTAB with enhancement to maintain convergence in finite precision arithmetic [ 21]. The matrix-vector mul- tiplications A˜ P are accomplished efficiently using FFTs [22]. Documentation for DDSCAT is available from ArXiv [23], with additional information available from [ 24]. In addition to differential scattering cross sections, DDSCAT reports dimensionless “efficiency factors” Qx /H11013Cx/H20849TUC/H20850/ /H9266aeff 2 for scattering and absorption, where Cx/H20849TUC/H20850is the total cross section for scattering or absorp- tion per TUC, normalized by /H9266aeff 2 , where aeff /H11013/H208493VTUC /4 /H9266/H208501/3 is the radius of a sphere with volume equal to the solid volume VTUC in one TUC. In the case of 1-D targets with periodicity Ly in the y direction, the absorption, scattering, and extinction cross sections per unit target length are dC x dL = 1 Ly Qx/H9266aeff 2 /H20849A1/H20850 for x = ext, sca, and abs, where Qx are the efficiency factors calculated by DDSCAT. In the case of 2-D targets with periodicities Lu and Lv, the absorption, scattering, and extinction cross sections per unit target area are dC x dA = Qx/H9266aeff 2 LuLv sin /H9258uv . /H20849A2/H20850 APPENDIX B: RELATION BETWEEN Si AND Ti FOR INFINITE CYLINDERS The analytic solution for infinite cylinders decomposes the incident and scattered radiation into components po- larized parallel and perpendicular to planes containing the cylinder axis and the propagation vector kˆ 0 or kˆ s. These polarization basis states differ from the choice that is usual for scattering by finite particles, where it is cus- tomary to decompose the incident and scattered waves into components polarized parallel and perpendicular to the scattering plane—the plane containing kˆ 0 and kˆ s. In the notation of Bohren and Huffman [ 16], the radia- tion scattered by an infinite cylinder can be written /H20875 Es · eˆ s/H20648 /H20849ck/H20850 Es · eˆ s/H11036 /H20849ck/H20850/H20876= i exp/H20849iks · r − i/H9275t/H20850/H20873 2i /H9266k0R sin /H9251/H20874 1/2 /H20875 T1 − T3 T3 T2 /H20876 /H11003/H20875 E0 · eˆ i/H20648 /H20849ck/H20850 E0 · eˆ i/H11036 /H20849ck/H20850/H20876, /H20849B1/H20850 where R is the distance from the cylinder axis, /H9251is the angle between k0 and the cylinder axis c, and superscript /H20849ck/H20850denotes polarization vectors parallel or perpendicular to planes containing the cylinder axis cˆ and either k0 or ks. The azimuthal angle /H9256is measured around the cylin- der axis cˆ , with /H9256= 0 for forward scattering. The scatter- ing angle /H9258= arccos/H20851kˆ 0 · kˆ s/H20852is /H9258= arccos/H208511− /H208491 − cos /H9256/H20850sin2 /H9251/H20852. /H20849B2/H20850 The scattering amplitude matrix elements Ti appearing in Eq. (B1) can be related to the matrix elements Si ap- pearing in Eq. (51): S = /H20873 2i /H9266sin /H9251/H20874 1/2 ATB−1 , /H20849B3/H20850 S /H11013/H20875 S2 /H208491d/H20850S3 /H208491d/H20850 S4 /H208491d/H20850S1 /H208491d/H20850/H20876, T /H11013/H20875 T1 − T3 T3 T2 /H20876, /H20849B4/H20850 A /H11013/H20875 eˆ s/H20648· eˆ s/H20648 /H20849ck/H20850eˆ s/H20648· eˆ s/H11036 /H20849ck/H20850 eˆ s/H11036· eˆ s/H20648 /H20849ck/H20850eˆ s/H11036· eˆ s/H11036 /H20849ck/H20850/H20876 = 1 sin /H9258/H20875 − cot /H9251/H208491 − cos /H9258/H20850sin /H9251sin /H9256 − sin /H9251sin /H9256− cot /H9251/H208491 − cos /H9258/H20850/H20876, /H20849B5/H20850 B /H11013/H20875 eˆ i/H20648· eˆ i/H20648 /H20849ck/H20850eˆ i/H20648· eˆ i/H11036 /H20849ck/H20850 eˆ i/H11036· eˆ i/H20648 /H20849ck/H20850eˆ i/H11036· eˆ i/H11036 /H20849ck/H20850/H20876 = 1 sin /H9258/H20875 cot /H9251/H208491 − cos /H9258/H20850sin /H9251sin /H9256 − sin /H9251sin /H9256cot /H9251/H208491 − cos /H9258/H20850/H20876, /H20849B6/H20850 2702 J. Opt. Soc. Am. A / Vol. 25, No. 11 / November 2008 B. T. Draine and P . J. Flatau

B−1 = sin /H9258 cot2 /H9251/H208491 − cos /H9258/H208502 + sin2 /H9251sin2 /H9256 /H11003/H20875 cot /H9251/H208491 − cos /H9258/H20850− sin /H9251sin /H9256 sin /H9251sin /H9256cot /H9251/H208491 − cos /H9258/H20850/H20876. /H20849B7/H20850 ACKNOWLEDGMENTS This research was supported in part by National Science Foundation (NSF) grant AST-0406883, and by the Office of Naval Research (ONR). We thank Dan Mackowski for providing his code for light scattering by infinite cylin- ders, H. A. Y ousif for discussions concerning scattering by infinite cylinders, and the referees for helpful comments. REFERENCES 1. E. M. Purcell and C. R. Pennypacker, “Scattering and absorption of light by nonspherical dielectric grains,” Astrophys. J. 186, 705–714 (1973). 2. B. T. Draine, “The discrete-dipole approximation and its application to interstellar graphite grains,” Astrophys. J. 333, 848–872 (1988). 3. B. T. Draine and P . Flatau, “Discrete-dipole approximation for scattering calculations,” J. Opt. Soc. Am. A 11, 1491–1499 (1994). 4. B. T. Draine, “The discrete dipole approximation for light scattering by irregular targets,” in Light Scattering by Nonspherical Particles: Theory, Measurements, and Applications, M. I. Mishchenko, J. W . Hovenier, and L. D. Travis, eds. (Academic, 2000), pp. 131–145. 5. R. Schmehl, B. M. Nebeker, and E. D. Hirleman, “Discrete- dipole approximation for scattering by features on surfaces by means of a two-dimensional fast Fourier transform technique,” J. Opt. Soc. Am. A 14, 3026–3036 (1997). 6. M. Paulus and O. J. F . Martin, “Green’s tensor technique for scattering in two-dimensional stratified media,” Phys. Rev . E 63, 066615 (2001). 7. P . Y ang and K. N. Liou, “Finite difference time domain method for light scattering by nonspherical and inhomogeneous particles,” in Light Scattering by Nonspherical Particles: Theory, Measurements, and Applications, M. I. Mishchenko, J. W . Hovenier, and L. D. Travis, eds. (Academic, 2000), pp. 173–221. 8. A. Taflove and S. C. Hagness, Advances in Computational Electrodynamics: the Finite-Difference Time-Domain Method (Artech House, 2005). 9. V . A. Markel, “Coupled-dipole approach to scattering of light from a one-dimensional periodic dipole structure,” J. Mod. Opt. 40, 2281–2291 (1993). 10. P . C. Chaumet, A. Rahmani, and G. W . Bryant, “Generalization of the coupled dipole method to periodic structures,” Phys. Rev . B 67, 165404 (2003). 11. P . C. Chaumet and A. Sentenac, “Numerical simulations of the electromagnetic field scattered by defects in a double- periodic structure,” Phys. Rev . B 72, 205437 (2005). 12. B. T. Draine and P . Flatau, “User Guide for the Discrete Dipole Approximation Code DDSCAT6.1,” http://arXiv .org/abs/ astro-ph/0409262 (2004). 13. B. T. Draine and J. Goodman, “Beyond Clausius-Mossotti— Wave propagation on a polarizable point lattice and the discrete dipole approximation,” Astrophys. J. 405, 685–697 (1993). 14. D. Gutkowicz-Krusin and B. T. Draine, “Propagation of electromagnetic waves on a rectangular lattice of polarizable points,” http://arXiv .org/abs/astro-ph/0403082 (2004). 15. M. Born and E. Wolf, Principles of Optics (Cambridge U. Press, 1999). 16. C. F . Bohren and D. R. Huffman, Absorption and Scattering of Light by Small Particles (Wiley , 1983). 17. D. Mackowski, private communication (2007). 18. Z. N. Utegulov, J. M. Shaw, B. T. Draine, S. A. Kim, and W . L. Johnson, “Surface-plasmon enhancement of Brillouin light scattering from gold-nanodisk arrays on glass,” Proc. SPIE 6641, 66411M (2007). 19. W . L. Johnson, S. A. Kim, Z. N. Utegulov, and B. T. Draine, “Surface-plasmon fields in two-dimensional arrays of gold nanodisks,” Proc. SPIE 7032 (2008) (to be published). 20. http://www.astro.princeton.edu/~draine/DDSCAT.html (2008). 21. M. A. Botchev, SUBROUTINE ZBCG2 , http://www.math.uu.nl/ people/vorst/zbcg2.f90 (2001). 22. J. J. Goodman, B. T. Draine, and P . J. Flatau, “Application of fast-Fourier transform techniques to the discrete dipole approximation,” Opt. Lett. 16, 1198–1200 (1990). 23. B. T. Draine and P . Flatau, “User Guide for the Discrete Dipole Approximation Code DDSCAT 7.0 ,” http://arXiv .org/abs/ 0809.0337 (2008). 24. http://ddscat.wikidot.com. B. T. Draine and P . J. Flatau Vol. 25, No. 11 / November 2008 / J. Opt. Soc. Am. A 2703