Wulff-Based Approach to Modeling the Plasmonic Response of Single Crystal, Twinned, and Core-Shell Nanoparticles.

The growing interest in plasmonic nanoparticles and their increasingly diverse applications is fuelled by the ability to tune properties via shape control, promoting intense experimental and theoretical research. Such shapes are dominated by geometries that can be described by the kinetic Wulff construction such as octahedra, thin triangular platelets, bipyramids, and decahedra, to name a few. Shape is critical in dictating the optical properties of these nanoparticles, in particular their localized surface plasmon resonance behavior, which can be modeled numerically. One challenge of the various available computational techniques is the representation of the nanoparticle shape. Specifically, in the discrete dipole approximation, a particle is represented by discretizing space via an array of uniformly distributed points-dipoles; this can be difficult to construct for complex shapes including those with multiple crystallographic facets, twins, and core-shell particles. Here, we describe a standalone user-friendly graphical user interface (GUI) that uses both kinetic and thermodynamic Wulff constructions to generate a dipole array for complex shapes, as well as the necessary input files for DDSCAT-based numerical approaches. Examples of the use of this GUI are described through three case studies spanning different shapes, compositions, and shell thicknesses. Key advances offered by this approach, in addition to simplicity, are the ability to create crystallographically correct structures and the addition of a conformal shell on complex shapes.

Introduction

Plasmonic nanoparticles (NPs) have gained much attention in the scientific community owing to their optical properties that can be exploited for a variety of applications, ranging from sensing[1] and photocatalysis,[2] to biomedicine[3] and optical circuits.[4] NPs of free electron metals confine light via collective electron cloud oscillations triggered by an incident oscillating electromagnetic field, giving rise to resonances known as localized surface plasmon resonances (LSPRs). LSPRs enhance light scattering and absorption whilst amplifying local electric fields at the NP’s surface. Commonly, plasmonic NPs are synthesized from Au and Ag[5] but novel plasmonic structures of earth-abundant materials such as Al, Cu, and Mg have recently been demonstrated theoretically and experimentally.[6−11] Other metals, such as Ga and In also present significant LSPR tunability[12,13] while Pd and Pt sustain rather weak and broad LSPRs.[14,15]

The energy and peak width of a LSPR can be tuned by controlling the composition, environment, size, and shape of the NPs, to name a few.[5,16,17] Shape is particularly appealing, as it easily and predictably controls the near-field distribution around a particle, creating for instance localization around either the corners or faces in a cube[18] or tip and shaft in a rod,[19] depending on the resonance frequency. Shape tuning can be supplemented by composition tuning by incorporating layers of different materials, either as a simple core–shell structure or complex multishell, egg-yolk, or other shapes.[20] Such core–shell structures are of particular interest not only because they introduce new parameters that affect LSPRs, such as shell composition and thickness,[21] but also because they can combine plasmonic and nonplasmonic materials,[22] thus providing further means to design functional NPs. Core–shell structures may also be used to prevent the oxidation of a core,[23] or occur spontaneously upon self-limiting oxidation of a metal.[7,9]

The shape of NPs, so critical for their optical properties, is dictated by the crystal structure and growth environment. At thermodynamic equilibrium (e.g., in vacuum or another scarcely interacting environment), NP shape can be predicted analytically from the (thermodynamic) surface free energy according to the Wulff construction.[24] Briefly, the distance normal to an (hkl) facet, h, is related to its surface energy γ aswhere Λ is a constant accounting for volume. While this is only valid for a free-floating particle during growth, extensions to the model including interaction with one or two interfaces with a substrate have been developed, named the Winterbottom and the Summertop constructions, respectively.[25,26] Similarly, the addition of internal rather than external boundaries to account for twinning leads to the modified Wulff construction.[27] Realizing that the thermodynamic shapes were rarely present in reaction products, the kinetic Wulff construction[28] was developed, where a growth velocity v is used instead of the thermodynamic surface free energy γ. The kinetic approach to the growth of twinned structures was then, recently, developed,[29] where twin boundaries, disclinations, and re-entrant surface kinetic effects help explain most experimentally obtained shapes for face-centered cubic (fcc) materials, which include most plasmonic metals (Cu, Ag, Au, Al).

