Light scattering by an ensemble of interacting dipolar particles with both electric and magnetic polarizabilities

We have studied the problem of light scattering by an ensemble of dipoles with both electric and magnetic polarizabilities. Using the coupled electric and magnetic dipole method as the formal base, we have generalized the eigenvector decomposition of the local dipole moments previously derived for purely electric particles to the case of both electric and magnetic dipoles. We have analyzed the properties of eigenvalues and eigenvectors in the most elementary case of two particles. In the purely electric case, the eigenvalues correspond to the resonance modes of the system due to the electromagnetic coupling of its components. For a two-dipole system with both electric and magnetic responses, purely electric, purely magnetic, and mixed states can be distinguished. The resonance spectrum is analyzed as a function of the magnetic permeability, and it is shown that the latter can be fitted quite accurately by the eigenmode decomposition.


I. INTRODUCTION
The study and development of metamaterials or, more specifically, of the left-handed materials ͑LHMs͒ which are capable of negative refraction ͓1͔ is currently a research field of enormous interest from both the basic and applied points of view.Many efforts have been made to bring the spectral range in which metamaterials exhibit the required properties into the visible part of the electromagnetic spectrum ͓2-5͔.
The key for obtaining such a negative refractive index is to achieve a magnetic permeability 1.For naturally occurring materials, = 1 in the optical range, but experiments as described in Ref. ͓4͔ show that effective values of the magnetic permeability different from 1 can be obtained in the visible for specially designed systems.Those systems are built from basic elements ͑split ring resonators ͓6,7͔, parallel nanorods ͓8͔, nano fishnets ͓9͔, etc.͒ made from classic materials having = 1, but due to their specific design, a magnetic response is obtained, and on a macroscopic scale one has an effective medium with eff 1.In order to be applicable the effective medium approximation, the artificial building blocks should be closely packed and much smaller than the incident wavelength.From the first condition, one deduces that interaction phenomena between the individual elements will be important and the second one suggests that the dipole approximation might give an appropriate description.The coupled electric and magnetic dipole method ͑CEMDM͒ introduced by Mulholland, Bohren, and Fuller ͓10͔ gives the right theoretical framework.It introduces in the equations of the coupled dipole method, next to the electric dipole fields, the contribution of magnetic dipoles generated by the presence of magnetic polarizable elements in the scattering structure.This is exactly the environment of a metamaterial.The magnetic polarizabilities generated by the artificial subunits are usually due to the nonsphericity of the latter.This means that generally tensorial magnetic polarizabilities should be used.This again can be handled in a straightforward manner by the CEMDM.
Furthermore, high values of eff are usually obtained by placing the subunits in resonance ͓6͔.Knowledge of the building blocks' resonance structure and how it is affected by the presence of neighboring elements is very important.A powerful tool to perform this analysis is the spectral decomposition of the interaction matrix.This method has been widely studied in case of purely electric interacting dipoles ͓11,12͔, while for electric and magnetic dipoles it has not been described yet.The objective of this paper is to fill in this gap by presenting the generalization of the eigenvalue decomposition to N interacting dipoles with both electric and magnetic polarizabilities.
On the other hand, an early theoretical study by Kerker et al. ͓13͔ showed that small spheres with 1 present unusual scattering behavior like the suppression of either the forward or backward scattering.These properties have interesting implications when one considers samples in which multiple scattering is present.Such effects were studied by Pinheiro et al. ͓14͔ who predicted, for instance, a reduction of the localization parameter kl * .Although similar studies of the resonance behavior of magnetic Mie particles were also performed ͓15͔, a detailed description of the precise interaction mechanism between magnetic scatterers is still lacking.In a previous publication ͓16͔, some effects related to the interaction of two dipolar particles with magnetic properties were discussed.For two purely electric interacting dipoles, there is an analytical solution ͓17,18͔.While for the general problem of N interacting dipoles there is no complete analytical solution, many approximations and numerical methods do exist, based on the spectral decomposition of the interaction matrix ͓11,12͔.The generalization of the eigenvalue decomposition for dipoles with both electric and magnetic polarizabilities we present is thus equally useful to study more hypothetical systems composed of small particles with 1.We apply this method to the most elementary case of N = 2 for which we give a closed mathematical form of the eigenvalues and eigenvectors.The case of two small particles is not only academic, but it is a practical situation for building nanometric rulers with interesting metrology applications in biology and nanotechnology in general ͓19͔.Fur-*merchieo@unican.es thermore, this constitutes the simplest model to analyze the physics of multiple-scattering phenomena found in more complex aggregate systems ͓16,20͔.
The paper is organized as follows.We first reformulate the coupled electric and magnetic dipole method in a more suitable form.We then explain the physical meaning of the resulting equations.Next, we present theoretical and numerical results for the case of two interacting dipoles and, finally, we draw some conclusions.

