1. Introduction - one dimensional discrete dipole approximation

The discrete-dipole approximation [3,4] is a general and exact method for computing scattering and absorption by a particle of arbitrary shape. In the last decade it allowed to routinely solve scattering problems by arbitrary non-spherical particles with size parameter (ratio of physical size to incomming wavelength) between x=-20. There is very strong interest in this method in variety of applications including atmospheric science, nanotechnology [9], chemistry, biology, electrical engineering, and astrophysics. The discrete dipole approximation (DDA) has surprisingly rich numerical and physical structure, mostly due to its simplicity, and to the fact that the boundary conditions are naturally included. In the limit of very many dipoles the method is exactly equivalent to integral form of Maxwell equations for light scattering problems. Some-what surprisingly, this simplicity of the DDA approach resulted in many exciting new extensions, for example, to periodic structures Chaumet [1] or particles placed on surfaces. There is also several of papers studying various aspects of the DDA numerics. There is no doubt that further improvements are possible. The discrete dipole approximation replaces the solid particle by an array of N point dipoles, with the spacing between dipoles small compared with the wavelength. Physical intuition suggest that the complexity of any solution algorithm can not be better than O(N). Current solvers indeed progressed to fast solvers with complexity proportional to O(N logN) but through iterative methods. Another outstanding problem is related to many incoming plane waves. Currently, for each of them one has to solve the iterative problem anew. Finally, for orientational averages over many incoming plane waves one would like to know full inverse of matrix describing the shape and composition of a particle. In this paper we attempt to show how to attack these outstanding problems in the one dimensional case of the discrete dipole approximation.

The most general DDA discretization [3] leads to set of 3N linear complex equations

where 3×3 matrix A_jk is defined as

and r _jk ≡|r _j-r _k| is the distance between respective position of dipoles, P _k has 3 components and it is polarization vector induced at dipole k, position of a dipole is defined by 3-component vector r _k and k is the wavenumber (certain care needs to be taken to differentiate index k and wavenumber k). For dipoles interacting with itself we define

where α_k is polarizability of dipole k. Equation (2) can be rewritten in the following form

where n _jk=r _jk/r_jk is the direction pointing from the position of dipole j towards dipole k. For special case of dipoles placed on x-axis vector n _jk=(1,0,0) for all dipoles, i.e.

Functions f (r_jk ) and g(r_jk ) are

and

and the sum of these two is

We will now consider simple example of just 3 dipoles placed on x-axis. According to Eq. (1) this leads to 9 equations

$\left(\begin{array}{ccccccccc}{\alpha }^{-1}& 0& 0& h_{12}& 0& 0& h_{13}& 0& 0\\ 0& {\alpha }^{-1}& 0& 0& f_{12}& 0& 0& f_{13}& 0\\ 0& 0& {\alpha }^{-1}& 0& 0& f_{12}& 0& 0& f_{13}\\ h_{21}& 0& 0& {\alpha }^{-1}& 0& 0& h_{23}& 0& 0\\ 0& f_{21}& 0& 0& {\alpha }^{-1}& 0& 0& f_{23}& 0\\ 0& 0& f_{21}& 0& 0& {\alpha }^{-1}& 0& 0& f_{23}\\ h_{31}& 0& 0& h_{32}& 0& 0& {\alpha }^{-1}& 0& 0\\ 0& f_{31}& 0& 0& f_{32}& 0& 0& {\alpha }^{-1}& 0\\ 0& 0& f_{31}& 0& 0& f_{32}& 0& 0& {\alpha }^{-1}\end{array}\right)\left(\begin{array}{c}{p}_{x,1}\\ p_{y,1}\\ p_{z,1}\\ p_{x,2}\\ p_{y,2}\\ p_{z,2}\\ p_{x,3}\\ p_{y,3}\\ p_{z,3}\end{array}\right)$

where we used (4) and definition (7). This can be re-arranged by combining x,y,z components separately to read

It is easy to see from Eq. (8) that x,y, z-components are decoupled and can be solved separately. Interactions defined by Eqs. (5) and (6) do not depend on position of dipoles but only on their distance h_jk ≡h(|r_j -r_k |). This is crucial simplification as we can write the interaction forces in a very simple way. For example, the x-component of the N-element polarization vector is given by