Here, this kinetic version of the regular (single crystal) and modified (twinned) Wulff construction[29] is used as the basis of a shape modeling code integrated in a user-friendly, standalone graphical user interface (GUI). Briefly, to derive a NP’s shape from the inputs of v and optional enhancements at re-entrant surfaces, twin boundaries, and disclinations, space is discretized in a cubic three-dimensional grid and growth velocities (vg) are calculated at each point of the grid, described as a vector p⃗ with respect to the center of the shape. This is mathematically implemented by calculating the following expression on each grid point p⃗where b is a smoothing factor and n⃗ = Fvn̂, where n̂ is the unit vector of the corresponding crystallographic facet, and F is the enhancement factor given by

The Wulff shape is then defined as an isosurface of growth velocities, because this is directly proportional to the distance from the geometric center of the particle to the facet (eq ). In the case of twinned NPs the created shape is mirrored along the twin plane for NPs with a single twin plane or rotated around the five-fold symmetry axis for NPs with five nonparallel twin planes.[30]

Modeling shape effects is the key to understand how geometry affects LSPRs (both far and near-field) and can support and inspire the design and synthesis of NPs for tailored light–matter interactions. To do this, one must solve Maxwell’s equations, which describe the electromagnetic interactions at play. Because analytical solutions to Maxwell’s equations are limited to a small number of simple geometries, such as spheres with the Mie solution[31] or ellipsoids with the Mie–Gans solution,[32] various numerical techniques have been developed to model arbitrary shapes. Prevailing approaches include the finite difference time domain (FDTD) method,[33] the discrete-dipole approximation (DDA),[34] the finite-element method (FEM)[35] and the boundary element method (BEM).[36] The first three approaches (FDTD, DDA, FEM) require discretization over the NP volume while for BEM the discretization is applied only to the NP surface.

In DDA,[34] particles are represented by an array of small cubic elements considered dipoles interacting with each other and with the incident electric field. These interactions result in a system of Maxwell’s equations that can be solved to obtain the polarization of each dipole and subsequently to calculate the absorption and scattering properties of the particle, as well as near-field effects such as field enhancement mapping and local charge distribution around nanostructures. The latter can give important information about LSPR modes such as their localization on the edges and corners of triangular plates.[37] In FDTD[33] the space and time derivatives that appear in Maxwell’s equations are replaced by finite differences, therefore requiring a discretization over both time and space, the latter achieved by a grid of cuboid elements; the problem is then solved iteratively until a steady-state solution is achieved, where the error is better defined than in DDA.[38,39] In the case of FEM, space discretization is achieved using elements, usually tetrahedral, for which the Helmholtz equation is satisfied along with appropriate conditions to ensure continuity and a consistent solution.[35]

The DDA is a hugely successful and popular method because in general, it requires comparatively low computational power, depending of course on the dipole number and interdipole distance.[39,40] Unlike other techniques, it uses a simple and straightforwardly physically meaningful discretization of space in dipoles. One downside to this simplicity is that equally sized cubic elements do not allow for a denser, better fitting grid for curved surfaces, making it difficult to model high aspect ratio structures for instance.[41] When modeling NP shapes, for DDA and for the aforementioned computational techniques, the appropriate geometry input must be created, which can prove difficult, especially for shapes with complex features or with many facets and angles such as a Marks decahedron. Here we present an approach to solve this struggle and facilitate simulation of the plasmonic properties of various NP shapes.

Acknowledging the advantages of the DDA technique and the already validated Wulff construction theory, we incorporate the modified kinetic Wulff construction code[29] in a GUI that creates a crystallographically correct NP shape and all the appropriate inputs for DDSCAT,[34] an open source code that uses the DDA method to calculate the optical properties of nanostructures. We show the modeling of single crystal and twinned fcc NPs with and without a shell for both concave and convex geometries, in both the kinetic and thermodynamic regime. Below we first describe the GUI, and then demonstrate its capabilities by modeling and calculating the absorption and scattering properties of Au, Ag, and Al NPs of various shapes as well as core–shell structures including Au@SiO2 decahedra, Ag@SiO2 cubes, and Al@Al2O3 bipyramids with various shell thicknesses.

Computational Details