A. Electric and magnetic polarizabilities
The coupled electric and magnetic dipole method was introduced by Mulholland et al. in ͓10͔ as a generalization of the coupled electric dipole method to improve the precision of the latter by introducing the magnetic contribution to the electric field and vice versa, for particles with purely electric properties.In this situation, for spherical and isotropic subunits, the first nonzero contributions to the Mie coefficients a 1 and b 1 are given by These expressions are valid only when =1, x Ӷ 1, and ͉m͉x Ӷ 1, where m is the refractive index, x = ka, k is the wave number, and a is the particle radius.For Rayleigh particles, x is very small and b 1 is usually not considered.However, when 1, b 1 becomes of the order of x 3 : namely, In the framework of the coupled dipole method, however, one considers usually the polarizabilities of the subunits instead of the Mie coefficients a 1 and b 1 .Therefore, we will use the electric and magnetic polarizabilities ␣ e and ␣ m respectively, which are related with a 1 and b 1 coefficients in the following way: These expressions for ␣ e and ␣ m are strictly valid only in the static limit-that is, when the frequency tends to zero or, equivalently, when the incident wavelength becomes infinite.When used in the nonstatic case, it can be shown that they are not valid anymore because the optical theorem ͑which expresses the energy conservation by relating the scattered and absorbed energy with that extinguished from the forward scattered beam͒ is not satisfied ͓21͔.The extinction cross section ͑ e ͒ becomes zero if and are both real, which is impossible because scattering attenuates the intensity in the forward direction.To overcome this inconsistency of a zero value for e , ␣ e and ␣ m have to be modified.Draine ͓22͔ proposed a formalism in which the radiation reaction is introduced in the electric polarizability: where ␣ 0 e is the static polarizability as given by Eq. ͑3a͒.Other methods have also been proposed ͑e.g., Doyle ͓23͔͒ to solve this problem; see Ref. ͓24͔ for a more detailed discussion.It is not unusual to find the polarizablity written in a more generic form as and analogously for ␣ m .By giving the coefficients x and y the adequate values, one can recover the different proposed corrections like Draine's.
Until now, we have considered the electric and magnetic polarizabilities to be scalar quantities.However, the development of metamaterials presenting strong magnetic responses when excited by an electromagnetic wave has been achieved by assembling nonspherical elements such as split ring resonators or pairs of nanorods.In those cases, the geometrical configuration of the scattering elements induces magnetic dipoles and hence gives the structure a nonzero magnetic polarizability.The description of those elements requires thus generally tensor polarizabilities.Here we shall consider that the tensorial character of the polarizabilities is due to the nonspherical shapes of the scattering elements and not due to the anisotropy of the dielectric constant .That is why from now on, in order to make the treatment as general as possible, we consider the electric and magnetic polarizabilties to be tensor quantities and we write them in an equivalent manner as in Eq. ͑5͒, ␣ e = − ͑x e + iy e ͒ −1 , ͑6͒ and analogously for ␣ m .␣ e is a symmetric tensor ͓25͔, and hence x e + iy e is symmetric too.Here, the quantities x e and y e are the tensor equivalents of the values x and y.For such polarizabilities, the relation between the applied electric ͑magnetic͒ field and the induced electric ͑magnetic͒ dipole moment is given by d = ␣ e E or m = ␣ m H.

