# Accounting for inhomogeneous broadening in nano-optics by electromagnetic modeling based on Monte Carlo methods

^{a}School of Engineering and Applied Sciences, Harvard University, Cambridge, MA 02138;^{b}State Key Laboratory of Supramolecular Structure and Materials, College of Chemistry, Jilin University, Changchun 130012, People’s Republic of China;^{c}Department of Chemistry, University of Toronto, Toronto, ON, Canada M5S 3H6;^{d}Department of Chemistry and Biochemistry, University of Maryland, College Park, MD 20742;^{e}Department of Chemical Engineering and Applied Chemistry, University of Toronto, Toronto, ON, Canada M5S 3E5; and^{f}The Institute of Biomaterials and Biomedical Engineering, University of Toronto, Toronto, ON, Canada M5S 3G9

See allHide authors and affiliations

Contributed by Federico Capasso, December 20, 2013 (sent for review July 25, 2013)

## Significance

The advent of nanotechnology has enabled the study of physical phenomena in structures with nanoscale dimensions. Pushing the limits of fabrication techniques inevitably leads to uncertainties—for example, an array of nanoscale resonators may be designed to be identical, but in fact has a distribution of sizes due to fabrication imperfections, leading to a distribution of responses. In optical experiments involving collections of nanostructures these effects are often acknowledged but rarely quantified. We demonstrate a technique that combines electromagnetic simulations with a Monte Carlo sampling technique to rigorously account for “inhomogeneous broadening” of optical resonances as a result of fabrication or synthesis imperfections. This approach has wide applicability to any experiments involving collections of structures designed to be identical.

## Abstract

Many experimental systems consist of large ensembles of uncoupled or weakly interacting elements operating as a single whole; this is particularly the case for applications in nano-optics and plasmonics, including colloidal solutions, plasmonic or dielectric nanoparticles on a substrate, antenna arrays, and others. In such experiments, measurements of the optical spectra of ensembles will differ from measurements of the independent elements as a result of small variations from element to element (also known as polydispersity) even if these elements are designed to be identical. In particular, sharp spectral features arising from narrow-band resonances will tend to appear broader and can even be washed out completely. Here, we explore this effect of inhomogeneous broadening as it occurs in colloidal nanopolymers comprising self-assembled nanorod chains in solution. Using a technique combining finite-difference time-domain simulations and Monte Carlo sampling, we predict the inhomogeneously broadened optical spectra of these colloidal nanopolymers and observe significant qualitative differences compared with the unbroadened spectra. The approach combining an electromagnetic simulation technique with Monte Carlo sampling is widely applicable for quantifying the effects of inhomogeneous broadening in a variety of physical systems, including those with many degrees of freedom that are otherwise computationally intractable.

In photonics experiments and applications, frequent use is made of ensembles of individual structures operating as a single whole; these include, for example, lithographically defined arrays of metallic nanostructures that form frequency-selective surfaces (1), metasurfaces (2⇓–4) or sensor arrays (5), colloidal solutions or suspensions (6, 7), randomly dispersed nanoshells, quantum dots or nanocrystals on a substrate (8), and many others.

Assuming that the elements in the ensembles are independent (i.e., they do not experience significant near- or far-field coupling), an assumption that can often be made in sparse, disordered systems, the optical response of these ensembles is simply the sum of the response of all of their constituents. In the case that such an ensemble is composed of many identical elements, its spectral response should be the same as that of each individual element. In real systems, however, the constituent elements are never precisely identical: Any fabrication or synthesis technique including top–down lithography and bottom–up self-assembly will introduce a distribution of geometrical parameters (also known as polydispersity) that leads to inhomogeneous broadening in the spectral features of the total ensemble (e.g., refs. 9⇓⇓⇓⇓⇓⇓–16). To avoid inhomogeneous broadening effects in experiments, complex techniques are sometimes used to measure the optical response of individual elements (17). Other times inhomogeneous broadening can be helpful, for example in situations where a broadband optical response is desired such as in photovoltaic applications (18).

Although full-wave electromagnetic simulations are often used to model and understand optical systems that cannot be described analytically (e.g., ref. 19), these methods cannot easily account for polydispersity that leads to inhomogeneous broadening. This issue is sometimes addressed by artificially increasing the damping constant of materials (15), but this approach yields only nonspecific, qualitative information, does not provide a way to distinguish between the various sources of polydispersity (geometrical or material), and is in general not physical.

