Point-dipole approximation for small systems of strongly coupled radiating nanorods

Systems of closely-spaced resonators can be strongly coupled by interactions mediated by scattered electromagnetic fields. In large systems the resulting response has been shown to be more sensitive to these collective interactions than to the detailed structure of individual resonators. Attempts to describe such systems have resulted in point-dipole approximations to resonators that are computationally efficient for large resonator ensembles. Here we provide a detailed study for the validity of point dipole approximations in small systems of strongly coupled plasmonic nanorods, including the cases of both super-radiant and subradiant excitations, where the characteristics of the excitation depends on the spatial separation between the nanorods. We show that over an appreciable range of rod lengths centered on 210 nm, when the relative separation kl in terms of the resonance wave number of light k satisfies \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$kl\gtrsim \pi /2$$\end{document}kl≳π/2, the point electric dipole model becomes accurate. However, when the resonators are closer, the finite-size and geometry of the resonators modifies the excitation modes, in particular the cooperative mode line shifts of the point dipole approximation begin to rapidly diverge at small separations. We also construct simplified effective models by describing a pair of nanorods as a single effective metamolecule.

Plasmonic nanorods are the most basic form of optical resonators. The scattering of light from any resonator has the ability to produce strong interactions that can result from the wave repeatedly scattering off the same resonator. The EM coupling between different resonators results in different eigenmodes of response, the strong interaction limits of which can be most easily achieved in microwave systems with low ohmic losses, but can also manifest in plasmonic systems 1,2 . The eigenmodes can also destructively interfere and manifest as Fano resonances in the transmitted field [3][4][5] , whose narrow resonances potentially make them useful in applications such as plasmonic rulers 6 and biosensors 7 . Designing material structures to support Fano resonances is difficult; not least due to the complex interactions of different modes, but variations in the resonators can affect the line shifts and widths of the interacting modes also 8 .
Apart from applications in plasmonics and nanophotonics, optical resonators (and nanorods in particular) that are much smaller than the wavelength of the driving light are now commonly used as the building blocks of metamaterials -artificial material composites that are designed to interact with light in ways no conventional materials can. Functionalities of metamaterials include perfect absorption 9 and optical magnetism 10 ; with potential applications ranging from cloaking [11][12][13] to perfect lenses [14][15][16] . An important part of understanding how metamaterials realize their functions is knowledge of the electromagnetic (EM) fields scattered by the constituent resonators.
Metamaterials that exhibit strong collective interactions are becoming increasingly popular with experimentalist 2,17-23 . However, modeling the EM interactions in large resonator systems is generally challenging. The interactions among the resonators can be simplified, e.g., by treating the array as an infinite lattice 24,25 or the resonators as point multipole sources 2,[26][27][28] . Point multipole descriptions in particular have been successful in modeling the cooperative response in planar arrays, e.g., developing electron-beam-driven light sources 29 and transmission properties 30,31 . Point dipole descriptions can also be extended to more complex metamolecules, such as those exhibiting toroidal dipoles 32 .
Here, we analyze the accuracy of the dipole approximation in more detail in small systems of plasmonic nanorods and show, qualitatively, at what separations the finite-size of a nanorod, and its near fields, must be accounted for. Our theoretical model does not require solving the full Maxwell's equations for a resonator ensemble; which is computationally demanding for more than a few resonators 33 . Rather, we utilize the formalism developed in ref. 27 to produce a system of coupled equations for the dynamics of the EM interactions of the scattered and incident EM fields. The method relies on capturing the fundamental physics of each resonator, e.g., its decay rate and resonance frequency, relevant for the radiative coupling between resonators. Our model is readily applied to more complex systems such as those whose resonators are distributed over two planes (one above another) with non-uniform orientations, e.g., a toroidal metamolecule 32 .
Finally, in our paper we also provide an alternative approach for treating each element of a nanorod as a separate meta-atom when we model closely spaced nanorods as a single effective metamolecule. This can notably reduce the number of degrees of freedom in the system.