B. System of 6N coupled equations
We consider N interacting small particles located at the spatial points r 1 , ... ,r N and an incident plane wave given by where E i 0 and H i 0 are the electric and magnetic fields at the point r i and k 0 is the incident wave vector.
The electric and magnetic dipole moments of the particles are induced by the incident field on each particle and the fields scattered by all the other particles.This can be expressed by the system of equations where 0 and 0 are the dielectric constant and magnetic permeability of the vacuum, respectively.Furthermore, we have used the operators and the coefficients A ij , B ij , and D ij are given by In Eqs.͑8͒ we have assumed that the particles have different electric and magnetic polarizabilities.In Eqs.͑9͒, I ˆis the identity matrix, n ji is the unit vector pointing from particle i to particle j, denotes the tensor product, and ϫ is the vector product.In Eqs.͑10͒, r ij is the distance between particles i and j.
The system of equation ͑8͒ is linear in d i and m i and can be solved by matrix inversion to obtain the induced electric and magnetic dipole moments of each particle.Once d i and m i are known for all i, we can calculate the scattered electric field in the far zone using the relation where s is the unit vector in the scattering direction and R is the distance from the center of mass of the particles to the detector with R ӷ ʈr i − r j ʈ and kR ӷ 1.The scattering amplitude is given by ͑12͒

C. Cross sections
Using the scattering amplitude f͑s͒, the extinction ͑ e ͒, scattering ͑ s ͒, and absorption ͑ a ͒ cross sections are given by e = 4 k where s is the unit vector in the scattering direction, s 0 is the unit vector pointing to the forward scattering direction, d⍀ is the solid angle element in the s direction, E 0 is the amplitude of the incident electric field, and I 0 is the intensity of the incident field.Using Eq. ͑12͒ we can calculate these cross sections.For e , substituting the scattering amplitude f͑s͒ ͓Eq.͑12͔͒ into Eq.͑13a͒ and taking into account the property where we have used k 0 = ks 0 and the definition ͑7͒ of the incident fields.
We now perform the calculation of the scattering cross section which is slightly more involved.We substitute f͑s͒ into the equation for s , Eq. ͑13b͒, to obtain From this equation, we see that the scattering cross section can be written as a sum of three contributions: where s e and s m are the contributions of the electric and magnetic dipole moments, respectively, and s em is the cross term.
Because the integration of the first two contributions is completely analogous to that performed by Markel ͓26͔ for the electric part, we will only show the final results for these quantities.However, the third contribution is a new one.It is computed in the Appendix.The diagonal terms are given by s,ii e = 8k 4 3 and the off-diagonal terms by After summation of s,ij e over i and j with i j, we obtain a real-valued number and, consequently, the coefficients in front of Im A and Im B are real numbers too.Therefore, we can place the symbol Im in front of the summation: where ͚ i j N represents a double sum over i and j with i j.As can be seen from their definitions ͑10a͒ and ͑10b͒, the factors A͑kr ij ͒ and B͑kr ij ͒ are symmetric under permutation of i and j.If in Eq. ͑19͒, we interchange the indices, we obtain the operator C ˆij defined in Eq. ͑9a͒ multiplied on the left side by d i * and on the right side by d j .Equation ͑19͒ reduces to Analogously, for the magnetic contribution s m we have for the diagonal terms and for the off-diagonal terms For the third contribution, s em , we have for the diagonal and off-diagonal terms, respectively, where We now rewrite Eqs.͑8a͒ and ͑8b͒ as We can now substitute the relations ͑25͒ into Eqs.͑20͒ and ͑22͒.We add the two terms in Eq. ͑23b͒ for s,ij em and, after some rearrangements, obtain To find the total electric and magnetic contributions of the scattering cross section, we add the diagonal and the offdiagonal terms: where I 3 is the identity matrix with dimension 3. Following Markel ͓26͔, we now introduce the Dirac notations for which the scalar product of two vectors u and v is defined as u * • v = ͗u ͉ v͘.We write the dipole moments as and the incident fields as Furthermore, we will use the definitions where O is the zero matrix of dimension 3. Using the property Im͗d ͉ E 0 ͘ =−Im͗E 0 ͉ d͘, the electric and magnetic cross sections can then be written more compactly as To obtain an even more compact notation, we unify the electric and magnetic components in a single vector.We define Using the polarizability notation introduced in Eq. ͑6͒ and due the symmetry of ␣ e and hence of x e , one can show that ͗d ͉ x e ͉ d͘ is a real value.We can thus write Im ͳ dͯ 2 3

