Plasmonic nanoclusters: a path towards negative-index metafluids

We introduce the concept of metafluids— liquid metamaterials based on clusters of metallic nanoparticles which we will te rm Artificial Plasmonic Molecules(APMs). APMs comprising four nanoparticles in a tetrahedral arrangement have isotropic electric and magne tic responses and are analyzed using the plasmon hybridization (PH) method, a n electrostatic eigenvalue equation, and vectorial finite element frequenc y domain (FEFD) electromagnetic simulations. With the aid of group theory, we identify the resonances that provide the strongest electric and magneti c response and study them as a function of separation between spherical nan oparticles. It is demonstrated that a colloidal solution of plasmonic tetrah edral nanoclusters can act as an optical medium with very large, small, or even ne gative effective permittivity,εeff, and substantial effective magnetic susceptibility, χeff = μeff − 1, in the visible or near infrared bands. We suggest paths for increasing the magnetic response, decreasing the damping, and developing a metafluid with simultaneously negative εeff andμeff. © 2007 Optical Society of America OCIS codes: (160.1245) Materials: Artificially engineered materials; (240.6680) Optics at surfaces: Surface plasmons; (260.2065) Physical optics: E ffective medium theory; (350.3618) Other areas of optics: Left-handed materials. References and links 1. H. Metiu, “Surface enhanced spectroscopy,” Prog. Surf. S ci. 17, 153 (1984). 2. G. C. Schatz and R. P. van Duyne, “Electromagnetic mechani sm of surface-enhanced spectroscopy,” in Handbook of Vibrational Spectroscopy , J. M. Chalmers and P. R. Griffiths, eds., pp. 1–16 (John Wiley , Chichester, 2002). 3. M. Moskovits, L. Tay, J. Yang, and T. Haslett, “SERS and the single molecule,” Top. Appl. Phys. 82, 215 (2002). 4. J. B. Jackson and N. J. Halas, “Surface Enhanced Raman Scat tering on tunable plasmonic nanoparticle substrates,” Proc. Nat. Acad. Sci. USA 101, 17,930 (2004). 5. D. R. Ward, N. K. Grady, C. S. Levin, N. J. Halas, Y. Wu, P. Nor dlander, and D. Natelson, “Electromigrated nanoscale gaps for surface-enhanced Raman spectroscopy,” Nano Lett.7, 1396–1400 (2007). 6. D. J. Anderson and M. Moskovits, “A SERS-active system bas ed on silver nanoparticles tethered to a deposited silver film,” J. Phys. Chem. B110, 13,722 (2006). 7. P. K. Jain, S. Eustis, and M. A. El-Sayed, “Plasmon couplin g in nanorod assemblies: Optical absorption, discrete dipole simulation, and exciton coupling model,” J. Phys. Ch em. B110, 18,243 (2006). 8. V. Manoharan, M. Elsesser, and D. Pine, “Dense Packing and Symmetry in Small Clusters of Microspheres,” Science301, 483 (2003). 9. G.-R. Yi, V. N. Manoharan, E. Michel, M. T. Elsesser, S.-M. Yang, and D. J. Pine, “Colloidal Clusters of Silica or Polymer Microspheres,” Adv. Mater. 16, 1204 (2004). 10. V. N. Manoharan and D. J. Pine, “Building Materials by Pac king Spheres,” Mater. Res. Soc. Bull. 29, 91 (2004). #85978 $15.00 USD Received 2 Aug 2007; revised 7 Sep 2007; accepted 8 Sep 2007; published 11 Oct 2007 (C) 2007 OSA 17 October 2007 / Vol. 15, No. 21 / OPTICS EXPRESS 14129 11. Y.-S. Cho, G.-R. Yi, S.-H. Kim, D. J. Pine, and S.-M. Yang, “Colloidal Clusters of Microspheres from Water-inOil Emulsions,” Chem. Mater. 17, 5006 (2005). 12. D. Psaltis, S. R. Quake, and C. Yang, “Developing optoflui dic technology through the fusion of microfluidics and optics,” Nature442, 381 (2006). 13. M. I. Stockman, K. Li, S. Brasselet, and J. Zyss, “Octupol ar metal nanoparticles as optically driven coherently controlled nanorotors,” Chem. Phys. Lett. 433, 130–135 (2006). 14. W. van Megen, T. C. Mortensen, S. R. Williams, and J. Müll er, “Measurement of the self-intermediate scattering function of suspensions of hard spherical particles near th e glass transition,” Phys. Rev. E 58, 6073 (1998). 15. O. N. Singh and A. Lakhtakia, eds., Electromagnetic Waves in Unconventional Materials and Str uctures(John Wiley, New York, 2000). 16. E. Prodan, C. Radloff, N. J. Halas, and P. Nordlander, “A h ybridization model for the plasmon response of complex nanoparticles,” Science 302, 419 (2003). 17. D. J. Bergman and D. Stroud, “Theory of resonances in the e lectromagnetic scattering by macroscopic bodies,” Phys. Rev. B22, 3527 (1980). 18. I. D. Mayergoyz, D. R. Fredkin, and Z. Zhang, “Electrosta tic (plasmon) resonances in nanoparticles,” Phys. Rev. B 72, 155,412 (2005). 19. E. Prodan and P. Nordlander, “Plasmon hybridization in s pherical nanoparticles,” J. Chem. Phys. 120, 5444 (2004). 20. H. Wang, D. W. Brandl, P. Nordlander, and N. J. Halas, “Pla smonic Nanostructures: Artificial Molecules,” Acc. Chem. Res. 40, 53–62 (2007). 21. E. Hao and G. C. Schatz, “Electromagnetic fields around si lver nanoparticles and dimers,” J. Chem. Phys. 120, 357 (2004). 22. Y. Urzhumov and G. Shvets, “Applications of Nanoparticl e Arrays to Coherent Anti-Stokes Raman Spectroscopy of Chiral Molecules,” Proc. SPIE5927, 1D–1 (2005). 23. D. Korobkin, Y. Urzhumov, B. Neuner III, C. Zorman, Z. Zha ng, I. D. Mayergoyz, and G. Shvets, “Mid-infrared metamaterial based on perforated SiC membrane: Engineerin g optical response using surface phonon polaritons,” Appl. Phys. A88, 605–609 (2007). 24. Y. Urzhumov, D. Korobkin, B. Neuner III, C. Zorman, and G. Shvets, “Optical Properties of Sub-Wavelength Hole Arrays in SiC Membranes,” J. Opt. A: Pure Appl. 9, S1–S12 (2007). 25. D. W. Brandl, N. A. Mirin, and P. Nordlander, “Plasmon mod es of nanosphere trimers and quadrumers,” J. Phys. Chem. B110, 12,302 (2006). 26. G. F. Koster, J. O. Dimmock, R. G. Wheeler, and H. Statz, Properties of the Thirty-Two Point Groups (MIT Press, Cambridge, Mass., 1963). 27. A. K. Sarychev, G. Shvets, and V. M. Shalaev, “Magnetic pl asmon resonance,” Phys. Rev. E 73, 036,609 (2006). 28. G. Shvets and Y. Urzhumov, “Negative index meta-materia ls b sed on two-dimensional metallic structures,” J. Opt. A: Pure Appl. Opt.8, S122–S130 (2006). 29. V. Lomakin, Y. Fainman, Y. Urzhumov, and G. Shvets, “Doub ly negative metamaterials in the near infrared and visible regimes based on thin film nanocomposites,” Opt. Exp ress14(23), 11,164 (2006). 30. G. Shvets and Y. Urzhumov, “Engineering the Electromagn etic Properties of Periodic Nanostructures Using Electrostatic Resonances,” Phys. Rev. Lett. 93, 243,902 (2004). 31. G. Shvets and Y. Urzhumov, “Electric and magnetic proper ties of sub-wavelength plasmonic crystals,” J. Opt. A: Pure Appl. Opt.7, S23–S31 (2005). 32. E. D. Palik, ed., Handbook of Optical Constants of Solids (Academic Press, Orlando, 1985). 33. I. D. Mayergoyz and Z. Zhang, “The Computation of Extinct ion Cross-sections of Resonant Metallic Nanoparticles Subject to Optical Radiation,” IEEE Trans. Magn. PF4-1, CEFC10,037 (2006). 34. Y. Urzhumov and G. Shvets, “Quasistatic effective mediu m theory of plasmonic nanostructures,” to appear in Proc. SPIE (2007). 35. P. Markos and C. Soukoulis, “Transmission properties an d effective electromagnetic parameters of double negative metamaterials,” Opt. Express 11, 649 (2003). 36. D. R. Smith, D. C. Vier, T. Koschny, and C. M. Soukoulis, “E lectromagnetic parameter retrieval from inhomogeneous metamaterials,” Phys. Rev. E 71, 036,617 (2005). 37. T. Koschny, P. Markos, D. R. Smith, and C. M. Soukoulis, “R esonant and antiresonant frequency dependence of the effective parameters of metamaterials,” Phys. Rev. E 68, 065,602 (2003). 38. A. L. Efros, “Comment II on “Resonant and antiresonant fr equency dependence of the effective parameters of metamaterials”,” Phys. Rev. E 70, 048,602 (2004). 39. C. G. de Kruif, E. M. F. van Iersel, A. Vrij, and W. B. Russel , “Hard sphere colloidal dispersions: Viscosity as a function of shear rate and volume fraction,” J. Chem. Phys. 83, 4717 (1985). 40. D. Kang and N. A. Clark, “Fast Growth of Silica Colloidal C rystals,” J. Kor. Phys. Soc. 41, 817 (2002). 41. P. N. Pusey and W. van Megen, “Phase behaviour of concentr at d suspensions of nearly hard colloidal spheres,” Nature320, 340 (1986). 42. A. Kassiba, M. Makowska-Janusik, J. Boucle, J. F. Bardea u, A. Bulou, N. Herlin, M. Mayne, and X. Armand, “Stoichiometry and interface effects on the electronic and optical properties of SiC nanoparticles,” Diamond and #85978 $15.00 USD Received 2 Aug 2007; revised 7 Sep 2007; accepted 8 Sep 2007; published 11 Oct 2007 (C) 2007 OSA 17 October 2007 / Vol. 15, No. 21 / OPTICS EXPRESS 14130 Related Materials11, 1243 (2002). 43. D. Korobkin, Y. Urzhumov, C. Zorman, and G. Shvets, “Far F ield Detection of the Super-Lensing Effect in Mid-Infrared: Theory and Experiment,” J. Mod. Opt. 52(16), 2351 (2005). 44. D. Korobkin, Y. Urzhumov, and G. Shvets, “Enhanced nearfield resolution in mid-infrared using metamaterials,” J. Opt. Soc. Am. B23, 468 (2006). 45. T. Taubner, D. Korobkin, Y. Urzhumov, G. Shvets, and R. Hi llenbrand, “Near-Field Microscopy Through a SiC Superlens,” Science 313, 1595 (2006).