Method
In analyzing the EM interactions between plasmonic nanorods and an incident EM field, we utilize the general theory derived in detail in ref. 27 , specifically, for the point electric dipole approximation. We regard nanorods as cylinder-shaped resonators and study their longitudinal polarization excitation; where the charge and current oscillations are assumed to be linear along the axis of the nanorod. Here, we provide a brief overview of our model, a more detailed description is provided in the Supplementary Material.
An incident electric displacement field D in (r,t) = D in ê in exp(ik in ⋅r − iΩ 0 t) and magnetic induction µ ε in 0 0 in i n , with frequency Ω 0 , polarization vector ê in and propagation vector k in drive each resonator j's internal charge and current sources. These source oscillations behave in a manner comparable to an LC circuit with resonance frequency ω j Here, L j and C j are respectively, an effective self-inductance and self-capacitance. A dynamic variable with units of charge Q j (t) and its time derivative, the current = I t Q t ( ) ( ) j j , describes the state of current oscillations within each resonator j. In order to analyze the coupled equations for the EM fields, we introduce the slowly varying normal mode oscillator amplitudes b j (t) 27 , with generalized coordinate the charge Q j (t) and conjugate momentum φ j (t), where In the rotating wave approximation the conjugate momentum and current are linearly-proportional 27 . We use Eq. (2) to describe a general resonator with sources of both polarization P j (r, t) and magnetization M j (r,t). The resonator's scattered EM fields result from the oscillations of Q j (t) and I j (t); which are proportional to P j (r, t) and M j (r, t), respectively 27 .
The collective interactions of N resonators with each other and an external field is described by the linear system of equations 27 Here, ḃ is the rate of change of b; a vector of N normal oscillator variables, and F in is a vector describing the interaction of resonator j with the incident EM field. The N×N interaction matrix  requires solving the scattered EM fields for the polarization and magnetization sources. The EM interactions between different resonators i ≠ j are described by the off-diagonal elements of . The strength of ≠ i j  depends on the separation and orientation of the resonators, and naturally accounts for the reduced coupling between elements at the edges compared to the center of a lattice, see Supplementary Material. The diagonal elements describe the EM interactions of a resonator with itself, resulting in the resonator's decay rate Γ j and resonance frequency; In our model we consider magnetization due to induced macroscopic currents. In a straight rod even though the induced current is non-zero, it is linear and therefore the corresponding magnetization is negligible; M j (r, t) ≃ 0. Thus, the scattered EM fields result from a nanorod's polarization sources P j (r, t) alone. This results in an effective accumulation of charge on the nanorod's ends. Analogous simulation methods can also be used to study collective responses of arrays of other resonant emitters, such as atoms 34 . Finite-size resonator model. For our finite-size model we consider a nanorod with a radius a and height H j , see Fig. 1. The polarization density is a uniform distribution of atomic point electric dipoles (with orientation vectors d j ) throughout the volume of the nanorod. For a single nanorod centered at the origin and aligned along the z axis, i.e., =d z, the spatial profile distribution of the polarization density is www.nature.com/scientificreports www.nature.com/scientificreports/ where Θ is the Heaviside function and ρ < a. An individual nanorod experiences radiative decay Γ E1,j , resulting from the interaction of the scattered field with the rod itself. The total decay rate of the rod Γ j also includes a phenomenological ohmic loss rate Γ O,j , where The orientation of the electric dipole is described by the unit vector d j and H j has dimensions of length. The rate at which the electric dipole radiates as a result of its self interactions is Γ E1,j . The radiative decay, ohmic losses and total decay rate are given in Equation (6).
A finite-size effective metamolecule. Thus far, we have only considered the interactions between individual nanorods and point electric dipole resonators. When there are a large number of resonators, one may look to optimize the model. Because resonators are often arranged in a lattice framework, it is possible to consider closely spaced parallel pairs of nanorods as a single effective metamolecule. Here, we extend our finite-size nanorod model to include this effective metamolecule with the same properties as its constituent nanorods, i.e., electric dipole properties.
We consider two parallel nanorods with location vectors , and polarization densities P ±,j (r, t), respectively. When l is small, we may approximate the pair as a single metamolecule with location vector r j = [x j , y j , z j ]. The metamolecule may be symmetrically excited, i.e., = or it may be antisymmetrically excited, i.e., P +,j (r, t) = −P −,j (r, t), where P +,j (r, t) is defined in Equation (5).
The resonance frequency and radiative decay rate of our effective metamolecule depend on the interactions of the individual nanorods. We calculate these properties later; by analyzing in detail a pair of parallel nanorods and point electric dipoles. In principle, if the nanorods are approximated as point electric dipoles, one may treat a closely spaced parallel pair as a single emitter. If both the electric dipoles are symmetrically excited, there is an effective single point electric dipole 35 . However, if both the dipoles are antisymmetrically excited, there is an effective emitter with both an electric quadrupole and a magnetic dipole moment 35 . In this work, however, we only consider the finite-size effective resonators.