In the present work, we demonstrate that a complex polydisperse ensemble of noninteracting elements can be fully modeled using a Monte Carlo approach (20, 21), using finite-difference time-domain (FDTD) simulations for the intermediate steps. Monte Carlo methods combined with electromagnetic calculations have previously been applied to problems in electromagnetics such as scattering from random rough surfaces (21⇓⇓–24) and light transport through tissues (25). Here we predict the extinction spectra of self-assembled gold nanorod chains (“nanopolymers”) suspended in a solution. This physical system has a large number of degrees of freedom (e.g., the lengths and widths of the individual nanorods comprising the chains, the total number of rods composing each chain, the gaps between the rods, and their orientation, etc.) and is therefore a particularly challenging model system.

## Model System: Gold Nanopolymers in Solution

Recent experiments have demonstrated that gold nanorods end tethered with polystyrene ligands can undergo self-assembly in solution and form linear (or bent) chains, in a process analogous to step-growth polymerization (26⇓–28). In this process, individual nanorods with an end tether on each end behave as monomers (Fig. 1*A*). In a colloidal polymer (a nanopolymer), the nanorods are the repeat units and the tethers between the nanorod ends act as bonds.

At a particular stage of self-assembly, the concentration of unreacted functional groups [*M*] is twice as large as the concentration of nanorod chains in the system (including individual unreacted nanorods), because each chain has two ends. The average degree of polymerization is defined aswhere [*NR*]_{0} is the initial concentration of nanorods in the solution. For the self-assembly time *t* = 0, , and the initial concentration of functional groups is [*M*]_{0} = 2[*NR*]_{0}.

The self-assembly occurs as follows: The first step is the reaction between two individual nanorods to form a dimer; the dimer can then react with a monomer to form a trimer or with another dimer to form a tetramer, and so forth (28). This process yields a mixture of chains comprising various numbers of nanorods *x*. The degree of polymerization of the entire mix of nanorods and nanorod chains at a particular moment in time can be quantified by the “conversion” parameter *p*, defined as the fraction of end tethers that have reacted, such thatThe conversion *p* is then related to byFor this type of step-growth polymerization, the concentration of chains containing *x* nanorods, *c*_{x,p}, can be predicted by the Flory (or “most probable”) distribution given a particular conversion *p* as (28, 29)A solution of gold nanorods and nanorod chains (*x*-mers) can be viewed as a “metamaterial fluid” or “metafluid” and will have different optical properties from that of the solvent alone as a result of the resonant scattering contribution of the *x*-mers. In the language of metafluids, the *x*-mers can be viewed as “artificial plasmonic molecules” (7) suspended in a liquid, and one can envision a characterization experiment in which the extinction spectrum of the fluid is measured using a broadband optical source and a spectrometer (Fig. 1*B*). If the solution is dilute and there is little agglomeration of *x*-mers, the extinction spectrum of the solution can be calculated as the sum of the extinction spectra of all of the individual *x*-mers, which can be predicted by a variety of full-wave electromagnetic simulation techniques (30). Because the solution is dilute, any multiple scattering effects can be neglected, and the spectra can be summed incoherently (i.e., neglecting phases) because the positions and orientations the individual nanorod chains are random and constantly changing throughout the solution via thermal motion. Nonetheless, the problem remains very challenging because of the large number of degrees of freedom: A solution can contain chains of nanorods of nearly any length, and every nanorod can differ in its geometrical parameters.

## Modeling Nanorod Chains Comprising Identical (Monodisperse) Nanorods

The extinction spectrum of a particular nanorod chain is determined by the length (*L*) and diameter (*d*) of the individual nanorods composing the chain, the inter-nanorod distance (*l*), and the number *x* of nanorods in the chain (see, e.g., ref. 31). For simplicity, we assumed that the chains remain linear (no bending). To model the spectrum of an individual nanorod chain, we first assumed that the values of *L*, *d*, and *l* are constant for all nanorods composing the chain and then examined the relationship between *x* and the normalized extinction spectrum *ε*_{x}(*λ*) of an individual chain with a particular, well-defined aggregation number *x*, as *x* increased from 1 to 10. *ε*_{x}(*λ*) is determined as *ε*_{x}(*λ*) *= σ*_{x}(*λ*)*/x*, where *σ*_{x}(*λ*) is the extinction cross section of the chain, and represents the extinction spectrum of a single chain normalized to its length.