͑33͒
We now define the tensorial electric absorption y a e as y a e = y e − 2 3 ik 3 I 3N .͑34͒ Similarly one has the tensorial magnetic absorption y a m .The total cross section s = s e + s m + s em can then be written as where we have used O y a m ͪ.

͑36͒
We can now write the extinction, scattering and absorbtion cross sections compactly as e = 4k I 0 Im͗F 0 ͉f͘, ͑37a͒ We see that a = 0 when y a e = y a m = O, as expected.It should be noted that the above equations are only valid when ␣ m O.In the opposite case when ␣ m = O, Markel's method should be used instead, starting from the equations of the coupled electric dipole method.In this case, one can write down the cross sections as e = 4k I 0 Im͗E 0 ͉d͘, ͑38a͒ We see that both sets of equations ͑37͒ and ͑38͒ have a completely similar structure.

III. EIGENMODE ANALYSIS
In this section we discuss the relation between the resonances arising in a system of N coupled electric and magnetic dipoles and eigenvalues and eigenvectors of the interaction matrix of our scattering system.By interaction matrix we refer to the matrix which transforms the incident field into the local dipole moments.We will see that the eigenvectors of the interaction matrix correspond to the oscillation modes of the electric and magnetic dipole moments.
The analysis consist in expanding the local dipole moments ͑electric and magnetic͒ in a basis of eigenvectors of the interaction matrix.We start by rewriting the system of equations ͑8͒ using the notation we have introduced in the previous section.We define the matrices K and as G was defined in Eqs.͑30͒, and C may be expressed as where C ˆij and G ˆij are defined in Eqs.͑9͒.Then we can write the system ͑8͒ as ͉f͘ = ͉F 0 ͘ + K͉f͘.͑41͒ By solving this equation for ͉f͘, we obtain ͉f͘ = ͑I − K͒ −1 ͉F 0 ͘.͑42͒ The matrix ͑I − K͒ −1 is what we called previously the interaction matrix of the scattering system.The matrix K is not symmetric in the most general case ͑usually ␣ e ␣ m ͒; how-ever, it can be shown that the diagonalization of the non symmetric matrix K = SSB ͓where = SS and B = ͑C ,−G ; G , C͔͒ can be replaced by the diagonalization of a symmetric one: W = SBS.Even though W is symmetric, it is complex and, therefore, non-Hermitian.Its eigenvectors form a complete set, but are not orthogonal in the usual sense.
The arguments above show that diagonalization is always possible and, therefore, a basis of eigenvectors exists, though not necessarily orthogonal.We can obtain from Eq. ͑42͒ a physical interpretation of the eigenvectors.We consider an eigenvector ͉v͘ which satisfies the equation where w is an eigenvalue of K. Comparing Eqs.͑43͒ and ͑42͒, we can say that when the incident field ͉ F 0 ͘ is proportional to an eigenvector, then the induced dipole moments ͑electric and magnetic͒ are given by the same eigenvector multiplied by the complex constant ͑1−w͒ −1 .More generally, if the incident field has a non zero projection on one or more eigenvectors of the system, oscillation modes corresponding to these eigenvectors are excited.We will see later that the associated eigenvalues describe the resonant enhancements of the scattering cross section.We will express this in a more rigorous way later.
For N particles, the interaction matrix is 6N dimensional.We start again with Eq. ͑42͒, and because we have assumed that ͑I − K͒ −1 diagonalizable, we can write where D is the matrix with the eigenvalues on its diagonal and P is the matrix in which each column is an eigenvector.Because ͑I − K͒ −1 is not symmetric, the eigenvectors are not necessarily orthogonal, but they do form a basis.
To obtain a consistent notation, we write both matrices P and P −1 as where the row vectors ͗y i ͉ are the left eigenvectors of ͑I − K͒ −1 and can be written as a linear combinations of the ͗x i ͉ such that the matrix they form is the inverse of P. Equation ͑44͒ can now be written as or, equivalently, as where i are the eigenvalues associated with the eigenvector ͉x i ͘ and † is the Hermitian conjugate of .In Eq. ͑46͒, ͉f͘ is written as a linear combination of the eigenvectors ͕͉x i ͖͘ which expresses the fact that these form a basis.Equation ͑46͒expresses mathematically what we have described in a more intuitive way above.Because the eigenvectors are not orthogonal, an eigenmode is excited when the projection of the incident field on an adequate linear combination of eigenmodes is non zero or, equivalently, when the projection on the left eigenmode is non zero.The next step is now to obtain the spectral representation for the scattering cross section s .To this end we substitute Eq. ͑46͒ into Eq.͑35͒.We do this in two steps.We begin with the term ͗F 0 ͉ f͘: where we have made the substitution i = ͑1−w i ͒ −1 with w i being an eigenvalue of the matrix K.The denominators 1 − w i are complex numbers and are functions of the elements of the tensor .For some combinations of those elements, the factor ͉1−w i ͉ −2 becomes maximum and we have a resonance.So if the parameters are well chosen, the scattering cross section can exhibit a strong enhancement.We will discuss such enhancements in detail in Secs.IV and V.This is completely analogous to the purely electric case where the scattering cross section could be written as the sum of 3N resonance terms.Now, because of the magnetic component, we have 6N eigenvalues.For the absorption cross section a , we calculate the term ͗f͉Y a ͉f͘.We then obtain ͗f͉Y a ͉f͘ = ͚ i,j i j * ͗x j ͉Y a ͉x i ͗͘y i ͉͉F 0 ͗͘F 0 ͉ † ͉y j ͘. ͑49͒ We have now all the necessary mathematical tools to interpret the results.