Results
In this section, we analyze the EM interactions of nanorods, both as point electric dipole emitters and accounting for their finite-size and geometry. We formulate the model for the nanorods by assuming all the nanorods are made of gold and have equal length, i.e., H j = H 0 for all j. As the nanorods are identical, they experience identical ohmic losses, radiative, and total decay rates, i.e., Γ O,j = Γ O , Γ E1,j = Γ E1 , and Γ j = Γ. We choose λ = .
  37 for the resonant scattering of light from plasmonic nanoparticles, and ohmic losses are accounted for in the Drude model. To simplify the comparison between our point dipole approximation and finite-size nanorod model, we use the same resonance frequency, radiative decay, ohmic losses and total decay rate in both models.
We analyze the characteristic response of the system in the absence of an incident EM field, studying the characteristic collective modes represented by the eigenvectors v n of the interaction matrix . The corresponding complex eigenvalues ξ n , describe the collective mode's characteristic linewidth (real part) and resonance www.nature.com/scientificreports www.nature.com/scientificreports/ frequency shift (imaginary part); ξ n = −γ n /2 − i(Ω n − ω 0 ). The number of modes is determined by the number of resonators, and the radiation may be suppressed γ n < Γ (subradiant), or enhanced γ n > Γ (superradiant). two parallel nanorods. As our first example, we consider two parallel nanorods (and two point electric dipoles) with location vectors r 1 = −r 2 = [0, l/2, 0]. The coupling matrix  has two eigenmodes of current oscillation. In the first mode, the current oscillations are in-phase (symmetric) with d d 1 2   = . In the second mode, the current oscillations are out-of-phase (antisymmetric) with = −d d 1 2 . The former we denote by a subscript 's'; and latter by a subscript 'a' .
In Fig. 2, we show the collective eigenmode's radiative resonance line shifts and linewidths for two finite-size parallel nanorods and two parallel point electric dipoles. Throughout the range of kl, the decay rates γ n of the point dipole model closely agree with the corresponding decay rates of the finite-size nanorod model. When kl ≈ π, the symmetric mode is subradiant with γ s ≈ 0.9Γ, and the antisymmetric mode is superradiant with γ a ≈ 1.1Γ; for both the nanorods and point electric dipoles. As the separation reduces, the symmetric mode becomes superradiant and the antisymmetric mode subradiant. When kl ≈ 2π/3, the decay rates of the point multipole approximation are: γ a ≈ 0.7Γ (subradiant); and γ s ≈ 1.3Γ (superradiant); the finite-size model shows decay rates: γ a ≈ 0.8Γ; and γ s ≈ 1.2Γ. As the separation becomes small kl ≈ π/6, the antisymmetric linewidths approach the ohmic loss rate, γ a ≈ 0.2Γ and the symmetric mode linewidths become more superradiant γ s ≈ 1.8Γ.
The lineshifts δω n , however, only agree when  π kl /2. As the separation becomes small  π kl /2, the line shift of the point electric dipole model begins to separate from the corresponding line shift of the finite-size resonator model. The line shift of the finite-size resonator model here is . For  π kl /2, the antisymmetric mode line shift of the point electric dipole model is blue shifted from Ω 0 , and the symmetric mode red shifted.
In Fig. 2, the nanorods have length λ = . , where λ  570 nm t . For both longer and shorter rod systems, the point dipole approximation becomes valid when the separation is  π kl /2, where k is the resonance wavenumber of the light. When we make small changes to the nanorod radius λ p /6 < a < λ p /4; we also find the point dipole approximation is valid for separations  π kl /2.

