## Abstract

We report a semi-analytical model for calculating the coupling effects between the dipolar surface plasmon nanoparticles of a periodic structure. This model involves real-valued frequencies only and is therefore applicable to periodic structures with arbitrary dipolar units and environments.

©2008 Optical Society of America

## 1. Introduction

Periodic metallic nanoparticle structures have important potential applications as miniaturized waveguides that carry signals by surface plasmons and for the enhancement and tuning of localized surface plasmon resonances for biological and chemical sensing purposes [1–13]. Study of the coupling interactions between the surface plasmon nanoparticles of periodic structures is critical in order to understand the signal transfer and resonance properties. Theoretical and experimental studies have been reported on such coupling interactions in a uniform dielectric environment and the resulting dispersion and loss properties of the periodically coupled surface plasmon modes [2,14–17]. In a previous work, we reported measurement and numerical modeling of the dispersion and extinction properties of transversely polarized surface plasmons in a 2-D periodic metallic nanoparticle array at the interface between air and glass half spaces [18]. The reported dispersion relation was significantly different from those of surface plasmons in a uniform dielectric environment [14,16,17].

Previous theoretical studies employed the coupled dipole method (CDM) and found the resonance properties by solving the equation **p**(ω)=α(ω)**E**(ω) for its complex eigenvalues ω [2,14,16]. **p** is the dipole moment of the surface plasmons in one particle, **E** is the electric field incident on this particle, *α* is the single particle polarizability, and ω is the angular frequency of the surface plasmon. The real parts of the eigenvalues are taken as the resonance frequencies; the imaginary parts of the eigenvalues are taken as the decay rates. In these previous works, the value of α(ω) is that of a nano sphere in a uniform dielectric environment, which can be analytically obtained at arbitrary complex frequencies (ω). However, experiments or first principle simulations, such as the finite difference time domain (FDTD) method, are often only efficient at obtaining α values at real frequencies. Therefore, an alternative theoretical approach is needed to model many other realistic situations. In this paper, we describe a different CDM method which only uses α values at real frequencies. The dispersion and extinction properties are obtained by this semi-analytical model for the surface plasmons in a 2-D periodic metallic nanoparticle array at an air/glass interface. The semi-analytical modeling results compare well with previous modeling carried out completely by FDTD [18]. We also report modeling results for a 1-D chain of metallic nanoparticles at an air/glass interface. All the modeling results reported in this paper are for one kind of transversely polarized surface plasmon, where the dipole moments are oriented along the air/glass interface and perpendicular to the wave vector; this statement won’t be repeated in the following part of the paper. The proposed method is expected to be applicable to the other surface plasmon modes as well.

## 2. The semi-analytical model

When studying the coupling interactions between different metallic nanoparticles, the surface plasmons in each particle are approximated by a point dipole. The dipolar polarizability of each nanoparticle, α, is obtained numerically by FDTD. The field from each point dipole is also numerically obtained by FDTD. The coupling between each point dipole is then calculated analytically. For situations in which the α values and dipole fields can be obtained analytically, like in the previous theoretical works [14,16], this model can be completely analytical. By this model, once the resonance properties and field profiles for a single unit cell have been obtained, the resonance properties at any wave vector, i.e., the whole dispersion relation, can be obtained analytically for any periodic arrangements of the same unit cell structure. By comparison, in our previous work carried out completely using the FDTD method [18], the resonance properties at each wave vector had to be calculated independently by FDTD.

#### 2.1. Point dipole approximation

We adopt the point dipole approximation, which has also been used by previous theoretical works [2,14,16]. It was found in ref [15] that when the center-to-center distance between nearest nanoparticles is larger than 1.5 times the particle diameter, the surface plasmons in each particle can be approximated by a point dipole when studying the coupling between particles. The structure we have modeled and experimentally characterized satisfies this requirement. Figure 1 is a scanning electron micrograph of the sample, with the definitions of the *(x,y,z)* axes also shown. It is an array of gold nanocylinder chains on a glass substrate. Each chain consists of gold cylinders spaced by 140 nm in the *x* direction. Each cylinder has a diameter of ~90 nm and a height of 55 nm. Adjacent chains are spaced by 300 nm in the y direction. Details of the fabrication processes can be found in Reference [19].