Introduction
The optical properties of metallic multi-nanoparticle structures have been of great theoretical and experimental interest in recent years due to biological and chemical sensing applications, including Surface Enhanced Raman Spectroscopy (SERS) and Localized Surface Plasmon (LSPR) sensing [1,2,3,4,5,6,7]. In the former, large electric field enhancements near the surfaces of particles or in the gaps of nanoparticle clusters near the plasmon frequencies lead to an increased Raman cross section. In the latter, a change of refractive index from a nearby molecule causes a red-shift of the plasmon frequencies. Plasmonic nanostructures have also attracted a great deal of attention as an approach to construct electromagnetic metamaterialsmedia with optical properties previously unavailable in nature.
Here, we introduce a new concept called a metafluid -a liquid metamaterial containing Artificial Plasmonic Molecules (APMs). APMs are geometrically ordered aggregates of plasmonic nanoparticles that typically consist of 2-15 individual "atoms" [8,9,10,11]. With the flexibility to engineer APM geometries, the optical properties of APMs can differ tremendously from those found in natural molecules. The size of an APM greatly exceeds that of a typical molecule yet may be considerably smaller than the optical wavelength. Due to the small spatial extent of the APMs compared to optical wavelengths, the resulting metafluid can still be viewed as an effective medium and characterized by its effective coefficients such as, for example, dielectric permittivity and magnetic permeability. By changing the size and arrangement of the constituent plasmonic nanoparticles inside an APM, the APM's optical response at the frequency of interest can be controlled in both magnitude (strong or weak) and character (electric or magnetic, scattering or dissipative). Recent interest in liquid-liquid optical waveguides [12] further motivates the development of metafluids.
The term "metafluid" in this paper is composed of two words: metamaterial and fluid. By "metamaterial" we mean an artificially created composite of regular materials that exhibits unusual electromagnetic properties, such as, for example, negative magnetic permeability or negative index of refraction. One can ascribe effective index of refraction to a composite medium, for example, when the structure is periodic (regardless of the distances between particles and their sizes). Such situation is known as the Bloch-Floquet regime or the photonic crystal regime. Assignment of refractive index is also possible in the effective medium regime, i. e. when the size of individual scatterers is much smaller than the wavelength in immersion medium. These two regimes are not mutually exclusive when the distance between particles is sub-wavelength; in Section 4 we take advantage of periodic boundary conditions to characterize the optical properties of a dense nanoparticle colloid. Since optical parameters of an effective medium depend mostly on the average distance between identical particles in the ensemble, and little on the locations of individual particles, period-independent spectral features of periodic ensembles must be shared by all random ensembles with the same particle number density.
In this paper we theoretically investigate an APM composed of four metallic nanospheres situated equidistant from one another at the vertices of a regular tetrahedron, as the first candidate for optical metafluids. We refer to this structure as the tetramer. This structure has recently attracted attention as a candidate for a coherently controlled nanorotor [13]. For electromagnetic metafluids, APMs with tetrahedral symmetry are attractive because their single-particle polarizabilities are orientation-independent. The effective dielectric tensor,ε, of most fluids is effectively a scalar because of the rapid rotation and high spatial density (∼ 10 23 cm −3 ) of the constituent molecules. However, when gigantic artificial plasmonic molecules described in this paper are part of the metafluid, their rotational frequency and concentration in the solution may not be sufficient to provide isotropization by temporal and spatial averaging. Therefore, the isotropic polarizability of tetrahedral plasmonic molecules becomes crucial for ensuring that the tensorsε andμ of a metafluid are spherical. Isotropy of particles also helps reduce the effects of diffusive (Brownian) motion on the optical properties of a liquid. Dynamic light scattering effects associated with fluctuations of particle positions [14] are not considered in this paper.
In effective medium regime, appropriate quantities to describe the propagation of a plane wave are dielectric and diamagnetic susceptibilities of the compound medium. In general, for a linear medium there are four such (tensor) quantities. In the so-called Tellegen representation, they relate electric and magnetic induction with electric and magnetic field intensity [15]: In the most general case where magneto-electric susceptibilitiesξ ,ζ do not vanish, the medium is called bianisotropic [15]. While a bianisotropic medium may be interesting on its own, in this article we focus only on regular media, which can be described with an effective permittivityε and permeabilityμ only. In Section 2.1, we prove that the tetrahedral group T d of the tetramer has sufficient symmetry to prohibit bianisotropy in the electromagnetic response of sub-wavelength plasmonic particles. Thus, an effective medium composed of tetramers is isotropic, non-chiral, and described by two scalar quantities, ε eff and µ eff . It should be mentioned that the tetramer is not the only metamolecule that forms metafluids with such properties. There are seven 3-dimensional point groups which guarantee a second-rank tensor to be spherical: three chiral groups (T , O, I) and 4 non-chiral groups (T d , O h , T h , I h ). We note that, in general, magneto-electric coupling termsξ andζ do not average to zero when accounting for the rotation of APMs. This means, for example, that a medium consisting of chiral APMs is also chiral. The remaining 4 groups can be utilized for the design of isotropic, non-bianisotropic optical metamaterials. The minimum number of identical spherical nanoparticles is 4 for T d symmetry (vertices of a regular tetrahedron), 6 for O h (octahedron), 12 for I h (icosahedron) and 20 for T h (pyritohedron), as summarized in Table 2. A tetramer is thus the minimal non-bianisotropic, fully isotropic "metamolecule".
Experimental routes exist to assemble colloidal nanoparticles into highly ordered clusters. One experimental route to the assembly of tetramers and larger symmetric structures, including isotropic 6-particle "octamers" and 12-particle "icosamers", is particle clustering in an oilin-water emulsion process [8]. Particles are first functionalized to be hydrophobic and then transferred to an oil solvent. The oil is then added to water with surfactant and sheared in a homogenizer, which yields surfactant-stabilized oil droplets in water. Next, the oil is evaporated from the emulsion, and particles are forced into clusters due to capillary forces and are held solidly together by van der Waals forces. Tetrahedral clusters are separated from clusters of other particle number by centrifugation in a density gradient. This technique for creating clusters is versatile and applies to all types of particles ranging from silica to PMMA [9]. In addition, clustering is possible for hydrophilic particles through a water-in-oil emulsion [11].
The remainder of the paper is organized as follows: in Section 2, we use the plasmon hybridization method [16] and a new finite-element implementation of the Electrostatic Eigenvalue surface integral formalism [17, 18] to find the plasmon modes of a tetramer. Using group theory, we classify the plasmon modes by their electric and magnetic properties. In Section 3, we examine the optical absorption spectra in the tetramer system using finite element frequency domain (FEFD) calculations. In Section 4, we characterize a tetramer colloid as an effective medium with isotropic dielectric permittivity and magnetic permeability, and discuss some experimental challenges concerned with the fabrication of negative index metafluids.