Au and Ag refractive indices (RI) were obtained from Johnson and Christy,[42] those of Al and Al2O3 from Palik,[43] and that of SiO2 from Rodríguez-de Marcos et al.[44] For Al, RI from Palik was available only up to 190 nm and the DDA extrapolation was used for 150–190 nm. The ambient RI was set to 1 (vacuum), and electron surface scattering corrections on the RI were not deemed necessary because all NPs are sufficiently small, that is, they have an effective diameter greater than 30 nm.[45] Scattering (Csca) and absorption (Cabs) cross sections were calculated by multiplying the corresponding DDSCAT output efficiencies with the effective area (πaeff2). The shapes were modeled with an interdipole distance range of 0.3–2.6 nm depending on the total size, shape, and material of the NP (Tables S1 and S2), with no fewer than 105 dipoles, except for cubes, to ensure accuracy.[39] The convergence of the results with the number of dipoles was investigated only for the smallest NPs (Figures S2–S5) as larger NPs are expected to require fewer dipoles to get results with the same accuracy.[39] In all cases, the incident light is modeled having two orthogonal polarizations, an approach commonly used to mimic unpolarized light, and propagates along the highest symmetry axis of the particle, that is, perpendicular to a face for the cube, along the direction of the 4-fold axis for the octahedron, along the direction of the 5-fold axis for the sharp and Marks decahedron, and perpendicular to the twin plane for the bipyramid and triangle. The contribution of the two polarization components was also investigated for Au and Ag decahedra and was found to mainly influence the LSPR intensities with little effect on the LSPR energy. These results (Figure S6) along with the GUI input parameters used for the shape modeling (Table S3) and further computational details can be found in the Supporting Information.

Results

User Interface

Figure shows the process of calculating a sharp decahedron and generating the dipole array with the aid of the Wulff construction GUI, and the subsequent calculation of the absorption and scattering spectra, as well as electric field distribution with the DDA. The Matlab-based GUI, deployed as a standalone application, features a main window (Figure a) with input panels and action buttons that guide the user through the steps from modeling the NP’s shape to creating the shape file (shape.dat) and then the parameter file (ddscat.par) which are inputs for the DDSCAT simulation. Additionally, a RI file is to be provided by the user (RI files of Au, Ag, Al, SiO2, and Al2O3 available in the Supporting Information).

Figure 1

Wulff construction tool. (a) Main Wulff construction GUI window, (b) resulting dipole representation, (c) scattering spectrum, and (d) near field (E⃗2) distribution at the NP’s mid-height obtained from the DDA calculation for a sharp Au decahedron.

Specifically, the surface growth velocities (v) of the {100}, {110}, and {111} planes as well as the kinetic growth re-entrant surface (φre-entrant), twin (φtwin), and disclination (φdisclination) enhancements are specified in the basic parameters panel. Thermodynamic surface energies can be used instead of growth velocities, with no enhancement, for thermodynamic shapes. The user can select no twin plane (single crystal), one {111} twin plane (monotwin), as common in fcc structures, or five nonparallel {111} twin planes (pentatwin). NPs with a shell can be calculated, providing information about the thickness of the shell (shell percentage) and the curvature of the shape (concave or convex) is supplied. The shell generated is conformal, that is, it has the same geometry as the Wulff NP without a shell.

After all the parameters are defined the shape isosurface is calculated and displayed; at this stage the shape is dimensionless. In the DDSCAT the size of the studied target is introduced in the parameter file as the effective radiuswhere V is the volume of the target. To calculate V and the effective radius (calculate effective radius), the user selects two arbitrary points on the shape and inputs the distance between them (Figure S1b). Conveniently, the two points can be on the same or different facets, or on the shell or core of the particle, allowing, for example, defining the size of the particle based on the plasmon length.

When performing the DDSCAT calculations the choice of number of dipoles or equivalently the interdipole distance is important for the accuracy of the obtained result. In DDSCAT the interdipole distance is defined by the total number of dipoles and the total volume of the studied structure. In the GUI, the user can specify the interdipole distance (define interdipole distance) through a dialog box (Figure S1c); the code then readjusts the number of dipoles to fit both the volume and interdipole distance requirements. To achieve a good accuracy the interdipole distance must be small compared to both any structural length of the target and the wavelength of the incoming radiation.[34] The convergence of the results can be checked by manipulating the interdipole length for a given shape.