The values of α for a single nanoparticle for a range of real-valued frequencies are obtained by comparing the field scattered by the particle with the field generated by a point dipole, under pulsed conditions. The particle scatters a normally incident plane wave. Figure 2 is a schematic explanation of the modeling. α is defined as the equivalent point dipole moment divided by the incident electric field: α(ω)=[scattered field (ω)/point dipole generated field (ω) × point dipole moment (ω)]/incident electric field (ω). The scattered field and point dipole generated field values are sampled at (*x,y,z*)=(0, 0, 1 µm), i.e. above the nanoparticle and the point dipole. The point dipole is positioned at (0, 0, 20nm), a short distance from the air/glass interface in order to more closely approximate the surface plasmons in the 55nm tall cylinder shaped nanoparticle. The incident electric field is taken at the position (*x,y,z*)=(0, 0, 0), and is the sum of the incoming field and the field reflected by the air/glass interface. The scattered field is obtained by subtracting the incident field, i.e. the field in the absence of the particle, from the total field. Figure 3 shows α calculated using this approach for a range of real-valued frequencies for the structure of Fig. 1. It should be noted that we account for the periodicity of the structure by considering it a one-dimensional periodic arrangement of unit cells spaced by 140 nm along the *x* direction. Each unit cell consists of a line of particles spaced by 300 nm along the *y* direction, and, for the equivalent point dipole moment calculation, is equated to a line of identical point dipoles. Therefore, the α values obtained are the polarizability values of the unit cell. In the following part of this paper, we will describe our model as if the unit cell is one single nanoparticle, for the simplicity of description. It should be noted that the method we introduce is applicable to any types of unit cells that consist of uniformly excited nanoparticles. In Section 3, we will apply the method to a unit cell consisting of a single nanoparticle.

To confirm that the point dipole approximation is valid, we compare the field scattered by the nanoparticle to the field generated by the point dipole by normalizing the field values at two positions along the chain to the field at (*x,y,z*)=(0, 0, 1 µm). It is appropriate to normalize to the field at this point as the polarizability is found by equating the scattered field at this point to the field at this point generated by a point dipole. The two other positions are (*x,y,z*)=(140nm, 0, 0), which corresponds to the nearest nanoparticle in the chain, and (*x,y,z*)=(980 nm, 0, 0), which corresponds to seven periods away along the chain. As shown in Fig. 4, the nanoparticle scattered field profile and the point dipole generated field profile are in good agreement, which confirms the accuracy of the point dipole approximation for the frequency range of interest.

#### 2.2. The model

Let the equivalent point dipole moment of the plasmon oscillation in a single nanoparticle be denoted by *P*. Writing the plasmon oscillation in the form of the equation of motion of a harmonic oscillator, we have:

where *M* is the total effective mass of all the electrons in the nanoparticle that participate in the plasmon oscillation, *t* is time, *K _{0}* is the Hooke’s coefficient for the restoring force and is a real number,

*Q*is the total charge of all the electrons in the nanoparticle that participate in the plasmon oscillation, and

*E*is the external incident electric field.

_{o}*ν*represents all loss mechanisms including ohmic and radiation losses. Though

*ν*should contain time derivatives of P to include the radiation loss, we will simply treat it as a real number under monochromatic condition. Coupling between particles is included in Equation (1) through the quantity

*E*, which is a sum of electric fields generated by plasmons in all the other nanoparticles at the position of the particle being considered.

_{i}*E*is a function of the oscillation angular frequency

_{i}*ω*and the wave vector

*k*of the surface plasmon mode (

_{x}*k*is along the direction of chain:

_{x}**x**). As we shall discuss, a real-valued artificial gain coefficient

*G*for the dipole moment P: dP/dt=GP, is introduced for the purpose of defining resonance. Denoting every oscillating quantity

*X*in Equation (1) by

*X*, we rewrite Equation (1) in the frequency domain:

_{0}e_{-iωt}where
${\omega}_{0}=\sqrt{\frac{{K}_{0}}{M}}$
is the surface plasmon resonance frequency of an isolated single particle, and K_{1}≡E_{i0}/P_{0}. Later we will calculate the coupling term K_{1} and assume a linear relationship between E_{i} and P at any single (ω, k_{x}) point. Thus in this model G is the gain coefficient for both the dipoles and the fields. Therefore G is the gain coefficient for the whole surface plasmon mode. Note that if we rewrite Equation (2) for an isolated single particle, i.e., with no coupling or artificial gain terms, we obtain the polarizability of a single metallic particle, P_{0}/E_{o0}, which has the same form as the polarizability of a single metallic nanosphere in a uniform dielectric environment under Drude’s model as have been used in References [14,16]. For metallic nanospheres in a uniform dielectric environment, radiation loss is included in ν by including a ω^{2} term [14,16]. The left hand side of Equation (2) is separated into real and imaginary parts.

In this paper, except in section 4, we define the resonance frequency *ω* of the surface plasmon mode to be a real value *ω _{r}* such that there exists an appropriate gain coefficient, called the threshold gain

*G*, for the mode to oscillate at

_{th}*ω*for an infinite length of time without being pumped by external fields (E

_{r}_{o}=0). This is the same definition as lasing. If we remove the artificial gain term in Equation (2) and look at the response of dipole moment P to external driving electric field E

_{o}, we see that this definition of resonance frequency is also the same as the peak frequency of the |P

_{0}|~|E

_{o0}| response, under the assumption that change in the imaginary part of Equation (2) can be ignored near resonance, which is true for high quality factor modes. However this definition is a little different from previous theoretical works where the resonance frequency is defined as the center oscillation frequency of a freely decaying surface plasmon resonance [2,14,16]. To understand this difference, the freely decaying case will be further discussed in section 4.

Let us now find the resonance conditions from Equation (2). Under zero external field *E _{o}*, the left hand side of Equation (2) is equal to zero. Thus the real part of

*K*determines the shift in resonance frequency

_{1}*ω*; and the imaginary part of

_{r}*K*changes the threshold gain

_{1}*G*required for a non-decaying oscillation. Because the value of

_{th}*K*is independent of the value of

_{1}*G*, the value of

*ω*is determined by the condition that the real part of

_{r}*E*/

_{o}*P*be zero, regardless of the value of

*G*. This can be physically explained by noting that when the real part of

*E*/

_{o}*P*is zero,

*E*is

_{o}*π*/

*2*out of phase with respect to

*P*, meaning that the external field does not change the strength of the restoring force and the surface plasmon oscillates at its intrinsic frequency. For this reason, in our modeling,

*G*is only incorporated to clarify the definition of resonance but not included in the actual calculation. Letting

*G*=0, we obtain the value of

*ω*by finding when the real part of

_{r}*E*/

_{o}*P*equals zero. The imaginary part of

*E*/

_{o}*P*is proportional to

*ωG*under the assumption that the number of electrons participating in the plasmon oscillation is invariant. In section 4 it will be shown that

_{th}*G*is twice the decay rate of the freely decaying surface plasmons. As an aside, the existence of

_{th}*Q*terms in Equation (2) indicates larger coupling effects for larger nanoparticles [9].

_{2}We now consider the calculation of *E _{o}*/

*P*. Because

*P*=

*α(E*=

_{o}+E_{i})*α(E*, we have

_{o}+K_{1}P)*E*/

_{o}*P*=

*1*/

*α-K*. The resonance frequency is therefore found by finding the frequency at which

_{1}*Re(1/α-K*=

_{1})*0*. The corresponding threshold gain

*G*is found from

_{th}*Im(1/α-K*. The value of

_{1})*α*is found as discussed in section 2.1. The value of

*K*is found by a weighted sum of the fields generated by single nanoparticles, with the weights dependent on the wave vector

_{1}*k*[18]. Three FDTD simulations are performed in order to find

_{x}*α*and

*K*, corresponding to the three plots in Fig. 2: the incident field, the sum of incident and scattered fields, and the point dipole field. For nanospheres in a uniform dielectric environment, it would be possible to replace the FDTD simulations by analytical calculations [14,16]. Figure 5 shows the values of K

_{1}_{1}calculated for k