IV. APPLICATION TO THE TWO-PARTICLE MODEL
We now discuss the case N = 2.We shall call this system a dipolar bisphere ͑DBS͒.We report explicit expressions for the eigenvalues and eigenvectors and describe how the eigenvalues are mathematically related to the scattering cross section, evaluated as a function of the interparticle distance and/or the electric and magnetic polarizabilities ͑the electricmagnetic interaction strength between the particles͒.Here, for the sake of simplicity, we will only consider scalar electric and magnetic polarizabilities.

A. Eigenvalues and eigenvectors
Because N = 2, the matrix dimension is 12ϫ 12.We have 12 eigenvalues and eigenvectors, some of them doubly degenerate.In Table I, we list all the eigenvalues and the corresponding eigenvectors.The first six components of the eigenvectors represent the electric dipole moments of particles 1 and 2, while the remaining six components are the magnetic dipole moments.This is consistent with the definition of the vector ͉f͘ in Eq. ͑32͒.As an example, we consider the eigenvector corresponding to the eigenvalue ͓1−␣ e ͑A + B͔͒: ͑0,1,0

͑50͒
Because the orientation of the eigenvectors depend on the geometrical orientation of the scattering system, we placed the DBS as shown in Fig. 1.This choice is completely arbitrary and does not alter the physical interpretation.In Fig. 2, we illustrate the relative orientations of the real parts of the electric ͑d i ͒ and magnetic eigenvectors ͑m i ͒ with respect to the DBS' orientation.The lengths of the arrows do not correspond to the magnitude of the vectors.
For comparison, we list in Table II the eigenvectors and eigenvalues in the purely electric case-i.e., when ␣ m =0.In this table, we see that the first two eigenvalues correspond to eigenvectors polarized parallel to the axis of symmetry, while the other two are eigenvectors polarized orthogonally to the axis.Consequently, we have two types of modes ͑a͒ longitudinal and ͑b͒ transverse, and each of them can be either symmetric or antisymmetric.On the other hand, when we look at the modes in Table I, we have purely electric and purely magnetic longitudinal modes, but we have also mixed electric and magnetic modes which are transverse.These eigenmodes correspond to transverse modes in the purely electric case, but they have an additional property: if their electric components are parallel, then their magnetic components are antisymmetric and vice versa.Thus, for transverse modes, we cannot speak anymore of purely symmetric or antisymmetric modes.
Comparing Tables I and II, we see that there are some common eigenvalues.These eigenvalues correspond to the longitudinal modes.The other eigenvalues in which both ␣ e TABLE I. Eigenvalues and eigenvectors of the matrix ͑I − K͒.To simplify the notation, the factors A 12 , B 12 , and D 12 are renamed, respectively, as A, B, and D.