Quasi-static analysis of the plasmon modes of tetramers
The Plasmon Hybridization (PH) method [16], along with several other methods such as the Electrostatic Eigenvalue (EE) approach [17,18], provides exact solutions for the plasmon resonances of a complex nanostructure in the quasi-static (non-retarded) limit.
In the PH approach, the plasmon modes of a multi-nanoparticle systems are expressed as a linear combination of the (primitive) plasmon modes of the individual particles [19]. The primitive plasmon modes interact with each other through the Coulomb forces induced by their surface charges. An appealing feature of the PH approach is that its eigenvalue problem is very similar to the eigenvalue problem for molecular orbitals in quantum chemistry. This analogy gives an insight into the relationship between the plasmon modes of a composite structure and the plasmons of its constituent particles, and encourages the use of the group theory for symmetry classification of these modes [20]. In contrast with the coupled-dipole approach [21,22], the PH method can account for the hybridization of primitive eigenmodes with arbitrary multipole order l.
The Electrostatic Eigenvalue approach [17, 18] is a general method for computing electrostatic eigenfunctions of arbitrarily-shaped particles or their ensembles. The EE method does not take advantage of the simplicity of primitive plasmon modes, which makes it applicable also to structures with non-spherical particles [23,24]. Our implementation of EE method is described in Section 2.3.