_{x}values spanning the first Brillouin Zone, for the 2-D metallic nanoparticle array at an air/glass interface shown in Fig. 1. The frequency is fixed at

*f*/

*c*=

*1.6 µm*. It presents a fast phase change between the light lines of air and glass, which has been reported and attributed to the phase change of total internal reflection [18]. In Fig. 6(a), we present calculations of |

_{-1}*Re(1/α-K*| over a range of frequencies and wavevectors spanning the first Brillouin Zone. As discussed, the resonance frequencies are found from the zeros of |

_{1})*Re(1/α-K*|. These are presented as the dispersion relation of Fig. 6(b), along with the threshold gain

_{1})*G*.

_{th}#### 2.3. Comparison with modeling carried out completely by FDTD

In order to evaluate the accuracy of our model, we employ it to calculate the extinction ratio of the 2-D metallic nanoparticle array at an air/glass interface as in Fig. 1, under plane wave illumination. We compare the results with previous modeling carried out using FDTD exclusively [18]. The extinction ratio is defined here as the ohmic loss divided by the incident power. The calculation starts with obtaining the equivalent excited point dipole moment *P* for the incident plane wave amplitude *E _{o}*, using

*E*/

_{o}*P*=

*1*/

*α-K*. The relationship between ohmic loss and |

_{1}*P*|

*is found by FDTD simulations of (*

^{2}*1-T-R*) under normal incidence illumination (

*k*=

_{x}*0*), where

*T*and

*R*are the power transmission and reflection coefficients, respectively. Fig. 7a shows the extinction ratio obtained with this model over a range of frequencies and wave vectors. Fig. 7b shows the results of the modeling carried out completely by FDTD for comparison. In the FDTD-only approach of Fig. 7b, FDTD simulations were done for each of ~30 k

_{x}values [18].

From Fig. 7 it can be seen that the results of the semi-analytical model are in very good agreement with the modeling carried out using FDTD exclusively. There are some potential sources of inaccuracy, however. Inaccuracy of this model can come from the fact that we have employed a universal α value, which is obtained for the plane wave normal incidence case as in Fig. 2, while the actual *α* value depends on the incident field’s direction and profile. In addition, the point dipole approximation may not be sufficiently accurate. As frequency and wave vector change, the plasmon profile within the metallic particle is also different, meaning that the number of electrons participating in the plasmon oscillation is different, which has been assumed to be invariant when calculating *G _{th}* in Fig. 6(b).

## 3. Discussion of modeling results

Comparing Fig. 5 and Fig. 6, it is observed that the fast decrease of resonance frequency between the light lines of air and glass corresponds to the fast increase of the real part of the coupling coefficient *K _{1}*. The fast decrease of

*G*, or twice the free decay rate, in the same region corresponds to the fast decrease of the imaginary part of

_{th}*K*.

_{1}*G*peaks up on the light line of air because of maximum coupling to the radiation mode; its fast decrease between the light lines of air and glass is because of decrease of the same coupling effect. The strength of this coupling depends on the amplitude of the radiation mode at the air/glass interface. This amplitude is determined by the plane wave transmission coefficient between glass and air, which peaks up on the light line of air and decreases to zero as

_{th}*k*increases to the light line of glass [18]. Below the light line of glass, radiation loss diminishes and

_{x}*G*is completely due to the metal ohmic loss.

_{th}It is worth pointing out that the fast phase change of *K _{1}* between the light lines of air and glass does not necessarily correspond to fast changes in both of its real and imaginary parts. Fig. 8(a) presents

*K*for a single chain of gold nanoparticles at an air/glass interface. This chain is the same as each chain in Fig. 1. Only the imaginary part of

_{1}*K*has a fast change between the light lines of air and glass. Figure 8(b) shows the modeled dispersion and

_{1}*G*behaviors. It is not obvious that the dispersion is significantly larger in the region between the light lines of air and glass.

_{th}In Fig. 6(b), the dispersion relation curve drops abruptly near the light line of air. This corresponds to the fast continuous increase of *Real(K _{1})* in Fig. 5, but not a jump in