We performed full-wave 3D FDTD simulations, using the total-field scattered-field (TFSF) formulation (32), implemented in a commercial software package (Lumerical). In the simulations, we set *L* = 52 nm, *d* = 13 nm, and *l* = 6.7 nm, the values corresponding to those in self-assembly experiments. The mesh size of the simulations was 0.5 nm such that all features were well resolved. We used a background index of refraction of 1.42, corresponding to typical solvents used in nanorod self-assembly [*N*,*N*-dimethylformamide (DMF)—water mixture with a water content of 15 wt% (33)], and an index of refraction of 1.57 for the polystyrene ligands (34). The incident light was set to be polarized with the electric field along the long axis of the chain as that is the orientation of the dominant dipole moment of the chains (further discussion in *SI Text*).

The simulated normalized extinction spectra *ε*_{x}(*λ*) of these nanorod chains with *x* from 1 to 10 are plotted in Fig. 2*A*. As *x* increased, the localized surface plasmon resonance (LSPR) peak experienced a gradual red shift from ∼890 nm to ∼1,020 nm. This type of red shift has been predicted and observed in plasmonic dimers, trimers, and longer chains and is generally attributed to a combination of capacitive near-field coupling between the neighboring nanorods and retardation effects that set in when the size of the chain becomes nonnegligible compared with the wavelength (31, 35⇓⇓–38). The effect was particularly strong as *x* increased from 1 to 4, as the majority of nanorods forming the chain acquired new nearest neighbors, but then quickly saturated for longer chains.

Next, we modeled the extinction spectra of a population of monodisperse nanorod chains at various conversions *p*. For dilute colloidal solutions, the total extinction, *Ext*_{tot,p}, is the sum of extinctions of the individual chains. Thus, for noninteracting, monodisperse nanorod chainswhere *b* is the interaction path length, *σ*_{x}(*λ*) is the extinction cross section of an individual nanorod chain; *c*_{x,p} is the concentration of the nanorod chains with a particular *x* at conversion *p*; and the subscripts 1, 2, … , *x* refer to the number of nanorods in the chain.

By inserting into Eq. **5** the extinction cross sections obtained from the FDTD simulations and *c*_{x,p} obtained from Eq. **4**, we calculated the extinction spectra of the entire population of nanorod chains at various values of *p* (corresponding to particular self-assembly times *t*). Extinction spectra for *p* of 0.3, 0.5, 0.7, and 0.9 are shown in Fig. 2 *B*–*E*, respectively. The colored lines show the contribution to the overall extinction spectrum by the nanorod chains with a particular *x*, that is, *bσ** _{x}*(

*λ*)

*c*

_{x,p}, and the black lines show the total extinction of all of the chains (

*Ext*

_{tot,p}(

*λ*)).

## Monte Carlo Modeling of Polydisperse Nanorod Chains

In an experimental setting, synthesized gold nanorods always exhibit polydispersity. In the present work, we assume that the nanorods have lengths and diameters with distributions of *L* = 52 ± 6.1 nm and *d* = 13 ± 1.6 nm, respectively, and the distance between the ends of nanorods in the self-assembled chains is *l* = 6.7 ± 1.4 nm [see *SI Text* for transmission electron microscopy (TEM) images from which these values are inferred]. Empirically (and as a consequence of the central limit theorem), these distributions are approximately Gaussian.