The first DDSCAT input file generated is the shape file (create shape file), which contains the array of dipoles that represent the NP as well as information about its orientation and composition. The second DDSCAT input file is the parameter file (create ddscat file), which includes details about the computational setup of the calculations, the material and effective radius of the target, the incident field, and the output files. A detailed description of the shape and parameter files and their variables can be found in the DDSCAT manual,[46] and a typical parameter file along with an extended description of the GUI can be found in the Supporting Information.

Case Study 1: NP Shape and Composition

Au and Ag NPs are dominant among plasmonic metals as they exhibit strong, tunable LSPRs throughout the visible and infrared region.[5] Unlike Ag, Au NPs are stable toward oxidation and biocompatible, enabling biomedical applications.[3] Alternatives to the rather expensive Ag and Au are becoming commonplace, Al being an example that is earth abundant and sustains LSPRs in the visible and ultraviolet (UV) region, providing opportunities for UV plasmonics[6] All three metals have well-established synthetic techniques[7,47,48] leading to a variety of single-crystal and twinned NPs[49] whose shape-dependent plasmonic properties can be predicted or confirmed via numerical simulations. Here, Au and Ag cubes, {100}-capped bipyramids, and decahedra are chosen as examples of single crystal, monotwin and pentatwin noble metal shapes, respectively, while cubes, octahedra, and {100}-capped bipyramids are chosen as representative Al shapes. The different structures studied, all generated with the Wulff construction tool (Table S3), are shown in Figure a. The effective radius of the NPs is 31 nm, corresponding to a cube with an edge length of 50 nm.

Figure 2

Wulff construction shapes where the black line shows the edge length (a) and the corresponding simulated scattering and absorption cross sections for (b) Au, (c) Ag, and (d) Al. Green, blue, and red solid lines correspond to cube, decahedron, and bipyramid, respectively; green and blue dotted lines correspond to the octahedron and Marks decahedron, respectively. The effective radius is 31 nm for all NPs.

The scattering and absorption cross sections of Au NPs as a function of wavelength are shown in Figure b. As expected, these rather small Au cubes and decahedra exhibit one, dipolar, LSPR peak in the region of 500–550 nm[50,51] while, given its higher anisotropy, the bipyramid features a red-shifted main peak and a second peak as a high energy shoulder.[52] The dominant peak for all shapes shifts toward higher wavelengths in the order: Marks decahedra, cubes, sharp decahedra, and bipyramids, a trend that reflects the combined effects of increasing plasmon length and anisotropy of the shapes. More specifically, given that the NPs have the same volume, the plasmon length, defined as the length over which the dipole oscillations take place,[17] increases in the order: cube, decahedron, and bipyramid. As the plasmon length increases, resonance occurs at higher wavelengths, that is, causing the noted red-shifts. This also explains the smaller peak wavelength of Marks decahedron compared to the sharp one. On the other hand, higher symmetry causes a blue-shift. The thicker Marks decahedron is more similar to a sphere further explaining its blue-shifted LSPR compared to the sharp decahedron and the cube. The peak intensities follow an increasing trend from marks decahedron, to cube, bipyramid, and sharp decahedron. Peak intensity increases with the plasmon length unless high anisotropy causes the appearance of a new peak or shoulder, consequently decreasing the highest peak intensity as observed in the case of the Au bipyramids. The near-field response for a Au sharp decahedron is also reported in Figure d.

The absorption and scattering cross sections for the Ag NPs are presented in Figure c. The six characteristic main LSPR modes of Ag cubes, well-identified in the literature,[53,54] span roughly 330–420 nm. The LSPR peaks are broader for the decahedra and the bipyramids as their symmetry leads to less mode degeneracy than the cube.[53,55] A red-shift trend in the peak wavelength follows the Au NPs pattern.

Finally, the LSPRs of Al NPs (Figure d) appear at lower wavelengths than Ag and Au of the same size, in the range of 150–250 nm. This is consistent with previous calculations for Al nanorods and spheres of comparable size.[56] Other calculations for Al octahedra have shown the presence of two peaks in the octahedron spectrum,[57] one of which in the 200–400 nm range shown here. Note that for Al the scattering and absorption profiles look different because the scattering and absorption peak intensity ratio changes with wavelength. This feature depends on the dielectric constant of the material and is less prominent as we move to Ag and Au. The ratio also depends on the shape and size of the NP.