Group theory analysis
There are two main components to understanding the plasmon modes of nanoparticle clusters as linear combinations of individual particle plasmons: multipolar mixing and geometric symmetry. The former is elaborated upon in Section 2.2. Here, we discuss the application of group theory for the determination of specific linear combinations of individual particle plasmons that make up the tetramer plasmon modes.
The plasmon modes of nanoparticle clusters have previously been shown to map onto the irreducible representations of a point group corresponding to the underlying symmetry of the cluster [25]. The plasmon modes of a planar equilateral trimer and quadrumer were classified according to the irreducible representations that make up the dipole representation of the point groups D 3h (trimer) and D 4h (quadrumer). A symmetry adapted basis for each representation could be generated using projection operators for each point-group, and the linear combinations of individual plasmons that compose a physical trimer or quadrumer plasmon mode could be expressed in these bases. The symmetry adapted linear combinations allowed for the separation of cluster plasmon modes into those polarized either in-plane or out-of-plane and were used to identify the modes that could be excited by light for each polarization. In this paper, we use this approach to analyze the plasmon modes of the tetramer.
The symmetry of the tetramer corresponds to the point group T d . It can be shown that a product of four dipole representations of the rotation group O(3) (which describes hybridization of four dipoles arranged tetrahedrally) splits into the following irreducible representations of T d group: This expansion is illustrated by Fig. 1. We note that all irreducible representations of T d , except A 2 , can be found amongst linear combinations of the single-particle dipoles (l = 1) arranged at the vertices of tetrahedron. Modes with A 2 symmetry occur only in decompositions of higher order multipoles, e.g. with l = 3. Note that each irreducible representation corresponds to an nfold degenerate plasmon energy where n is the dimension of the representation. Thus, A 1 , a onedimensional representation, corresponds to one non-degenerate plasmon energy, while each instance of T 2 , a three-dimensional representation, corresponds to a separate, triple-degenerate energy.
Visualization (side) Visualization (corner) Irreducible representation In the quasi-static limit, the strongest coupling of light is to the dipole modes. Therefore, we will focus mostly on the linear combinations of dipole modes. Projection operators may be used on the irreducible representations to obtain the dipolar symmetry adapted basis set associated with each representation in Eq. (3). The resulting linear combinations are shown in Fig. 1. Of these basis functions, only T 2 has a non-vanishing net electric dipole moment.
To prove this fact and to provide additional insight into the electromagnetic properties of tetrahedral clusters, we perform a systematic decomposition of all plasmon modes into electric and magnetic multipoles based on their rotational properties and inversion parity. Such decomposition arises from the observation that all 3-dimensional point groups, including the point group T d of tetrahedron, are subgroups of the full group of the sphere O(3) = SO(3) C i . Irreducible representations of the latter group are well-known [26]; they are characterized by the angular momentum J = 0, 1, 2, . . . and inversion parity P = ±1. An electric 2 J -pole is a rank-J tensor with inversion parity (−1) J . For example, the electric dipole is a vector that changes sign upon inversion. A magnetic 2 J -pole is a rank-J tensor of parity (−1) J+1 . For instance, the magnetic dipole does not change sign upon inversion and is a pseudovector. This identifies all irreducible representations of O(3) as either electric or magnetic multipoles. The complete multipolar decomposition of T d group is presented in Table 1. Its electric multipole part is found in a variety of books on point groups, e. g. [26].
The most important conclusion of this analysis is that only T 2 (T 1 ) modes of tetramers can have non-zero electric (respectively, magnetic) dipole moment. Additionally, Table 1 proves that T d symmetry does not allow plasmon modes to have both electric and magnetic dipole moment; therefore, magneto-electric coupling terms ξ , ζ vanish in tetramer colloids. Considering that clusters of other highly symmetric shapes, ranging from octamers to icosamers, have already been synthesized from various dielectric materials [8], we have performed a similar multipolar analysis for all non-chiral cubic groups. Results that are relevant to electromagnetic homogenization of colloids are condensed in Table 2. We emphasize here that even though magnetism of plasmon resonance is a phenomenon of order η ≡ ω 2 R 2 /c 2 [27], it would be incorrect to assume that it is always associated with electric quadrupole resonances. Though the latter is generally true for lower-symmetry particles [28, 29], higher symmetry groups may enforce nullification of electric quadrupole moment of magnetic dipole modes. For example, magnetically active eigenmodes of a tetramer are electric octupoles. Another example of this sort was reported earlier for square-lattice (C 4v ) plasmonic crystals, where optical magnetism was apparently caused by two-dimensional octupole (

Plasmon energies from PH model
In this section, we use plasmon hybridization to calculate the plasmon modes of the tetramer in the quasi-static limit as a function of the distance between spheres. The metals are implicitly described by a Drude form of the dielectric function, ε p (ω) = ε ∞ − ω 2 B /ω 2 , where ε ∞ is the background permittivity of the metal and the bulk plasmon frequency is ω B = 4πn 2 0 /m e . Theoretical calculations are presented here for tetramers composed of gold nanospheres, and multipolar orders of up to l max = 24 are used. Gold is simulated using ω B = 8.95 eV and a frequency independent ε ∞ = 9.5. This parameterization provides an accurate description of the dielectric properties of Au above 500 nm. To simplify the problem a bit, we have assumed that the solvent and dielectric coating of nanospheres have the same permittivity ε s = n 2 s , where n s = 1.4. Organic solvents with such index are readily available; for instance, ethanol has index of refraction n = 1.36 that matches almost perfectly the index of chemically grown silica (n SiO 2 ≈ 1.39) shells. Solvents with higher molecular weight (e.g., iso-octane with n s = 1.39) can provide even better index matching with silica when necessary; in our calculations the index matching assumption is a good approximation but not a condition of applicability. This setup is assumed throughout the rest of the paper. Fig. 2 shows the plasmon modes of a tetramer as a function of the distance between spheres separated according to the irreducible representation of the tetramer point group. Only the lowest multipolar plasmons l = 1, 2, 3 are included in the graph. For large separations, the interaction between the spheres is weak, and the plasmon energies are essentially the 4(2l + 1)-fold degenerate energies of the 2 l -polar plasmons of a single sphere. As the separation decreases, the primary interaction that determines splitting between the modes is the electric dipolar (1/r 3 ) interaction. As the spheres move closer, interactions of the electric dipole and higher order (l ≥ 2) multipoles with other higher order multipoles (1/r l+l ′ +1 ) become more important.

First-principle electrostatic simulations of Artificial Plasmonic Molecules
The surge of recent interest in the optical properties of plasmonic nanoparticles originates from the unique property of negative-permittivity interfaces to support source-free excitations known as surface plasmons. These excitations exist even for particle sizes much smaller than the wavelength of light at which they occur, which suggests that they are electrostatic in nature. Consequently, they can be found as solutions of the electrostatic Laplace equation with no external field or charge: For homogeneous negative permittivity particles (ε p < 0) in a uniform transparent immersion medium (ε s > 0), this equation can be recast as a linear generalized eigenvalue problem [17] in which the electrostatic permittivity of plasmonic particles plays the role of an eigenvalue: where θ (x) equals 1 inside the particle(s) and zero elsewhere, and s = 1/(1 − ε p /ε s ).
If the boundary of plasmonic particles is sufficiently smooth, differential equation (5) can be reduced to a linear integral equation [18] for electrostatic surface charges σ : where G(x, x ′ ) is the electrostatic Green's function, n(x) is the outward normal to the surface of plasmonic particle, and λ = (ε p − ε s )/(ε p + ε s ) is the electrostatic eigenvalue [18]. Numerical discretization of both versions of the Electrostatic Eigenvalue (EE) method is straightforward and was implemented using FEM software package COMSOL Multiphysics. The differential equation (5) method may be preferable for periodic systems [23, 24], where periodicity is easily imposed as boundary conditions for potential φ , whereas in the surface integral approach (6) periodic boundary conditions must be embedded into the Green's function G(x, x ′ ). On the other hand, the volumetric equation (5) has many more degrees of freedom for the same number of mesh elements on the particle surface than the equivalent equation (6). In general, this leads to a large number of unphysical solutions. The surface integral approach does not suffer from this problem, at least for very smooth surfaces such as spheres.  Table 1) except triplets (T 1 , T 2 ) is presented. Gap-to-diameter ratio in the cluster is 1/10. Triplets are shown separately in Fig. 4. Using the surface charge equation (6), we have performed finite element method (FEM) calculations using experimentally relevant parameters. The lowest-lying resonances of each symmetry type are plotted in Fig. 3. Since the electrostatic spectrum is scale-invariant, the only dimensionless structural parameter in the problem is the ratio of a sphere diameter D to the gap h between their surfaces. Another dimensionless parameter, the dielectric contrast defined as the ratio of the dielectric constant of the particles (ε p ) and the dielectric constant of the solvent (ε s ). The dielectric contrast influences the energies of all electrostatic resonances. The vacuum wavelength λ vac is not a parameter, but rather a label, related unambiguously to the dielectric contrast ε p /ε s .
To understand how electrostatic resonances are excited by incident radiation and how they contribute to optical extinction and absorption spectra, one may start with the quasi-static approximation, as suggested by the sub-wavelength nature of electrostatic resonances.
In the quasi-static approximation, the strongest interaction between incident light and particles is the coupling of a nearly uniform electric field with the induced electric dipole moment of the particles. The strength of this interaction is characterized by the normalized dipole moment of an eigenmode [33, 24, 34]: where σ n (x) is the charge eigenfunction of n th resonance, φ n (x) is its potential, ∂ /∂ n is the normal derivative evaluated on the plasmonic side of the surface, and V p = θ dV is the volume of metal in a cluster. In addition to the strong excitation of electric dipole resonances, which remain strong even in the non-retarded limit, weakly inhomogeneous electric and magnetic fields of an incident electromagnetic wave also induce various electric and magnetic multipoles. Though lots of non-dipolar modes are excited by inhomogeneous fields, only some of them carry a magnetic dipole moment. This induced magnetic moment can be calculated in quasi-static approximation from the total current J = J c + ∂ P/∂t ≡ ε−1 4π ∂ E/∂t in plasmonic particles: After simple transformations, the energy-normalized magnetic dipole moment of an eigenmode can be expressed in terms of surface integrals [34]: Expression (9) demonstrates that only eigenmodes that transform as a pseudovector may have a non-vanishing magnetic dipole moment in the lowest order to retardation parameter η. According to Table 1, pseudovectors (such as the rotation operator R α ) transform under the T 1 representation of the group T d .
Resonant magnetic response due to a plasmon eigenmode is a product of two factors: (i) the resonantly enhanced coupling coefficient between the incident electromagnetic field and the field of a plasmon eigenmode, and (ii) the magnetic dipole moment carried by the mode. Though calculation of the first factor is out of the scope of this paper, our estimate of the second factor rules out all non-T 1 modes from analysis of optical magnetism. The electrostatic nature of surface plasmons suggests that magnetically active plasmon modes are predominantly excited by an inhomogeneous electric field as electric octupoles. The magnetic response of tetramers can be seen as a consequence of the fact that they possess both an electric octupole and a magnetic dipole moment.
In the remainder of this article, we pay attention only to the T 2 and T 1 modes of tetramers. Amongst all T 2 (or T 1 ) modes, the ones with the lowest negative resonant permittivity eigenvalue (thus, the lowest frequency) should have the strongest electric (respectively, magnetic) response. This is because higher eigenmodes of a particular symmetry type have additional sign changes in their surface charge eigenfunction, σ n , resulting in a smaller coupling to the incident electromagnetic field. In addition, the higher-frequency resonances experience stronger damping due to the increase of resistive losses in metal. Figure 4 shows a plot of plasmon resonance positions of the lowest two modes of the tetramer, which happen to be T 1 and T 2 , plotted against the gap-to-diameter ratio of gold tetramers in a dielectric environment of index n s = 1.4. Gold is assumed as the plasmonic material through the remainder of this article; the dielectric function of gold is modeled using interpolation of optical constants measured by various authors and compiled in reference [32]. The right vertical axis is labeled by vacuum wavelength λ vac corresponding to the real part of this dielectric function. The graphical insets of Fig. 4 show electric fields in cross-sections of these resonances. We note that these plasmon modes can also be constructed as linear superpositions of the modes in the bases for the irreducible representations of the T 1 and T 2 modes in Fig. 1.

