Clustering in particle chains - summation techniques for the periodic Green's function
Yarden Mazor Yakir Hadad Ben Z. Steinberg



School of Electrical Engineering, Tel Aviv University, Ramat-Aviv, Tel-Aviv 69978  Israel

Clustering in particle chains - summation techniques for the periodic Green's function
    Yarden Mazor[yardenm2@mail.tau.ac.il], Yakir Hadad, Ben Z. Steinberg
    
======================================================================================


1D lattice summations of the 3D Green's function are needed in many applications such as photonic crystals, antenna arrays, and so on. Such summations are usually divided into two cases, depending on the location of the observer: Out of the summation axis, or on the summation axis. Here, as a service for the community, we present and summarize the summation formulas for both cases. On the summation axis, we use polylogarithmic functions to express the summation, and Away from the summation axis we use Poisson summation (equivalent to the expansion of the field to cylindrical harmonics)


This text is not meant to be a comprehensive overview of the literature in this topic. We have included several references to selected works that incorporate parts of this overview, or other related methods. If someone feels that we have missed or did not credit his work in related matters, please do not hesitate to approach us, and we would gladly revise the bibliography. 


§ INTRODUCTION

The most general form of a cluster in a particle chain is shown in Fig. <ref>.

The theoretical "Construction" process of such a structure may be visualised as taking a N basic chains composed of particles with polarizability α_n, (n=1,...,N) and inter-particle distance D, and placing them in some general configuration, so all the chains have colinear chain axes (in this case visualised by the parallel dashed gray lines). If we assume the particles are small enough with respect to the inter-particle distance and the wavelength we may use the Discrete Dipole Approximation. In addition we assume that if as small particle with polarizability α is subject to an electric field whose local value in the absence of the particle is E^L it will respond by forming an electric dipole moment p=α E^L. if we mark the dipoles induced in each particle by p_nm where n enumerates the particles inside the cluster (n=[1,...,N]), and m will be the cluster index. 

It is worth noting before we continue, that particular cases of the formulation presented in this summary can be found in many works. <cit.> treated simple cases of 2 parallel chains, <cit.> dealt with general longitodinal particle clusters. Cases involving various summation techniques in 2D can be found in <cit.>,<cit.>,<cit.>. 

We will first formulate the most general case, and then continue by treating each of the basic "building blocks" separately. 


§ FORMULATION

 The local electric field value for each of the particles may then be written in the most general form as
	
    E^L=∑_n=1^N∑_m=-∞^∞G(r_0,r_nm)p_nm

where G is the dyadic Green's function in free space, r_0 is the location of the particle under examination, and r_nm is the location of the p_nm dipole.
Since the structure is periodic with a period D, we may express each dipole along the chain as 
	
    p_nm=p_n0e^iβ mD

Since the system is periodic, we may treat only the particles that reside in cluster 0 (noted in darker blue in Fig. <ref>). If we mark the particle we are examining in cluster 0 as the n'th particle, taking the summation in eq1 and multiplying from the left by α_n'^-1α_n' will result in the equation
	
    α_n'^-1p_n'0= ∑_n=1^N∑_m=-∞^∞G(r_0,r_nm)p_n0e^iβ mD

or in a more intuitive form
	
    α_n'^-1p_n'0= ∑_n=1^N(∑_m=-∞^∞G(r_0,r_nm)e^iβ mD)p_n0