Case Study 2: Core–Shell NPs

Silica (SiO2) shells are commonplace in nanoscience because silica is an inert material that helps increase the stability of the NPs while its thickness can be used to control the LSPR characteristics of the core material.[48] The optical properties of Al NPs which develop self-limiting alumina (Al2O3) shells are also attractive as the oxide can passivate and protect the NP.[7,58] Here we demonstrate the capabilities of our approach by calculating the scattering and absorption properties of Au@SiO2 sharp and Marks decahedra, Ag@SiO2 cubes, and Al@Al2O3 bipyramids for varying oxide thickness t. The core geometry has an edge length of 50 nm for cubes, sharp decahedra, and bipyramids and 20 nm for Marks decahedra and is as illustrated for all shapes in Figure a.

The scattering and absorption of Au@SiO2 sharp decahedra are shown in Figure a. As the oxide shell increases the plasmon peak red-shifts gradually by 50 nm at t = 10 nm because of the higher RI of SiO2 compared to that of vacuum. This value fits between the ∼20 nm redshifts reported for silica-coated Au spheres of various sizes[59] and larger ∼100 nm shifts for silica-coated Au triangles.[60] This trend is consistent with the observation that higher wavelength LSPRs, found in more anisotropic shapes, exhibit higher RI sensitivity.[16] The scattering and absorption for silica-coated Marks decahedra (Figure b) follow a similar pattern.

Figure 3

Simulated scattering and absorption cross sections for core–shell Au@SiO2 sharp decahedra (a), Au@SiO2 Marks decahedra (b), Ag@SiO2 cubes (c), and Al@Al2O3 bipyramids (d) with varying oxide thickness. Spectra offset for clarity.

Silica-coated Ag cubes (Figure c) demonstrate a slightly higher red-shift than Au for the high wavelength peaks and a smaller red-shift for lower wavelength peaks. The high wavelength LSPR shifts more than what has been reported for SiO2-coated Ag spheres owing to the larger RI sensitivity of cubes.[61]

Figure d shows the scattering and absorption of Al@Al2O3 core–shell bipyramids with different Al2O3 shell thickness. The significant decrease in intensity for an oxide layer as small as t = 10 nm is consistent with reports of similar LSPR suppression for other shapes.[56,58] We also note that Al2O3 causes a notable red-shift of the LSPR positions, following the trend of computational findings for cylinders which, similarly to bipyramids, are highly anisotropic shapes.[56]

Case Study 3: Triangular Plates

Thin nanoplates, including triangles,[62−64] disks,[65] and hexagons,[66] are another interesting group of plasmonic NPs. Triangular plates are quite attractive as their high RI sensitivity, stemming from their sharp corners, makes them suitable for sensing applications[64] while their high anisotropy creates strong local fields.[37] In this last case study, we demonstrate the applicability of our approach by modeling Au and Ag triangular plates and calculating their scattering, absorption, and near-field properties. The modeled NPs (Figure a) have an edge length of 75 nm and thickness of 10 nm.

Figure 4

Wulff modeled triangular plate (a) and calculated (b) E⃗2 field distribution at the NP mid-height and scattering and absorption cross sections for (c) Au and (d) Ag of edge length 75 nm and thickness 10 nm.

Figure b shows the E⃗2 field distribution for a Au triangular plate, calculated at a peak wavelength of 656 nm and shown at the NP mid-height. Here, light is polarized vertically, that is, along the height of the triangle including the top corner. The plasmon-enhanced field is localized at the three sharp corners of the plate, with a higher intensity around the top corner owing to the polarization.[37] Additional, weaker enhancement is present along the NP’s edges. Figure c,d shows the scattering and absorption cross sections for the Au and Ag triangular plates, respectively. Au has a dominant peak at 656 nm with a shoulder at 740 nm, while Ag features at least six distinct peaks ranging from 460 to 630 nm. As expected, for both compositions peaks are red-shifted with respect to the previous studied shapes following the anisotropy trend discussed in the first case study.

Discussion