Electromagnetic spectra of tetramer colloids
In the previous section we provided some insight into the electromagnetic properties of symmetric tetrahedral clusters using the quasi-static plasmon hybridization and electrostatic eigenvalue theories. In particular, we described plasmonic resonances that may have the right properties to provide enhanced electric and magnetic susceptibilities. However, in finite-sized clusters (not too small compared with the wavelength of light λ ≡ 2πc/ω), retardation effects become important. These include a shift in the resonant frequencies [18] and the excitation of resonant modes that do not possess an electric dipole moment. To predict the exact frequencies of these resonances and their strength, we have made finite-element frequency-domain (FEFD) electromagnetic simulations of tetramers using the commercial software package COMSOL. Extinction and absorption cross-section are measured as functions of frequency in the following fashion. A single tetramer is placed in a rectangular domain, on lateral sides of which either periodic or mirror-symmetry boundary conditions are applied. In effect, a doubly-periodic rectangular array (with periods L y , L z ) of identical tetramers is simulated. As long as individual tetramers interact only weakly via their near-field, and far-field interactions are not resonantly enhanced [21,22], spectra of ordered arrays are close to those of randomly distributed/oriented tetrahedral clusters. These conditions are fulfilled by allowing sufficient separation between tetramers, and by using wavelengths sufficiently longer than the largest of the two periods, ruling out Wood's anomalies. The array is illuminated by a monochromatic, linearly polarized plane wave of unit intensity, with E||ŷ and H||ẑ, incident normally ( k = kx) on the yz-plane of the array.  Complex amplitudes of transmitted (t) and reflected (r) waves are measured in the far field and interpreted as forward and backward scattering amplitudes, respectively. This allows one to define extinction σ ext = (1 − T )S 0 and absorption σ abs = (1 − T − R)S 0 cross-sections, where S 0 = L y L z is the cross-sectional area of one unit cell and T = |t| 2 and R = |r| 2 are energy transmission and reflection coefficients respectively. In the limit of small extinction, σ ext is related to the decay constant κ (with dimensions of inverse length) through the usual formula κ = σ ext n 0 , where n 0 = 1/V 0 is the number density of tetramers, and V 0 = S 0 L x ≡ L x L y L z is the specific volume per cluster. Indeed, if the distance between consecutive layers of scatterers is L x , then the wave intensity is damped by factor Optical spectra of tetramers made of 90 nm and 120 nm gold spheres are presented in Fig. 5  and 6. In both cases, two strong resonances are observed in the extinction and absorption spectra. The positions of these resonances are in rough agreement with electrostatic predictions (see Fig. 4) with the gap-to-diameter ratios extrapolated to 2/90 and 2/120, although red-shifted notably. This electromagnetic red shift phenomenon is well understood in terms of the sur-  The magnetic dipole resonance is related predominantly to the sextupole (if viewed twodimensionally) resonance of the three spheres lying in the plane orthogonal to the incident magnetic field H 0 = H 0ẑ . This mode is similar to the A 2 plasmon of the trimer studied in a previous paper using plasmon hybridization [25]. The electric field in this resonance circulates around the center of the trimer, causing a non-vanishing magnetic flux Φ = B z dx dy through the z = const plane. Since sextupolar symmetry corresponds to an azimuthal number of M J = 3, it is clear that this resonance is also an electric octupole (J = 3, M J = 3) in three dimensions. The field picture of the electrostatic T 1 -symmetric resonance in the inset of Fig. 4 supports our conclusion that the lowest-frequency resonance in optical spectra (Fig. 5, 6) is related to the lowest T 1 resonance.