*K*or an anti-crossing split between the plasmon dominated and photon dominated polariton modes as in ref [16]. In our model, only the plasmon dominated mode will be observed because all the electric fields that have been considered are generated from plasmons. At this abrupt yet continuous drop, |

_{1}*dω*/

_{r}*dk*| is larger than the speed of light (fast light), which may seem paradoxical. This is because above the light line of glass the photon part of the polariton does not propagate along the

_{x}*x*-axis and

*dω*/

_{r}*dk*is not a correct representation of the group velocity

_{x}*dω*/

_{r}*d*[20].

**k**## 4. Decay rate of periodically coupled surface plasmons

In the previous theoretical studies, by solving the equation **p**(ω)=α(ω)**E**(ω) for its complex eigenvalues ω, the center oscillation frequencies and the decay rates of freely decaying surface plasmon modes were obtained [2,14,16]. In this section, we solve for the same properties using our model and compare with the preceding sections where we use a gain term for non-decaying oscillations.

Denoting every oscillating quantity *X* in Equation (1) by *X _{0}e^{-(iω+s)t}*, and ignoring the gain G and the external driving field E

_{o}, we rewrite Equation (1) in the frequency domain with separated real and imaginary parts:

The resonance frequency ω is obtained from equating the real part of Equation (3) to 0; the free decay rate s is obtained from equating the imaginary part of Equation (3) to 0. Compared to Equation (2), the resonance frequency thus defined is different from preceding sections by an amount that shifts with the decay rate s. By comparing Equations (2) and (3), it can be seen that the decay s=G/2. Here we have assumed that the mode profile has not changed, i.e., the parameters in Equation (2) and (3) are the same, which is not accurate if the decay rate is high. Another factor that we have ignored is that *ν* is no longer a real number due to radiation loss at complex frequencies, which will significantly change the result if the decay rate is high. As long as the decay rate, or the threshold gain, is small enough, the results obtained from the two different definitions are close to each other.

## 5. Summary

We have demonstrated a semi-analytical model to calculate the periodic coupling between dipolar surface plasmons of metallic nanoparticles. In this model the single nanoparticle polarizability values and the point dipole field values are obtained by FDTD simulations; then the periodic coupling effects is calculated analytically. The model only involves real-valued frequencies and is appropriate for general periodic structures composed of dipolar units, while it requires only a few number of FDTD simulations to obtain the whole dispersion diagram. The modeled dispersion and loss properties have been reported and physically explained. The modeled extinction properties agree well with full FDTD simulation results.

## Acknowledgements

We thank the Defense Advanced Research Projects Agency (DARPA), the Charles Stark Draper Laboratory and the Harvard Nanoscale Science and Engineering Center (NSEC) for the financial support of this work. The Harvard NSEC is supported by the National Science Foundation (NSF).

## References and links

**1. **M. Quinten, A. Leitner, J. R. Krenn, and F. R. Aussenegg, “Electromagnetic energy transport via linear chains of silver nanoparticles,” Opt. Lett. **23**, 1331–1333 (1998). [CrossRef]

**2. **M. L. Brongersma, J. W. Hartman, and H. A. Atwater, “Electromagnetic energy transfer and switching in nanoparticle chain arrays below the diffraction limit,” Phys. Rev. B **62**, R16356–16359 (2000). [CrossRef]

**3. **S. A. Maier, P. G. Kik, H. A. Atwater, S. Meltzer, E. Harel, B. E. Koel, and A. A. G. Requicha, “Local detection of electromagnetic energy transport below the diffraction limit in metal nanoparticle plasmon waveguides,” Nature Mater. **2**, 229–232 (2003). [CrossRef]

**4. **N. Félidj, J. Aubard, G. Lévi, J. R. Krenn, A. Hohenau, G. Schider, A. Leitner, and F. R. Aussenegg, “Optimized surface-enhanced Raman scattering on gold nanoparticle arrays,” Appl. Phys. Lett. **82**, 3095–3097 (2003). [CrossRef]

**5. **J. R. Krenn, A. Dereux, J. C. Weeber, E. Bourillot, Y. Lacroute, and J. P. Goudonnet, “Squeezing the optical near-field zone by plasmon coupling of metallic nanoparticles,” Phys. Rev. Lett. **82**, 2590–2593 (1999). [CrossRef]