The results above have shown the applicability of the code to a variety of NPs, some well-studied and some novel; all being crystallographically correct. To our knowledge, current shape-generating tools, used to provide the geometry for various electromagnetic simulation open source codes or commercial packages do not take into account crystallographic directions. For example, built-in 3D drawing platforms are used to generate arbitrary geometry inputs for FEM calculations performed with the COMSOL package and for FDTD packages such as CST Microwave or Lumerical. Some crystallographically correct shapes can be imported as 3D CAD structures created with the open-source FORTRAN code SOWOS[67] which performs Wulff construction modeling, although without including kinetic enhancements or twin planes. Other FDTD software like the open source MEEP[68] use manually created geometries defined through a variety of options including vectors or equations. For the DDA it is possible to find tools, such as DDSCAT Convert,[69] that create the dipole array from a file containing three-dimensional information about the geometry, but again crystallographic orientation is ignored.

Unlike these approaches, with the Wulff construction GUI we provide a facile way to create intrinsically correct NP shapes with correct angles between the NP facets and a consistent area for the facets of the same type. This is encoded in the crystallographic directions that are considered for the Wulff construction, that is, the facet angles are the angles formed between the well-defined crystallographic planes. Currently only the three most stable fcc facets, {111}, {110}, and {100}[70] are considered but the addition of more facets, when required, is trivial. Note that other Wulff shape modelling tools[67,71] are either restricted to shape visualization or the output data needs to be processed in a nontrivial way before used in any electromagnetic simulation software. They also do not offer the required control over the critical parameters for DDA such as the interdipole distance. Conveniently for convergence studies, our tool can be used to systematically alter the interdipole distance, by controlling its value through the appropriate action button, or to maintain the same number of dipoles for varying shape sizes, by fixing the step size value, the latter being applicable only for particles with the same shape.

Further, with the second case study we have specifically demonstrated the use of the Wulff construction tool to calculate the properties of coated NPs. This is a useful feature when investigating optical trends for deliberate or spontaneously formed conformal shells. It therefore, applies very conveniently to most oxide layers, which are of increasing interest given the increasing importance of non-Au plasmonics.

It is important to note that the Wulff construction tool tends to create slightly rounded shapes because of the adopted discretization process; this conveniently happens to mirror the typical experimental shapes. Yet this is not perfect, and one must be careful to recognize that small shape changes can influence the plasmonic behavior and thus care should be taken to choose an appropriate step size for each shape in order to eliminate the so-called staircase effect. As the number of dipoles decreases, deviations from the user-set dimension and interdipole distance are more apparent because the numerical volume calculation becomes inevitably less accurate. The shell is conformal, as appropriate for relatively thin shell layers on shapes of varying complexity.[72] The GUI does not currently support core and shell of different shapes[73] or allow for thickness variations at the NP edges and tips.[60] Another limitation is the minimum thickness that can be modeled: as the code (and DDSCAT) uses the same dipole density for the core and shell, creating a very thin (<4 nm) shell requires a dense array that leads to long computational time.

As a newly developed tool, the Wulff construction GUI has a generous range of future developments including for instance further crystal structures, addition of a substrate,[74] and the ability to model NPs with multiple shells.[20] While the current output is tailored to the DDSCAT, the isosurface can easily be used to generate a CAD shape and in principle provide a crystallographically accurate geometry input for the FDTD and FEM techniques, or for 3D printing.

Conclusions

We described a MATLAB-based standalone GUI that models the shape of fcc NPs, based on the modified kinetic Wulff construction theory, and creates the required input files for the DDSCAT simulations. The range of accessible shapes includes, but is not limited to, cubes, octahedra, bipyramids, stars, plates, pentagonal rods, and multiple decahedra-related structures. All structures modeled have crystallographically correct angles. To demonstrate the capabilities of the GUI we modeled the plasmonic properties of Au, Ag, and Al NPs of various shapes. Next, the effects of oxide shells, including SiO2 and the native Al2O3 on Al, on the optical response of NPs was used to display the capability of the GUI to add a conformal shell on a complex NP. Finally, the near-field and far-field optical properties of triangular plates were also calculated. The results, consistent with the literature when available, show the simplicity and power of the coupled use of the GUI and DDSCAT to predict the plasmonic response of metallic nanomaterials. This GUI is therefore, expected to be an advantageous tool for facilitating the studies of nanoplasmonics, with interesting future extensions.