Effective permittivity and permeability of tetramer colloids
Periodically-arranged cluster arrays are easier to characterize in electromagnetic simulations than random suspensions, because a simple procedure exists for retrieving effective medium parameters from the amplitude and phase of reflected and transmitted waves scattered off a metamaterial slab [35,36,37]. Yet, in effective medium regime and for sub-wavelength interparticle distances we expect that spectral features related to single-cluster resonances are common for random and periodic ensembles. To complete our identification of electric and magnetic resonances of a tetramer, we have applied that standard procedure and evaluated ε eff and µ eff of a cluster array.
A slightly different orientation of a tetramer has been chosen for these calculations: it was determined that orienting the two opposite, orthogonal edges of a tetrahedron along the major axes of rectangular cluster array minimizes splitting of single-cluster resonances due to the nearest-neighbor interactions between clusters of the lattice. Such choice of orientation makes two of the tetramer symmetry planes compatible with those of rectangular lattice. Since a tetrahedron does not have three mutually orthogonal planes, this array lacks a central symmetry plane. Fortunately, homogenization methods for asymmetric structures have been recently developed [36].  Effective medium parameters of slabs of periodic metamaterials, extracted using the standard homogenization method, are known to exhibit some unusual (non-Lorentzian) behavior [37], the origin of which had caused some debate in the past [38]. This behavior is believed to be associated with spatial dispersion in periodic structures [38]; examples of non-Lorentzian behavior include bands with negative Im ε eff in the vicinity of a magnetic resonance, or negative Im µ eff around an electric resonance. In addition, real parts of ε eff (µ eff ) show "reversed" Lorentz-shaped kinks in magnetic (electric) resonances. These "anti-resonances" do not violate any laws of physics: for example, the overall medium response remains passive [38]. Another consequence of dealing with a periodic array is the fact that the effective permeability determined using the scattering procedure [35] is forced to go to zero near an electric resonance [37]. Indeed, the refractive index n eff = √ ε eff √ µ eff equals k Bloch c/ω, but k Bloch is limited by the size of the Brillouin zone [37]. A large rise of ε eff (see Fig. 8) near electric resonance thus inevitably causes µ eff to go down, as seen in Fig. 9. For the purpose of demonstrating optical magnetism in metafluids, we focus on the vicinity of the magnetic resonance frequency and disregard the above-mentioned artifacts.  For illustration purposes, we have made numerical simulations with an extremely small (1 nm) gap between spheres. Such a small gap serves to maximize the frequency separation of the electric and magnetic resonances and to minimize the effect of the former on the weak magnetic resonance. From Fig. 9, it is clear that there is a regular, Lorentz-shaped magnetic resonance at λ = 890 nm, somewhat distorted by an adjacent anti-resonance band associated with the strong electric dipole resonance, which begins at 850 nm. To make this magnetic resonance perfectly clear, we have also done these simulations with losses in gold reduced tenfold with respect to their true values. The dash-dotted curve on Fig. 9 demonstrates that with low losses, negative permeability would be possible in a metafluid with the cluster number density of 1/(0.0115 µm 3 ), in which clusters occupy 13% of the volume.
Although suspensions of solid gold particles in water or ethanol aggregate due to strong van der Waals forces at much lower concentrations, typically 0.001 − 0.01% vol., colloids with ∼ 100 nm-sized silica particles are stable with volume fractions ∼ 10% and even up to 50% [39, 40], which is close to the liquid-solid phase transition in a hard-sphere system [41]. We expect that silica-covered gold nanoshells can be concentrated to 1% vol. and above. Ex-perimental studies of such nanoclusters, their formation and stability, are being conducted and will be reported elsewhere. Note that if fluidity of a negative permeability metamaterial is not a requirement, one can simply condense tetrahedral APMs to densities approaching the closepacking density (≈ 64% for random packing, ≈ 74% for hexagonal or cubic packing) to create a solid metamaterial with isotropic permittivity and permeability. Achieving µ eff < 0 or ε eff < 0 at such high volume fractions is a much easier task than in liquid form. Thereupon plasmonic metafluids can be used as precursors to isotropic doubly-negative metamaterials operating in effective medium regime; fabrication of photonic crystals and other non-subwavelength metamaterials using colloid condensation has been previously demonstrated [40,10].
Potentially, the strength of the magnetic resonance, characterized by the normalized magnetic moment, (9), can be increased by utilizing nanoparticles with more complicated shapes, including, for instance, non-concentric and non-spherical nanoshells [20]. In addition, recent progress in fabrication of crystalline SiC nanoparticles [42] provides hope that negative-ε nanoparticles with losses an order of magnitude smaller than those in gold [43,44,45,24] can be utilized to create metafluids with optical magnetism in the mid-infrared range. We observe from Fig. 9 that this 10-fold reduction of losses is sufficient for achieving µ eff < 0 even without further enhancement of the magnetic response through particle shape engineering or demanding higher cluster density.
Resistive damping and colloid stability are not the only problems that an experimental demonstration of negative permeability in metafluids will face. Additional challenge comes from strong sensitivity of resonant frequencies to the gaps between particles comprising an APM. The width of plasmon resonances in plasmonic nanostructures is typically on the ∆λ E ∼ 100 nm scale for electric and ∆λ M ∼ 10 nm scale for magnetic resonances, as seen from Figures 8,9. From Fig. 4 we can now estimate the allowed width of gap distribution in weakly polydisperse colloids. For example, if the gap-to-diameter distribution is centered at 0.06 (which corresponds to a magnetic resonance at λ ∼ 650 nm if retardation red-shift can be ignored), then in order to have magnetic resonances of most particles within the ±10 nm range around that wavelength, the gap/diameter ratios must be contained in the 0.05 − 0.07 interval. For 100 nm spheres, this translates into a ±1 nm requirement for gap distributions. Although such unprecedented uniformity of gaps is almost impossible to achieve with solid gold particles, it may be possible if thin gold nanoshells are grown on silica nanoparticles, which can be created in almost perfectly spherical shapes and with very smooth surfaces [9]. Numerical simulations indicate that effective medium parameters of metafluids with gold nanoshells differ very little from those with solid spheres, as long as the shell is thicker than the skin depth in metal (∼ 20 − 25 nm).
In this Section, we have shown theoretically that colloidal solutions of plasmonic nanoclusters can have both negative dielectric permittivity and negative magnetic permeability. However, the most exciting applications, such as the negative index metafluid, demand ε eff and µ eff simultaneously negative. Achieving ε eff < 0 and µ eff < 0 at some frequency requires one more design step: the positions of the strongest electric (ED) and magnetic dipole (MD) resonances must be engineered such that ω MD > ω ED . In that case, MD resonance could be placed within the narrow ε eff < 0 band above ω ED . This could be accomplished by adding to the metafluid plasmonic core-shell particles that have been shown [16] to exhibit red-shifted ED resonance. Alternatively, re-ordering of ED and MD resonances in a tetramer can be accomplished by letting plasmonic spheres touch each other in the cluster. More sophisticated choices of particle and cluster geometry and topology are possible and may be utilized in negative-index metafluid engineering.

Conclusions
The concept of metafluids has been introduced. Metafluids consist of Artificial Plasmonic Molecules (APMs) suspended in a fluid. We have concentrated on APMs in the shape of tetrahedral plasmonic nanoclusters. We have investigated the properties of plasmonic tetramers using electrostatic methods aided by group theory and fully electromagnetic finite-element simulations. We found that in the quasi-static approximation, electric response is dictated only by T 2 -symmetric plasmons and magnetic response only by T 1 plasmon states, which are excited by electric field as electric octupoles due to retardation effects. The electric and magnetic response of the tetramer allows one to construct an effective medium with a completely isotropic electric and magnetic response. Electromagnetic simulations indicate that achieving ε eff < 0 and µ eff < 0 in colloidal solutions of "artificial molecules" should be possible using either sufficiently high concentrations of gold clusters or materials with low-loss negative permittivity.