**6. **B. Lamprecht, G. Schider, R. T. Lechner, H. Ditlbacher, J. R. Krenn, A. Leitner, and F. R. Aussenegg, “Metal Nanoparticle Gratings: Influence of Dipolar Particle Interaction on the Plasmon Resonance,” Phys. Rev. Lett. **84**, 4721–4724 (2000). [CrossRef] [PubMed]

**7. **S. A. Maier, M. L. Brongersma, P. G. Kik, and H. A. Atwater, “Observation of near-field coupling in metal nanoparticle chains using far-field polarization spectroscopy,” Phys. Rev. B **65**, 193408 (2002). [CrossRef]

**8. **S. A. Maier, P. G. Kik, and H. A. Atwater, “Observation of coupled plasmon-polariton modes in Au nanoparticle chain waveguides of different lengths: Estimation of waveguide loss,” Appl. Phys. Lett. **81**, 1714–1716 (2002). [CrossRef]

**9. **Q.-H. Wei, K.-H. Su, S. Durant, and X. Zhang, “Plasmon Resonance of Finite One-Dimensional Au Nanoparticle Chains,” Nano Lett. **4**, 1067–1071 (2004). [CrossRef]

**10. **C. L. Haynes, A. D. McFarland, L. L. Zhao, R. P. Van Duyne, G. C. Schatz, L. Gunnarsson, J. Prikulis, B. Kasemo, and M. J. Kall, “Nanoparticle Optics: The Importance of Radiative Dipole Coupling in Two-Dimensional Nanoparticle Arrays,” Phys. Chem. B **107**, 7337–7342 (2003). [CrossRef]

**11. **S. Zou, N. Janel, and G. C. Schatz, “Silver nanoparticle array structures that produce remarkably narrow plasmon lineshapes,” J. Chem. Phys. **120**, 10871–10875 (2004). [CrossRef] [PubMed]

**12. **E. M. Hicks, S. Zou, G. C. Schatz, K. G. Spears, R. P. V. Duyne, L. Gunnarsson, T. Rindzevicius, B. Kasemo, and M. Kall, “Controlling Plasmon Line Shapes through Diffractive Coupling in Linear Arrays of Cylindrical Nanoparticles Fabricated by Electron Beam Lithography,” Nano Lett. **5**, 1065–1070 (2005). [CrossRef] [PubMed]

**13. **F. J. García de Abajo, “Colloquium: Light scattering by particle and hole arrays,” Rev. Mod. Phys. **79**, 1267–1290 (2007). [CrossRef]

**14. **W. H. Weber and G. W. Ford, “Propagation of optical excitations by dipolar interactions in metal nanoparticle chains,” Phys. Rev. B **70**, 125429 (2004). [CrossRef]

**15. **S. Y. Park and D. Stroud, “Surface-plasmon dispersion relations in chains of metallic nanoparticles: An exact quasistatic calculation,” Phys. Rev. B **69**, 125418 (2004). [CrossRef]

**16. **A. F. Koenderink and A. Polman, “Complex response and polariton-like dispersion splitting in periodic metal nanoparticle chains,” Phys. Rev. B **74**, 033402 (2006). [CrossRef]

**17. **A. F. Koenderink, R. de Waele, J. C. Prangsma, and A. Polman, “Experimental evidence for large dynamic effects on the plasmon dispersion of subwavelength metal nanoparticle waveguides,” Phys. Rev. B **76**, 201403(R) (2007). [CrossRef]

**18. **T. Yang and K. B. Crozier, “Dispersion and Extinction of Surface Plasmons in Gold Nanoparticle Chains: Influence of the Air/Glass Interface,” Opt. Express **16**, 8570–8580. [PubMed]

**19. **K. B. Crozier, E. Togan, E. Simsek, and T. Yang, “Experimental measurement of the dispersion relations of the surface plasmon modes of metal nanoparticle chains,” Opt. Express **15**, 17482–17493 (2007). [CrossRef] [PubMed]

**20. **Think of a plane wave propagating in z direction. It has an infinite dω/dk_{x} value.