two interacting pairs of nanorods.
In this section, we analyze two parallel pairs of horizontal nanorods.
Firstly, we treat the interacting pairs as four discrete finite-size nanorods that have a non-vanishing polarization density. Secondly, we optimize the model by treating the two pairs as two effective metamolecules. In each case, we compare the model to four discrete point electric dipoles.
Four discrete nanorods. In general, we position the jth pair of nanorods at r ±,j = [x j , y j ± l/2, z j ]. In our example, we set y j = z j = 0, kx j = π for all j, and vary the parameter l. When the interactions between individual nanorods are considered, there are four collective modes, see Fig. 3. When each parallel pair of nanorods are symmetrically excited, they can be approximated as out-of-phase (E1a) and in-phase (E1s) effective electric dipoles, see Fig. 3(a,c), respectively. When the nanorods in each pair are antisymmetrically excited, they can be likened to two resonators with both electric quadrupole and magnetic dipole moments, where each pair are out-of-phase (E2a) or in-phase (E2s), see Fig. 3(b,d), respectively. In Fig. 4, we show the collective mode resonance line shifts and linewidths of four point electric dipole resonators and those of four interacting finite-size nanorods. Again, we find the collective mode decay rates of the  . At π  kl 9 /10, the E1s mode becomes noticeably subradiant with γ . Γ  0 9 E1s . Here also, the E2a mode becomes noticeably superradiant with γ . Γ  1 1 E2a . The E1a mode remains superradiant while the E2s mode remains subradiant throughout the range.
Two effective metamolecules. When we approximate each pair of nanorods as an effective single resonator; if the nanorods' current oscillations are in-phase there is an effective electric dipole resonator, if the current oscillations are out-of-phase there is an effective resonator with both electric quadrupole and magnetic dipole responses 35 . In  www.nature.com/scientificreports www.nature.com/scientificreports/ principle, cross coupling can occur between the effective resonators whereby an in-phase pair and an out-of-phase pair get mixed due to interactions. However, here we consider only two interacting in-phase pairs, and separately, two out-of-phase pairs; neglecting such processes. There are two modes of collective oscillation for each effective resonator system; symmetric and antisymmetric. We also designate these as E1s and E1a, respectively, when the nanorods within each pair are in-phase, and E2s and E2a, respectively, when the nanorods within each pair are out-of-phase.
In Fig. 2, we calculated the collective mode decay rates for two parallel nanorods. These decay rates, γ s and γ a , provide us the total decay rate for in-phase and out-of-phase pairs of nanorods, respectively. Also in Fig. 2, we calculated the line shifts of the collective modes, this allows us to determine the resonance frequencies for in-phase and out-of-phase pairs of nanorods; Ω s , and Ω a , respectively. In general, Ω s , Ω a ≠ ω 0 , this means that in our effective resonator model, the diagonal elements of  also contain the imaginary component; s,a , where δω s,a are the line shifts of two in-phase and out-of-phase parallel nanorods, respectively, see Fig. 2.
In Fig. 5 , δω ] E2s (2a) of the collective modes of oscillation of the N = 2 effective interacting pairs of nanorods (denoted by the superscripts (2s) for in-phase and (2a) for out-of-phase pairs) vary with the parameter l, and compare to the corresponding modes' line widths and shifts of N = 4 point electric dipole system (denoted by the superscript (1)).
The linewidths resulting from the interacting N = 4 point electric dipoles closely agree with those of the N = 2 effective metamolecule models, both when in-phase and when out-of-phase. Also, comparing Fig. 5(a,b), with Fig. 4, the effective metamolecules' linewidths closely match those of the finite-size nanorod model.
The point electric dipole model's line shifts in Fig. 5(c,d) begin to separate from the corresponding effective metamolecule line shifts when π  kl /2. As kl reduces, the point electric dipole model's line shifts for the E2a and E2s (E1a and E1s) modes blue shift (red shift) from the out-of-phase (in-phase) effective metamolecule's corresponding modes. Comparing Fig. 5(c,d) with the corresponding line shifts in Fig. 4, the line shifts of the finite-size nanorod model closely agree with those of the effective metamolecule model.

Conclusions
Understanding the complex EM interactions in resonator ensembles is important for the design of metamaterials. Our circuit element resonator model provides an efficient way of understanding the dynamics of the system without having to fully solve Maxwell's equations. We have studied the strong collective modes of current oscillations resulting from the EM interactions in closely spaced resonator systems. These collective modes have an associated radiative response that can be either superradiant or subradiant, and together with the resonance line shift, is strongly influenced by the spatial separation of the resonators. Though in this work, we have considered all the resonance frequencies of the resonators to be equal, variation in the resonances (inhomogeneous broadening) can generally suppress the collective radiation interactions 38 .
We have analyzed, in detail, the validity of the point electric dipole approximation of interacting resonators in small systems; demonstrating how we can model plasmonic nanorod systems both as point electric dipole resonators and accounting for their finite-size and geometry. In particular, we have determined how interacting discrete nanorods with an appreciable range of lengths centered on λ .   H 0 24 210 nm 0 0 can be approximated as www.nature.com/scientificreports www.nature.com/scientificreports/ interacting point electric dipoles, especially when their separation is greater than  π kl /2. For closely spaced resonators  π kl /2, their finite-size and geometry becomes increasingly important. An alternative approach for treating each resonator as a separate meta-atom is to model closely spaced resonators as a single effective metamolecule, reducing the number of degrees of freedom. In principle, this could be extended to other more complex effective metamolecules, e.g, toroidal metamolecules 32 .