B. Modeling the scattering cross section using the eigenvalues
The theoretical formulation of Sec.II can be used to model the spectrum of the scattering cross section as a function of the interparticle distance and electric polarizability for different values of ␣ m .The method is based on Eqs.͑37͒, ͑48͒, and ͑49͒.We first rewrite Eqs.͑48͒, and ͑49͒,

͑52͒
The scattering cross section ͓Eq.͑37͔͒ can then be separated in two parts: a first sum with the coefficients i i * and a second one with i j * .We will show that in some cases the scattering cross section can be fitted very accurately by using only the first contribution.In other cases, however, part or all of the cross terms i j * have to be used.Using Eqs.͑51͒ and ͑52͒ and after orientational averaging, the mean scattering cross section can be rewritten as We present here the results for the scattering cross section as a function of the interparticle distance r for different val-ues of and .To obtain a visible resonance spectrum, we choose = −2.013.This value corresponds to a transverse symmetric resonance mode in the case of a purely electric bisphere ͑␣ m =0͒ with r = .We keep the value of constant and change , starting with = 1 and continuing to the negative values from −1.88 down to −2.013.We chose negative values of for the same reason as we did so for ͑␣ m becomes infinite when =−2͒.This produces higher scattered intensities for a single particle as can be seen from one of the transverse electric field components in the far field region ͑kR ӷ 1͒ for a magnetic particle with x Ӷ 1 and ͉m͉x Ӷ 1 ͓13͔: where E 0 is the amplitude of the incident electric field and and are the spherical coordinates of the scattering wavevector.High scattered intensities imply a high degree of multiple scattering, which is what we are interested in.
In this study, we consider purely real optical constants and and the polarizabilities are defined by the relations ͑3͒.We use this approximation because, it has been previously used successfully to reproduce experimental results ͓24,27-29͔.