where entries on diagonals are determined by the first row and column only. Such matrix is called Toeplitz and leads to very efficient solution method which we will discuss shortly. Not

only is this matrix Toeplitz but it is also symmetric, i.e. h_jk =h_kj . This is because the dipole-dipole interactions depend only on the absolute distance between them. Thus, (9) is fully determined by its first row. I.e.

where notation of the left hand side indicates that Toeplitz matrix A is generated by column {α ^-1, h ₁₂, h ₁₃, …, h ₁ N}. Such structure holds for the y and z components. Equation (10) is rather simple case of DDA which nevertheless exhibits its interesting numerical aspects. To the best of our knowledge all current fast algorithms to solve (9) employ iterative technique of guessing polarization vector and gradually improve on it. The algorithm requires numerous matrix times vector multiplications which can be efficiently done using Fast Fourier Transform because the matrix has Toeplitz structure. However, one wonders if the inverse problem has simple structure. This is indeed the case in one dimensions as the following discussion shows.

2. Displacement structure of Toeplitz matrices and its relationship to Fast Fourier Transform.

In this section we summarize recent results related to numerical algorithms to invert and multiply Toeplitz matrix. The notation and most results follow work of Gohberg and Olshevsky [6,7]. The algorithms will be presented in terms of matrix multiplications and so we need to define several matrices of special structure. The most general Toeplitz

matrix is defined by its first column a _c=[a ₁, a ₂, …,a _n-1] and its first row a _r=[a ₀, a _-1, …,a _-n+1] . In the case of symmetric matrix the first column is the same as the first row a _c=a _r. The φ-cyclic circulant has form

where vector r=(r_i ${)}_{i=0}^{n-1}$. Circulant (φ=1) is completely defined by its first column r and each subsequent column is right shifted. We are particularly interested in φ-circulants because they are closely related to the Fast Fourier Transform (FFT). Multiplication of φ-circulant matrix by vector is fast and requires three FFTs as the following result [2] shows

where

here diag(r) is diagonal matrix, ξⁿ =φ, and ℱ is the Fourier matrix and the asterix ^⋆means conjugate transpose. For strictly circulant matrix φ=1 and all elements of D_φ are 1, thus

and multiplication of Circ(r) by a vector b is done by taking FFT of b, i.e. ℱ(b), multiplying this vector term by term with ℱ(r), and finally taking inverse FFT of the resulting vector. Equation (13) follows the same pattern of three FFTs and multiplications except that each element of vectors b and r is first scaled by elements of vector D_φ =[1, ξ, ξ ², …, ξ ^n-1]. Even though (13) is written in terms of matrix times matrix multiplications in practice this is never done. Rather, it is a series of matrix times vector multiplications which are performed. Because it is fast to multiply by cyclic circulants our plan is to reduce multiplication and inversion of a Toeplitz matrix in terms of circulants.

We begin with simple derivations. Let us write matrix (11) in terms of circulant (φ=1)and skew-circulant (φ=-1). For this we introduce two vectors which are almost reverse of the vectors (a_c,a_r ) which generate a Toeplitz matrix. These are defined as b_c =[0, a _-n+1, … a _-2, a _-1] and b_r =[0, a _n-1, … a ₂, a ₁]. By adding and subtracting A(b) we can rewrite A(a) as

Notice that matrix A(a)+A(b)=A(a+b) and A(a)-A(b)=A(a-b). By direct inspection we see that A(a+b)=Circ₁(a+b) and A(a-b)=Circ_-1(a-b). Therefore A(a) multiplied by a vector d is given by