A more accurate way to write this system of equations is by separating n=n' from the rest of the possible values of n which gives the form
	
    (∑_m=-∞
     m≠ 0^∞G̅(r_0,r_n'm)e^iβ mD-α_n'^-1)p_n'0+ 
    	∑_n=1 
     n≠ n'^N(∑_m=-∞^∞G̅(r_0,r_nm)e^iβ mD)p_n0=0

Where G̅=6πϵ_0/k^3G. Each value of n' defines a certain particle to be examined, and in fact defines a certain Base Chain we are treating. The different Base Chains from which we construct a more general Clustered Chain are noted in Fig. <ref> on the right side as C1,C2,C3,C4,C5. The last equation defines a 3N× 3N matrix equation where each value selected for n' essentially defines the lines 3n'-2, 3n'-1, 3n' of it. The matrix representing the entire system may also be described as N× N block matrix 𝔐, where each block is 3× 3  and defines the interaction between a certainbase chain to another. Using the block-matrix notation the equation may be written as 𝔐·p=0 or 
	
    [ M_11 M_12    ⋯ M_1N; M_21    ⋱         ⋮;    ⋮              ⋮; M_N1    ⋯    ⋯ M_NN ][ p_10; p_20;    ⋮; p_N0 ]
    	=0

The value of M_n',q' depends mostly on the geometrical positioning of the base chain p' in relation to the base chain of reference n'. This positioning is defined by 2 parameters  - the longitudinal shift d and the transverse shift d_0. These parameters are illustrated in Fig. <ref> for C2 and C4.
A non-trivial solution for this system exists only if
	
    det(𝔐)=0




§ THE DIAGONAL TERMS OF 𝔐

The diagonal terms arise from the left brackets in eq5 and are well-known from work done on simple particle chains (not clustered), Therefore the diagonal block M_n'n' may be written as
	
    M_n'n'=
    	[ T 0 0; 0 T 0; 0 0 L ]

where
	
    [                                                                           T                                                                           = 3/2[ 1/kDf_1(kD,β D)+ i/(kD)^2f_2(kD,β D)- 1/(kD)^3f_3(kD,β D) ] - α̅_n'^-1;                                                                            ;                                                                           L                                                                           =                  3[ - i/(kD)^2f_2(kD,β D)+ 1/(kD)^3f_3(kD,β D) ] - α̅_n'^-1;                                                                             ]

where

    f_s(x,y)=Li_s[e^i(x+y)]+Li_s[e^i(x-y)]

This block represents the interaction between a certain base chain with itself.



§ TERMS OF 𝔐 WHICH REPRESENT D_0=0

These terms represent the interaction of 2 chains that are completely co-linear with each other (Share a common chain axis). In Fig. <ref> the base chains that posses such property are C2,C3. This case is "isolated" in Fig. <ref>. 
	
 The M_n',q' term corresponding to this case has the form
	
    M_n',q'= 6π/k^3∑_mA(mD+d)e^iβ mD

Where
	
    A(z)=e^ikz/4πz[k^2A_1+(1/z^2-ik/z)A_2]

Is the dyadic Green's function simplified for the case of particles residing all on the Z-axis and A_1=(1,1,0),A_2=(-1,-1,2). In order to give a simple solution to the summation given in eq11 we are required to assume that the given inter-particle distance d is a rational fraction of the chain period D. This requires that the distances d,D are both integral multiplications of some basic distance δ
meaning
	
    D=Lδ,   d=ℓδ

Under this assumption the summation in eq11 may be rewritten as
	
    M_n',q'=h_1[ℓ]A_1-ih_2[ℓ]A_2+h_3[ℓ]A_2

where 
	
    h_s[ℓ]=3/21/(kδ)^s∑_m e^ikdmL+ℓ+iβδ mL/mL+ℓ^s.




 §.§ Evaluation of h_s[ℓ]

For 1≤ℓ≤ L-1, h_s[ℓ] can be re-written as
	
    h_s[ℓ] = 3/2e^-iβδℓ/(kδ)^s∑_m=0^∞[e^i(k+β)δ]^mL+ℓ/(mL+ℓ)^s  +  3/2e^-iβℓδ/(kδ)^s∑_m=1^∞[e^i(k-β)δ]^mL-ℓ/(mL-ℓ)^s.

We concentrate on the first sum above. It has the form
	
    σ=∑_m=0^∞(e^ix)^mL+ℓ/(mL+ℓ)^s.

This sum can be re-written as
	
    σ=∑_n=1^∞(e^ix)^n/n^s· a_n(ℓ)

where a_n(ℓ) is an auxiliary periodic sequence of period N, satisfying
	
    a_n(ℓ)={[              0, n=1,… , L,   nℓ;              1,             n=ℓ ]. .

Clearly, the series periodicity implies the recurrence relation a_n+L=a_n, whose characteristic polynomial p(λ)=λ^L-1 has N distinct roots
	
    λ^L-1=0 ⇒ {λ_r}_r=0^L-1=e^i2π r/L,

hence the infinite sequence a_n(ℓ), ∀  n can be generated by the finite sum
	
    a_n(ℓ)=∑_r=0^L-1 C_r(ℓ)λ_r^n=∑_r=0^L-1 C_r(ℓ) e^i2π r n/L.

The coefficients C_r(ℓ) can be determined using the L-initial conditions of the recurrence relation [given by eq17] in eq19. The result is the L× L Vandermonde matrix equation a(ℓ)=ΛC(ℓ). Here a(ℓ) is a vector of L entries, whose elements are given by eq15, C(ℓ) is the vector of unknown coefficients, and Λ is a Vandermonde matrix, whose n,r entry is Λ_nr=e^i2π rn/L. Due to its specific structure, Λ is also a unitary transformation from the Euclidean basis a_n(ℓ) (roam with ℓ) to a Fourier basis. Its inverse is its adjoint (normalize first). Hence
	
    C_r(ℓ)   =    L^-1e^-i2π rℓ/L ⇒
    
    	a_n(ℓ)   =   1/L∑_r=0^L-1e^i2π r (n-ℓ)/L, ∀ n.

Substituting this result into eq16 and exchanging the order of summation, we can express σ as a finite sum of Polylogarithms
	
    σ=L^-1∑_r=0^L-1e^-i2π rℓ/LLi_s(e^ix+i2π r/L).

Likewise, we may repeat essentially the same procedure for the second sum in eq14 (note the lower summation bound; shift the index by 1, and at the end change r↦ L-r'). The final result for h_s[ℓ],  1≤ℓ≤ L-1, is
	
    h_s[ℓ]   =   3/2Le^-iβδℓ/(kδ)^s∑_r=0^L-1 e^-i2π r ℓ/L f_s(kδ,βδ + 2π r/L)

and f_s are given in eq9b. Though developed for the case of 1≤ℓ≤ L-1 the expression given in eq22 is in fact valid for all values of ℓ that satisfy 1≤|ℓ|≤ L-1.
As a last remark for this section we mention that even though these formulas were developed here from "scratch", one could find many similarities to relations from signal processing where conversions of sampling rate take place.



§ TERMS OF 𝔐 WHICH REPRESENT D_0≠ 0

This case corresponds to interactions between chains as shown in Fig. <ref>.
	
This case presents us with a different challenge since the transverse shift between the base chains couples different polarizations such as X and Z polarizations. In order to properly treat this type of interactions we will write the corresponding term M_n',q' as
	
    M_n',q'=
    	[ M_n',q'^xx M_n',q'^xy M_n',q'^xz; M_n',q'^yx M_n',q'^yy M_n',q'^yz; M_n',q'^zx M_n',q'^zy M_n',q'^zz ]

Since we are discussing only 2D planar clustering of particles (the XZ plane) M_n',q'^xy=M_n',q'^yx=M_n',q'^yz=M_n',q'^zy=0.
In addition this block has to be symmetrical due to the symmetry in Green's dyad and therefore M_n',q'^xz=M_n',q'^zx. This enables to write down the block as
	
    M_n',q'=
    	[ M_n',q'^xx          0 M_n',q'^xz;          0 M_n',q'^yy          0; M_n',q'^xz          0 M_n',q'^zz ]

In the following analysis we will use the following 
	
    R_m=√(d_0^2+(mD-d)^2), r_m=(d_0,0,mD)

The main tool we are going to use to evaluate the terms in the M_n',q' block is the Poisson summation formula. this will allow us to "convert" the algebraic summations  to summations over a series of Hankel functions. since for large arguments we may obtain exponential-like decay in Hankel functions, this allows us to sum a finite number of elements in the series and still obtain a very good approximation for the terms in the matrix.


 §.§ Evaluation of M_n',q'^xx

The value of the x component of E^L that the red chain creates at the origin
	
    E_x^L(r=(0,0,0))=∑_m=-∞^∞G_xx(0,r_m)p_0xe^iβ mD

where p_0xe^iβ mD are the X components of the dipole moments induced on the m'th particle in the red chain.
substituting G_xx we may obtain an explicit form for M_n',q'^xx
	
    M_n',q'^xx=3/2∑_m=-∞^∞e^ikR_m[1/kR_m-(kd_0)^2/(kR_m)^3-1/(kR_m)^3+i1/(kR_m)^2+3(kd_0)^2/(kR_m)^5-3i(kd_0)^2/(kR_m)^4]e^iβ mD

For easier evaluation this may also be written as
	
    M_n',q'^xx=3/2∑_m=-∞^∞[e^ikR_m/kR_m+∂^2/∂ (kd_0)^2e^ikR_m/kR_m]e^iβ mD

We start by evaluating the summation over the left term in eq29. Substituting R_m we get
	
    3/2∑_m=-∞^∞e^ikR_m/kR_m=3/2∑_m=-∞^∞e^ik√(d_0^2+(mD-d)^2)/k√(d_0^2+(mD-d)^2)

To evaluate this sum we use Poisson summation formula
	
    ∑_m=-∞^∞f(m)=∑_n=-∞^∞∫_-∞^∞f(x')e^-2π inx'dx'

Substituting our expressions into eq31 we obtain
	
    3/2∑_m=-∞^∞e^ik√(d_0^2+(mD-d)^2)/k√(d_0^2+(mD-d)^2)e^iβ mD=3/2∑_n=-∞^∞∫_-∞^∞e^ik√(d_0^2+(x'D-d)^2)/k√(d_0^2+(x'D-d)^2)e^iβ x'De^-2π inx'dx'

Evaluation of such an integral is possible using integral tables
	
    ∫_-∞^∞e^ik√(d_0^2+(x'D-d)^2)/k√(d_0^2+(x'D-d)^2)e^iβ x'De^-2π inx'=
    	-e^(β/k-2π n/kD)kd/kDπ/iH_0^(1)[kd_0√(1-(2π n/kD-β/k)^2)]

we may define the normalized parameters for easier notation
	
    d̅=kd,  d̅_0=kd_0,  D̅=kD,  β̃_n=2π n/D̅-β/k

and finally obtain
	
    3/2∑_m=-∞^∞e^ikR_m/kR_m=-3/2∑_n=-∞^∞e^-iβ̃_nd̅/D̅π/iH_0^(1)[d̅_0√(1-β̃_n^2)]

We may treat the summation over the right term in eq29 using the same method and applying the needed derivatives we obtain
	
    3/2∑_m=-∞^∞∂^2/∂ (kd_0)^2e^ikR_m/kR_m=3/2∑_n=-∞^∞e^-iβ̃_nd̅/D̅π/i(1-β̃_n^2)·
    ·1/2(1-β̃_n^2){H_0^(1)[d̅_0√(1-β̃_n^2)]-
    	H_2^(1)[d̅_0√(1-β̃_n^2)]}

Adding eq35 and eq36 we obtain the final expression
	
    M_n',q'^xx=3/2∑_n=-∞^∞e^-iβ̃_nd̅/D̅π/i·{-1/2(β̃_n^2+1)H_0^(1)[d̅_0√(1-β̃_n^2)]+
    	1/2(β̃_n^2-1)H_2^(1)[d̅_0√(1-β̃_n^2)]}

Where H^(1)_n are Hankel functions of the first kind, of order n. We will repeat a very similar process for the other M terms in the block.



 §.§ Evaluation of M_n',q'^yy

The evaluation of M_n',q'^yy is somewhat less tedious since y polarization is not coupled with other polarizations. The value of the y component of E^L that the red chain creates at the origin due to the y components of the dipoles in the red chain
	
    E_y^L(r=(0,0,0))=∑_m=-∞^∞G_yy(0,r_m)p_0ye^iβ mD

where p_0ye^iβ mD are the y components of the dipole moments induced on the m'th particle in the red chain.
substituting G_yy we may obtain an explicit form for M_n',q'^yy
	
    M_n',q'^yy=3/2∑_m=-∞^∞e^ikR_m[1/kR_m+i/(kR_m)^2-1/(kR_m)^3]e^iβ mD

This may also be written as 
	
    M_n',q'^yy=3/2∑_m=-∞^∞[e^ikR_m/kR_m+1/d̅_0∂/∂d̅_0e^ikR_m/kR_m]e^iβ mD

Summation over the left term gives the result presented in eq35. After performing the derivatives summation over the right term gives 

    3/2∑_m=-∞^∞1/d̅_0∂/∂d̅_0e^ikR_m/kR_m=
    3/2∑_n=-∞^∞e^-iβ̃_nd̅/D̅π/i√(1-β̃_n^2)/d̅_0
    H_1^(1)[d̅_0√(1-β̃_n^2)]

Adding eq35 and eq41
	
    M_n',q'^yy=3/2∑_n=-∞^∞e^-iβ̃_nd̅/D̅π/i·{-H_0^(1)[d̅_0√(1-β̃_n^2)]+
    	√(1-β̃_n^2)/d̅_0H_1^(1)[d̅_0√(1-β̃_n^2)]}




 §.§ Evaluation of M_n',q'^zz

The value of the z component of E^L that the red chain creates at the origin due to the z components of the dipoles in the red chain 
	
    E_z^L(r=(0,0,0))=∑_m=-∞^∞G_zz(0,r_m)p_0ze^iβ mD

where p_0ze^iβ mD are the z components of the dipole moments induced on the m'th particle in the red chain.
substituting G_zz we may obtain an explicit form for M_n',q'^zz
	
    M_n',q'^zz=3/2∑_m=-∞^∞e^ikR_m[1/kR_m-(d̅-mD̅)^2/(kR_m)^3-1/(kR_m)^3+i1/(kR_m)^2.
    .+3(d̅-mD̅)^2/(kR_m)^5-3i(d̅-mD̅)^2/(kR_m)^4]e^iβ mD

This may also be written as 
	
    M_n',q'^zz=3/2∑_m=-∞^∞[e^ikR_m/kR_m+1/d̅_0∂^2/∂(d̅-mD̅)^2e^ikR_m/kR_m]e^iβ mD

Again, Summation over the the left term is given by eq35. Using Fourier transform properties, we may write an expression for the summation over the right term
	
    3/2∑_m=-∞^∞[e^ikR_m/kR_m+1/d̅_0∂^2/∂(d̅-mD̅)^2e^ikR_m/kR_m]e^iβ mD=
    	3/2∑_n=-∞^∞e^-iβ̃_nd̅/D̅π/iβ̃_n^2H_0^(1)[d̅_0√(1-β̃_n^2)]

Substituting the results into eq45 we get
	
    M_n',q'^zz=3/2∑_n=-∞^∞e^-iβ̃_nd̅/D̅π/i(β̃_n^2-1)H_0^(1)[d̅_0√(1-β̃_n^2)]




 §.§ Evaluation of M_n',q'^xz

The value of the z component of E^L that the red chain creates at the origin due to the x components of the dipoles in the red chain 
	
    E_z^L(r=(0,0,0))=∑_m=-∞^∞G_xz(0,r_m)p_0xe^iβ mD

where p_0xe^iβ mD are the z components of the dipole moments induced on the m'th particle in the red chain.
substituting G_xz we may obtain an explicit form for M_n',q'^xz
	
    M_n',q'^xz=3/2∑_m=-∞^∞e^ikR_m[-d̅_0(d̅-mD̅)/(kR_m)^3-3id̅_0(d̅-mD̅)/(kR_m)^4+3d̅_0(d̅-
    	mD̅)/(kR_m)^5]e^iβ mD

eq49 can also be presented as
	
    M_n',q'^xz=3/2∑_m=-∞^∞[∂/∂d̅_0∂/∂(d̅-mD̅)e^ikR_m/kR_m]e^iβ mD

Using previous results obtained we get
	
    M_n',q'^xz=3/2∑_n=-∞^∞e^-iβ̃_nd̅/D̅π/i(iβ̃_n√(1-β̃_n^2))H_1^(1)[d̅_0√(1-β̃_n^2)]



 §.§ How many Hankel function terms shuold we sum?

The Hankel function summations we have received are composed of terms of the general form

    A_n=CQ(β̃_n)H_ν^(1)[d̅_0√(1-β̃_n^2)]

Where Q may be a polynomial of order up to 2 or a combination of a polynomial with the square root function (Still, the maximal power of β̃_n in Q is 2), and C is some complex constant. For the cases examined in sections  <ref> through  <ref> the "worst-case scenario" for the decay of the terms in the Hankel function series is when Q is a polynomial of order 2. for this case the series will have terms such as
	
    A_n∝β̃_n^2H_ν^(1)[d̅_0√(1-β̃_n^2)]

or more explicitly
	
    A_n∝(2π n/D̅-β/k)^2H_ν^(1)[d̅_0√(1-(2π n/D̅-β/k)^2)]


We may define the criteria for the number of terms to use for the summation as n_0 such that the tail of the summation ∑_n=n_0^∞A_n will be significantly smaller then a threshold constant C.
For large values of n we may use the large argument approximation for Hankel's function

    H_ν^(1)(z)∼√(2/π z)e^i(z-1/2νπ-1/4π)

If we subsitute this approximation into the required sum we obtain

    |∑_n=n_0^∞A_n|≤∑_n=n_0^∞|A_n|∼∑_n=n_0^∞π/D̅n^3/2/√(2πd_0/D)e^-2d_0/Dπ n<C

from eq56 it is easy to see that the decay rate of the summation is strongly dependant on the ratio d_0/D. Evaluating the given summation may prove challenging, but if we don't mind giving a more strict criteria we may evaluate the needed n_0 from the following

    ∑_n=n_0^∞π/D̅n^3/2/√(2πd_0/D)e^-2d_0/Dπ n<∑_n=n_0^∞π/D̅n^2/√(2πd_0/D)e^-2d_0/Dπ n<C

And the required value of n_0 may be extracted from the relation

    π/D̅1/√(2πd_0/D)e^-2πd_0/D(n_0+1)/(e^-2πd_0/D-1)^3[n_0^2e^-4πd_0/D+(-2n_0^2-2n_0+1)e^-2πd_0/D+(n_0+1)^2]<C



§ EXAMPLE - EXAMINATION AND VERIFICATION OF THE MAGNETIC MODEL FOR PARTICLE RINGS


In <cit.> a magnetic polarizability model for particle rings is developed. Each ring is composed of N spherical particles of radius a, positioned on the circumference of a ring with radius R with uniform angular shift from each other. The electric dipoles induced on each particle in this case are in a direction tangent to the ring. Example for the setup for such a particle ring with 6 particles which we will use for the example is given in Fig. <ref>. 

When exposed to y-polarized magnetic field H_0ŷ, an equivalent magnetic dipole moment will be induced on the particle ring that satisfies the relation α_mH_0=m_0 where the inverse magnetic polarizability is given by

    4π/k^3α_m^-1=8/3N(kR)^2[6πϵ_0/k^3α_p^-1]
    -1/4N(kR)^5∑_m=1^N-1e^2ikRsinmπ/N/sin^3(mπ/N)Q_m

Where

    Q_m={3-3(kR)^2+(1+4(kR)^2)cos(2mπ/N)-.
    .-kR(kRcos(4mπ/N) +i[5sin(mπ/N)+sin(3mpi/N)])}

And α_p is the electric polarizability of the particles composing the ring. The dispersion curve for a chain of such rings may be calculated in two ways. The first, is by taking the equivalent magnetic polarizability of such a ring, and assuming the chain is a simple chain composed of magnetic particles that have the given polarizability. This way is very simple and streightforward. The second, by taking a particle-by-particle model of such a ring, using the formulation developed here previously. Ofcourse the second method contains many more degrees of freedom and the dispersion curve will most certainly contain many branches, but we expect one of these branches to be consistent with the model in <cit.>, and moreover, the nullspace vector for that branch will exhibit tangent vectors, that represend the tangent dipoles. The chain setup is demonstated in Fig. <ref>.

Dispersion curves for such chain of rings is given in Fig. <ref>.  Its very clear that the 2 methods provide dispersion curves that overlap. calculating the nullspace vector for some selected solution from the dispersion curves shows dipole vectors tangent to the particle ring, which proves that this curve represents magnetic mode that's propagating along the chain. the parameters used for the examined chain: D=λ_p/20, R= λ_p/200, a=λ_p/800.

99
AluParallel
Andrea Alu, Pavel A Belov and Nader Engheta, "Coupling and guided propagation along parallel chains of plasmonic nanoparticles", New Journal of Physics, Volume 13, (2011)

LongCh
Yarden Mazor, Ben Z. Steinberg, "Longitudinal chirality, enhanced nonreciprocity, and nanoscale planar one-way plasmonic guiding", Phys. Rev. B 86 (4), 045120



TretyakovBook
S. Tretyakov Analytical modeling in applied electromagnetics, Artech House (2003).

Cappolino_PRE_2011 
S. Steshenko, F. Capolino, P. Alitalo, and S. Tretyakov, Phys. Rev. E, 84, 016607 (2011).

Vitaly
D. Van Orden and V. Lomakin, "Rapidly Convergent Representations for Periodic Green's Functions of a Linear Array in Layered Media," IEEE Transactions on Antennas and Propagation, vol. 60, no. 2, pp. 870-879, (2012).



EnghetaRings
Andrea Alu and Nader Engheta, "Dynamical theory of artificial optical magnetism produced by rings of plasmonic nanoparticles", Phys. Rev. B, 78, 8 085112, (2008)