To model the extinction spectrum of a collection of nanorod chains with geometrical variations, one could in principle perform an exhaustive set of FDTD simulations, sweeping over all possible nanorod lengths and diameters, as well as over all of the possible gap lengths between adjacent nanorods in the chains, and weigh the resulting spectra appropriately to predict the expected LSPR spectrum of the ensemble (in *SI Text*, we show this type of calculation for nanorods with just two polydisperse parameters using a semianalytical approach). However, even for a modest number of nanorod constituents of each chain, this parameter space explodes, making this computational problem intractable. To overcome this, we used a strategy that combines the Monte Carlo method with FDTD simulations, which is graphically described in Fig. 3. For a chain comprising *x* nanorods, we assumed that the geometrical parameters *L*, *d*, and *l* for each nanorod and gap are independent and identically Gaussian distributed throughout the chain. Accordingly, we selected the geometrical parameters for each nanorod and each gap stochastically from the appropriate empirical Gaussian distribution (Fig. 3*A*) and then used FDTD simulations to calculate the normalized LSPR spectrum for the nanorod chain (Fig. 3*B*). We iterate this process until a relatively smooth, invariant distribution emerges from the average of the simulated spectra and then fit this average spectrum to a Gaussian distribution to obtain an estimate of the average extinction spectrum *σ*_{x}′(*λ*) (Fig. 3*C*). For this work, we performed 250 simulations for the monomers, 150 simulations for the dimers, 90 simulations for the 3-mers, and 60 simulations each for chains comprising 4–10 nanorods. This number of simulations was enough to demonstrate the effects of inhomogeneous broadening and was a compromise between accuracy and total computation time. The resulting Gaussian fits to the normalized extinction spectra are shown in Fig. 4*A*.

Note that we expect this spectrum of an ensemble of polydisperse nanorod chains to resemble a Gaussian more than a Lorentzian distribution (as would be expected for a single, isolated resonance) as a consequence of a general correlation between resonance peak wavelengths and the overall lengths of their corresponding nanorod chains (39). Because the nanorod chain lengths are Gaussian distributed, the resonance peak wavelengths tend to be Gaussian distributed as well. Because the widths of these Gaussian distributions of resonance peaks tend to be greater than the widths of the individual resonances, we expect that the ensemble spectrums will have more Gaussian than Lorentzian character. This is verified for the special case of unreacted polydisperse nanorods in *SI Text*.

As *x* increased from 1 to 10 we again see a red shift of the LSPR peak from ∼850 nm to ∼1,050 nm (Fig. 4*A*). The normalized extinction spectra of individual chains comprising polydisperse nanorods are substantially broader than their monodisperse counterparts. This broadening is not a result of increased absorption or scattering losses that would also lead to a broader peak in extinction due to decreased quality factors, but is instead a result of inhomogeneous broadening. Although it is evident from Fig. 4*A* that the normalized LSPR peak heights and locations have not fully stabilized (more simulations would be necessary for the results to fully converge), we can still clearly see the trends in normalized LSPR extinction spectra as *x* increases. For example, whereas in monodisperse chains the normalized LSPR peak heights tend to decrease slightly as *x* increases (Fig. 2*A*), the opposite is true for polydisperse chains. This is because the polydispersity has a stronger effect on the chains comprising fewer nanorods than on the longer chains: In the longer chains the small variations in the individual rods tend to cancel out; this means that the normalized spectra of the ensemble of shorter polydisperse chains tend to be broader and correspondingly smaller in amplitude compared with ensembles of longer chains. In this way, the Monte Carlo approach yields a qualitatively different prediction from that obtained from a monodisperse approximation.

By inserting the averaged and fitted values of *ε** _{x}*(

*λ*)′ (obtained from the Monte Carlo simulations) and

*c*

_{x,p}(obtained from Eq.

**4**according to the Flory distribution) into Eq.

**5**, we calculated the extinction spectra of the entire population of polydisperse nanorod chains at various values of conversion

*p*. Extinction spectra for

*p*values of 0.3, 0.5, 0.7, and 0.9 are shown in Fig. 4

*B–E*, respectively, along with the contributions of individual populations of nanorod chains of various

*x*. The total extinction spectrum of all chains for

*p*ranging from 0 to 0.9 is shown in Fig. 5

*B*, compared with the same calculation using monodisperse chains shown in Fig. 5

*A*. In the monodisperse case, the total extinction spectra of intermediate

*p*values are distinctly bimodal due to the relatively narrow LSPR peak widths associated with individual nanorod chains and the relatively large LSPR peak red shifts associated with increasing

*x*. In the polydisperse case, the relatively broad LSPR peak widths that result from inhomogeneous broadening wash out much of this bimodal spectral feature. Instead, the polydisperse spectra each feature a single broad peak that slowly red shifts with increasing self-assembly time (and hence conversion

*p*).