where we introduced definition of convolution of two vectors in terms of multiplication of circulant by a vector Conv(p,q)≡Circ(p)q. Equation (16) and (13) shows that multiplication of a Toeplitz matrix by a vector can be done with 6 FFTs and scaling (for skew-circulant part). In iterative techniques, it is required to multiply A(a)d for many d′s the FFTs of vectors (a+b) and (a-b) can be precalculated and the complexity of multiplication by a vector reduces to 4 FFTs of length n. Equations (16) forms a basis for almost all current implementations of the discrete dipole approximation. These implementations depend on iterative calculations of A(a)d with vector d changed according to conjugate gradient scheme. There are variants of (16) which rely on embedding matrix A(a) to larger 2n dimensional matrix and zero padding but they are equivalent to the scheme given above. Now, we will report an important result related to the inverse of Toeplitz matrix. To this end we need to introduce φ-cyclic lower shift matrix

$Z_{\phi }=\left(\begin{array}{ccccc}0& \dots & \dots & \dots & \phi \\ 1& 0& & & ⋮\\ 0& 1& \ddots & & \\ & \ddots & & \ddots & ⋮\\ 0& \dots & 0& 1& 0\end{array}\right)$

and assume that we can calculate two vectors x and y defined by

where e ₀=(1, 0, …, 0, 0) and e _n-1=(0, 0, …, 0, 1) and x ₀≠0 then inverse of matrix A is given by

This result is modification of celebrated Gohberg-Semencul [6,7] formula expressed in terms of circulants. Equation (18) provides very fast algorithm for inverse matrix times vector multiplication A ^-1 d in terms of 4 convolutions. It shows that one needs just two auxiliary vectors to form the Toeplitz matrix inverse.

3. Proposed algorithms

We will now propose two algorithms to solve discrete dipole approximation in a prototype, one-dimensional case.

3.1. Algorithm 1. Conjugate gradient (iterative)

Consider Eq. (10) A(h)P=E which is set of equations for unknown polarization vector P and known right hand side E and Toeplitz matrix defined by N-element vector h=(α ^-1 h ₁₂ h ₁₃ … h _1N) associated with function (7). Equations (16) and (13) show that multiplication of a Toeplitz matrix by a vector can be done with 6 FFTs and scaling (for skew-circulant part). Using conjugate gradient iterations we have to multiply A(h)P for many incomming plane parallel waves E. In that case the FFTs of vectors (h+b) and (h-b), where b is defined in Equation (15), can be precalculated and the complexity of each multiplication by a vector reduces to 4 FFTs of length N. This algorithm is similar to that implemented in the DDSCAT code [4]. However, the algorithm we propose here requires only N-dimensional FFT rather than 2N-dimensional version used in the DDSCAT code [8]. There is also economy in memory use because one doesn’t have to form extended matrix of dimension 2N.

3.2. Algorithm 2. Two-phase inversion (the main result).

Consider Eq. (10) A(h)P=E which is a set of equations for unknown polarization vector P and known right hand side E and Toeplitz matrix defined by N-element vector h=(α ^-1, h ₁₂, h ₁₃, …, h _1N) defined by function (7). There are two phases of the algorithm.

Phase 1. In the first phase one has to calculate and store two additional vectors defined in (17). This can be done iteratively by using conjugate gradient method (Algorithm 1).

Phase 2. The second phase applies A ^-1(h)E_inc using expression (18). This requires 4 convolutions and each of them can be performed by using 3 FFTs (see Eq. 13) with scaling. Each subsequent case of the incomming plane wave can be solved by employing Phase 2. For many right hand sides of (10) this algorithm is faster in comparison to conjugate gradient Algorithm 1.

4. Discussion

In this paper we propose new algorithms for solution of light scattering on non-spherical particles using one-dimensional variant of the discrete dipole approximation. We discuss recent advances in algorithms for matrices with structures and show that it is possible to apply these advances to use FFT and non-iterative solvers. Algorithm 1 can be extended to 3D for arbitrary non-spherical particles and should be competitive with current implementations. There are variants [10] of Algorithm 2 which do not require iterative solution for auxiliary vectors x and y defined by (17). Future challenge is an extension of Algorithm 2 to the three-dimensional case, i.e. block-Toeplitz matrices [5]. Extension to 3D parallelepiped for homogeneous and isotropic refractive index should be straightforward. Apparently, the way non-spherical particles are extended to parallelepiped requires introduction of vacuum dipoles with refractive index equal to 1 which ruins the Toeplitz structure on main diagonal.