A. Scattering cross section s as a function of the interparticle distance
The results we present below have been obtained by means of the Monte Carlo simulation method.We consider a DBS freely floating in space.We generate n random orientations of the axis for a given interparticle distance and fixed optical properties, and compute the scattering cross section for each case, obtaining a mean value ͗ s ͘.We then repeat the same procedure for different interparticle distances r / and plot ͗ s ͘ as a function of r / for different values of ͑Figs.3-7͒ together with a fit of the averaged ͗ s ͘ obtained from Eq. ͑53͒.As we have mentioned above, the coefficients L i and L ij are a function of r, but to compute the fitting curves, we will consider them constant.
Due to the definition ͑34͒, we see that for real polarizabilities, we have y a e 0 and y a m 0 ͑we recall that we consider now scalar polarizabilities͒.From Eq. ͑53͒ we see that, in general, the fit will be obtained by using 36 basis functions: 8 contain single eigenvalues ͑coefficients L i ͒ and 28 contain products of the eigenvalues ͑coefficients L ij ͒.It is not always necessary to use these 36 basis functions to obtain a good fit.In many cases, we have obtained good results by considering only the 8 basis functions.These already reproduce the positions of the resonance peaks with a high degree of accuracy.
In the following results we introduce a magnetic response-i.e., 1.We start with = −1.88͑see Fig. 4͒.Again, there is a succession of resonance peaks, the one centered at r / Ϸ 1.5 being the dominant one.In the same figure, we represent the fit using 8 and 35 basis functions, which was the maximum value that could be handled by the fitting routine.The fit with 35 functions cannot be distinguished from the result obtained by Monte Carlo simulation.On the other hand, we see that for the 8 terms, the resonance peaks are very well reproduced, although differences can be clearly seen for distances r / Ͻ0.5.
The above behavior can easily be explained if one makes a plot of the individual resonance terms as functions of the interparticle distance.Then in the same way as we did for the case = 1, we can identify the terms which correspond to each of the observed peaks.For = −1.88, the peaks centered in r Ϸ n with n integer are described by the term 1+A͑−␣ e + ␣ m ͒ /2−S / 2. The peaks centered in r Ϸ͑n +1/2͒ are are described by the term 1 + A͑␣ e − ␣ m ͒ /2−S /2.Both groups of resonances have electric and magnetic contributions as can be seen from Table I.Finally, the one centered at r Ϸ 0.7 is due to the term 1 + ␣ e ͑A + B͒ which corresponds to a purely electric eigenvector.Only three terms are necessary to reproduce accurately the results of the simulation.Finally, to find a complete agreement, the 28 crossed terms have to be also inserted.In practice, however, not all of the 36 terms have always to be used.
Next, we consider the case = −1.94͑see Fig. 5͒.The observed behavior is very similar to the previous case for = −1.88 with some slight differences.The peak centered in r Ϸ 1.5 is again very prominent which indicates that the mixed eigenmode ͑with electric and magnetic components͒ corresponding to the eigenvalue 1 + A͑␣ e − ␣ m ͒ /2−S /2 is dominating.
The next value used is = −1.98 ͑see Fig. 6͒ and the highest peak has now shifted to the interparticle distance r Ϸ 2. This implies that the dominant mode is that corresponding to the eigenvalue 1 + A͑␣ m − ␣ e ͒ /2−S / 2. This mode, as can be seen from Table I, has symmetric electric components and antisymmetric magnetic components.It should be noted that, for = −1.98,we see that the fit with 8 basis functions reproduces the results of the simulation much better than for the other values of 1 we have studied.
Finally, for = = −2.103͑see Fig. 7͒, we have a behavior completely different from what has been discussed previously.It can also be observed that in this case the fit does not reproduce accurately the results obtained by simulations.In Ref. ͓13͔ it was shown that when = , the backscattered intensity is zero.Due to this phenomenon, in one special orientation of the DBS ͑axis of symmetry collinear with k 0 ͒, multiple scattering can exist only in one direction ͓16͔.
If we look now at the individual resonance terms as functions of r, we find that four terms contribute to the dominant peak centered at r = 0.7.This can be understood when we look at the mathematical structure of the eigenvalues.We have = and thus ␣ e = ␣ m , which means that the eigenvalues 1 + ␣ e ͑A + B͒ and 1 + ␣ m ͑A + B͒ are identical.The same holds for ͓1+A͑␣ m − ␣ e ͒ /2−S /2͔ and ͓1+A͑−␣ m + ␣ e ͒ /2 − S /2͔, both equal to 1 − S /2.