In this paper we intentionally make no comparison with experimental optical data. Although our simulation method effectively accounts for many key physical effects contributing to the optical response of self-assembled nanopolymer solutions including the polydispersity of nanorod widths and lengths and gaps between nanorods, it does not account for chain bending or retardation effects for chains that are not oriented perpendicular to the incident light. Our model system of self-assembled nanopolymers in solution is a particularly challenging system to simulate due to the overall number of degrees of freedom and the computational resources required for every full-wave 3D simulation of large nanorod chains with small features (such as gap sizes) that must be well resolved. Despite this, we believe that the Monte Carlo approach combined with electromagnetic simulations (or analytical calculations) that we demonstrate here may be the most efficient method that provides meaningful information about inhomogeneous broadening in optical systems, especially when applied to slightly simpler systems such as lithographically defined or self-assembled nanostructures on a substrate.

## Conclusion

In conclusion, we used FDTD simulations combined with a Monte Carlo approach to study the effects of inhomogeneous broadening on the extinction spectra of populations of gold nanorods formed in solution. We found that inhomogeneous broadening due to dispersion in the geometrical parameters of the nanorods (lengths and widths) and the gaps between neighboring rods significantly affected the shape and bandwidth of the resonance spectra of solutions of nanorod chains. More generally, we conclude that in systems involving large collections of independent resonant elements, inhomogeneous broadening introduces significant differences between the resonant responses of individual elements and the ensemble. To account for such differences, it is possible to run separate calculations or simulations for every possible set of geometrical parameters and then perform a weighted average; however, as in the present demonstration, this is often an intractable problem, especially for structures with many degrees of freedom and resource-expensive numerical techniques such as FDTD. This, however, can be overcome by using a Monte Carlo approach consisting of iterated stochastic sampling from the entire parameter space combined with numerical simulations. In the present demonstration, 3D full-wave FDTD simulations are used, but in principle any analytical, semianalytical, or fully numerical method can be applied.

## Acknowledgments

M.A.K. acknowledges helpful discussions with N. Yu, R. Blanchard, and J. Fan. The Lumerical FDTD simulations in this article were run on the Odyssey cluster supported by the Harvard Faculty of Arts and Sciences Division Research Computing Group. This work was supported in part by the Air Force Office of Scientific Research under Grant FA9550-12-1-0289. E. K. thanks the Natural Sciences and Engineering Research Council (Canada) for supporting this work by a Discovery grant and a Strategic Network for Bioplasmonic Systems Biopsys grant. M.A.K. was supported by the National Science Foundation through a Graduate Research Fellowship.

## Footnotes

↵

^{1}H.G. and M.A.K. contributed equally to this work.- ↵
^{2}To whom correspondence should be addressed. E-mail: capasso{at}seas.harvard.edu.

Author contributions: H.G. and M.A.K. designed research; H.G., M.A.K., and K.L. performed research; K.L., Z.N., and E.K. contributed new reagents/analytic tools; H.G., M.A.K., K.L., and F.C. analyzed data; and H.G., M.A.K., K.L., and F.C. wrote the paper.

The authors declare no conflict of interest.

This article contains supporting information online at www.pnas.org/lookup/suppl/doi:10.1073/pnas.1323392111/-/DCSupplemental.

## References

- ↵
- Munk BA

- ↵
- ↵
- Kats MA,
- et al.

- ↵
- ↵
- Yanik AA,
- et al.

- ↵
- Elghanian R,
- Storhoff JJ,
- Mucic RC,
- Letsinger RL,
- Mirkin CA

- ↵
- ↵
- ↵
- Hanarp P,
- Käll M,
- Sutherland DS

- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- Robert CP,
- Casella G

- ↵Sadiku MNO (2009)
*Monte Carlo Methods for Electromagnetics*(CRC, Boca Raton, FL). - ↵
- ↵
- ↵
- ↵
- ↵
- Liu K,
- et al.

- ↵
- ↵
- Odian G

- ↵
- Rubinstein M,
- Colby RH

- ↵
- ↵
- ↵
- Taflove A,
- Hagness SC

- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵

## Citation Manager Formats

## Article Classifications

- Physical Sciences
- Applied Physical Sciences