VI. CONCLUSIONS
In this work we present a theoretical study of the electromagnetic coupling of an ensemble of dipoles having both electric and magnetic polarizabilities.In particular, we have studied optical resonances produced by mutual interaction as a function of interparticle distances and looked at their evolution by changing magnetic properties of the particles.
In the first part of the paper we have analyzed how the formulation of the coupled dipole method, as introduced by Markel, could be generalized to electric and magnetic scatterers.The extinction, scattering, and absorption cross sections can be written in terms of the eigenvalues and eigenvectors of the generalized interaction matrix.
The fact that generally the eigenvectors are not orthogonal, but do form a basis, implies that the scattering and absorption cross sections can be written as sums of two contributions.The first is a sum of the individual eigenvalues, and the second is the sum of the double products between different eigenvalues.
We calculated the scattering cross section for a freely floating DBS as a function of the interparticle distance r.For values of and close to −2, which is the resonance value for an individual dipole, a succession of equally spaced resonance peaks was observed.We fixed the value of the electric permittivity at = −2.013,so that a symmetric transverse mode of the purely electric DBS is excited and changed the magnitude of .By changing the permeability, the resonance structure was significantly modified.
For each combination of and , a fit of the observed resonance structure was obtained, using the eigenvalue decomposition.We compared the results obtained for the fit using only the 8 eigenfunctions and the fit with the maximum number of basis functions the fitting routine could handle.The quality of the fit depends on the number of included basis function.However, for certain combinations of and , we saw that only 8 basis functions provide a very good approximation.On the other hand, for the = = −2.013, a fit as accurate as in the previous cases could not be obtained.This is probably due to the dependence of the fitting coefficients on the interparticle distance r.
The whole study was performed for configurations in which the interparticle distances are fixed.These results can be applied to situations where the distance is variable, provided that we know the distribution function of the interparticle distance and orientation.

͑A3͒
where we applied the properties of the vector product and the definition of the unit vector n ji = ͑r j − r i ͒ / r ij .The vector s is the unit vector in the scattering direction and in spherical coordinates is given by s = ͑sin cos ,sin sin ,cos ͒. ͑A4͒ We could put the factor containing d i and m i outside the integral sign because these are independent of s.We shall now consider that n ji = e z , where e z is the unit vector parallel with the z axis.This choice does not reduce the generality of the calculation, because the orientation of the system is arbitrary as long as the incident field has not to be fixed and the integration is done for each term separately.Using this reference frame, we can write s • n ji = cos .We have thus to compute the integral ͵ s exp͑ikr ij cos ͒d⍀.

͑A5͒
In reality we have three integrals because s is a vectorial magnitude.In the case i = j, all three integrals are zero.When i j, then the x and y components of this integral are zero, while the z component is nonzero and given by ͵ s z exp͑ikr ij cos ͒d⍀ = 4iͩ− cos͑kr ij ͒ kr ij + sin͑kr ij ͒ ͑kr ij ͒ 2 ͪ.

͑A6͒
Close inspection of this result shows that the expression between parentheses is the real part of k −3 D ij defined in Eq. ͑10c͒.Hence, we can write s,ij em as s,ij em = 4k

͑A7͒
s em is then obtained by summing s,ij em over i and i j.We now have to make again the substitution e z = n ji : s em = 4k

͑A8͒
We know that s em is real; thus, from the above equation we see that each factor ͓−d i ϫ m j * + m i ϫ d j * ͔ should be purely imaginary.Applying the properties of complex numbers we can rewrite this last equation as In this equation we can substitute the definition of the operator G ˆij given in by Eq. ͑9b͒.In matrix notation we then finally get s em = 4k ͑A11͒ FIG. 1. Schematic representation of the dipolar bisphere.
FIG.2.Scheme of the relative orientations of the real parts of the eigenvectors.

FIG. 3 .
FIG.3.Evolution of the mean scattering cross ͗ s ͘ section as a function of the interparticle distance r / .͗ s ͘ is the mean value for 1000 realizations.

TABLE II .
Eigenvalues and eigenvectors in the purely electric case ͑␣ m =0͒.