Abstract
Particle simulation has been widely used in studying plasmas. The technique follows the motion of a large assembly of charged particles in their self-consistent electric and magnetic fields. Plasmons, collective oscillations of the free electrons in conducting media such as metals, are connected to plasmas by very similar physics, in particular, the notion of collective charge oscillations. In many cases of interest, plasmons are theoretically characterized by solving the classical Maxwell’s equations, where the electromagnetic responses can be described by bulk permittivity. That approach pays more attention to fields rather than motion of electrons. In this work, however, we apply the particle simulation method to model the kinetics of plasmons, by updating both particle position and momentum (Newton–Lorentz equation) and electromagnetic fields (Ampere and Faraday laws) that are connected by current. Particle simulation of plasmons can offer insights and information that supplement those gained by traditional experimental and theoretical approaches. Specifically, we present two case studies to show its capabilities of modeling single-electron excitation of plasmons, tracing instantaneous movements of electrons to elucidate the physical dynamics of plasmons, and revealing electron spill-out effects of ultrasmall nanoparticles approaching the quantum limit. These preliminary demonstrations open the door to realistic particle simulations of plasmons.
1 Introduction
Plasmons, collective oscillations of a free electron gas, are usually theoretically investigated by electromagnetic field simulations through solving Maxwell’s equations. The motions of electrons are indirectly represented by the bulk permittivity of metal. Such treatment overlooks real-time motions of free electrons, especially those electrons in the vicinity of the metal surface. For example, as the size of metal nanostructures shrinks below 10 nm, the motions of the surface electrons start to play an essential role, leading to electron spill-out and quantum tunnelling effects [1], [2]. Electromagnetic field simulations often introduce a fictitious layer with modified permittivity to represent these surface electrons for various effects, such as non-local responses [3], [4] and charge transfer plasmons [5], [6]. Although computationally efficient, these treatments might oversimplify the underlying physics and sometimes rely on fitting parameters to describe certain effects. Electromagnetic field simulations also encounter difficulties in modelling the dynamic interaction between a single-electron source and the sample in a typical electron energy loss spectroscopy (EELS) experiment [7]. For example, the electron source is usually modeled as an infinitely long line-like electron beam with zero radius using the boundary element method (BEM) [8], or an infinitely long broadband current density with a radius of 1 nm using the finite element method (FEM) [9], which are both in the frequency domain and do not capture the dynamics when a single electron is passing by the sample. Even if the finite-difference time-domain (FDTD) method is applied, the electron source is approximated by a series of optical dipoles [7], which might not encapsulate all the physics of the electrons in motion at half the speed of light. A microscopic computational tool that traces the freely moving electrons with relatively low computational cost is thus desired.
In fact, plasmons share similar physics and equations of motion as plasmas, one of the four fundamental states of matter. A plasma consists of a gas of charged ions and free electrons with motion–dynamics that are easily affected by the external electromagnetic fields. Among many modelling methods, particle simulation, or the particle-in-cell (PIC) method [10], [11], has become the most powerful simulation tool in plasma community since its development more than half a century ago [12], [13]. The PIC method is a kinetic simulation method, which tracks the motion of a large assembly of charged particles in their self-consistent electric and magnetic fields, and therefore works particularly well for ultrafast (fs-scale) and ultrasmall (nm-scale) problems. It introduces the concept of a macroparticle [10] (also referred to as superparticle, or finite-size particle), a computational particle that represents many real electrons or ions. This key feature makes particle simulation efficient and successful in plasma studies, especially for practical plasma systems with many degrees of freedom, for example, experimental astrophysics [14], plasma-based accelerators [15], [16], [17], inertial confinement fusion [18], attosecond pulse generation [19], and THz radiation sources [20], [21], [22], [23].
Inspired by plasma modelling that properly tracks the motions of free electrons, we start our journey with the application of particle simulation to model plasmons. After carefully examining the essential differences between field and particle simulations, we identify two classic plasmon problems to demonstrate the applicability of particle simulation. In the first example, we simulate single-electron excitation of plasmons on a gold (Au) nanoribbon in a typical EELS experimental setup [9]; and in particular, we analyze the electron dynamics. In the second example, we investigate the electron spill-out effects in sodium (Na) nanowires [24] by making use of the soft boundary condition of particle simulation. Finally, we list the limitations and possible extensions of our current model. Going further, we are interested in applying particle simulation to solve problems such as: how many free electrons indeed participate in the plasmon oscillation, whether surface electrons escape from the metal surface during a plasmon oscillation, i. e., how plasmonic hot-electrons are generated, nonlinear behavior, as well as stochastic and transport phenomena in plasmonics.
2 Particle-in-cell (PIC) simulation
Particle simulation of plasmons is substantially different from classical electromagnetic field simulation methods, as illustrated in Figure 1. The link between them can be established by the free electron model (or plasma model) [25] that was developed in 1927 by Arnold Sommerfeld, where a gas of free electrons of number density n moves against a background of positive ion cores. The details of the lattice potential and electron interactions are not taken into account. Instead, one simply assumes that some aspects of the band structure are incorporated into the effective optical mass m of each electron. The electrons oscillate in response to the applied electromagnetic field, and their motion is damped via collisions occurring with a characteristic collision frequency γ. The dielectric function of the free electron gas is therefore described by [25]:
where

Schematics of field and particle simulation methods. (A) In field simulation of plasmons, the metal is described by dielectric function ε. (B) In particle simulation of plasmons, the metal is represented by number density n of free electrons.
Particle simulation (or the PIC method) simultaneously solves both the Maxwell’s equations:
and the Newton–Lorentz equation:
Here,
Every particle moves and constitutes a local current density
Here, W (x) represents a shape function of the macroparticles [28], which are similarly used in the calculation of
Compared to field simulation, the extra component accessible with particle simulation is the ability to trace the particles, which leads to the second key feature, i. e., soft boundary conditions and auto-inclusion of non-local effects. As shown in Figure 1B, the electromagnetic fields are defined at fixed mesh nodes in PIC, similar to the conventional FDTD method, leading to a hard-wall boundary for the solutions of the fields (represented by the cloud shape around the medium). However, the positions of the particles are not necessarily at the mesh nodes, and the particles are allowed to move anywhere including escaping from the boundary around the medium as shown in Figure 1B. In other words, particle simulation has an intrinsic soft boundary condition for particle motion, which allows us to explore the electron spill-out in quantum plasmonics [2], [24] or other optical charge transport problems [31], [32], [33]. On the other hand, the motion of the electrons is related not only to the field applied at a local position, but also depends on the fields at other positions. This defines the non-local effect [34], [35]. By self-consistently solving the fields and tracing the particles, our PIC method automatically takes into consideration the non-local effects.
3 Results and discussions
In this section, we will elaborate two case studies by particle simulations. The first example, single-electron excitation of plasmons in EELS, will evidently illustrate free electron tracing in particle simulation. It enables us not only to numerically excite plasmons by a single fast electron, but also to efficiently capture the dynamics starting from the very first oscillating event to the steady-state plasmon oscillation. In the second example, leveraging on the implementation of the soft boundary condition, we revisit the famous problem of studying electron spill-out effects of quantum-sized nanoparticles. The results will be benchmarked to hydrodynamic models.
3.1 Single-electron excitation in EELS
In this first case study, we focus on the excitation of plasmons in EELS experiments, where a sample is exposed to a beam of electrons with a known, narrow range of kinetic energies. To model the excitation process, the beam of electrons is often assumed as a current source term in conventional electromagnetic field simulations [8], [9]. Microscopically, this ignores the kinetic nature of the focused electron beam. In scanning transmission electron microscopy (STEM), only < 1% of the electrons emitted by the electron gun are eventually focused onto the sample, with most of them being filtered out to form a monochromatic and focused electron beam. We could optimistically assume that the electron beam generates a current of 0.1 nA on the sample (note: more realistic values could be only a tenth of this number), which is equivalent to 1.5 × 108 electrons per second. Considering a typical lifetime of 100 fs for a surface plasmon resonance [36], only 1.5 × 10−5 electrons on average will fly past the sample before the plasmons dampen-out entirely. We will need an unattainable electron flux 105 times higher to get more than one electron on average to arrive at the sample within the plasmon decay time. Therefore, in an EELS experiment, we almost exclusively generate independent single-electron excitation events that occur sequentially. Here, with the PIC method, we are able to realistically simulate single-electron excitation of plasmons.
We study both bright and dark plasmon modes in an Au nanoribbon with dimensions of 4 × 840 × 20 nm by using single-electron excitation in a setup similar to our previous EELS experiment [9]. As schematically depicted in Figure 2A, dark and bright plasmon modes have zero and nonzero net dipole moments, respectively. Owing to their vanishing dipole moment, dark modes have longer lifetimes and also suffer from suppressed radiative losses [37], [38], making them more efficient for storing electromagnetic energy than bright modes. These attractive properties have made them ideal for applications such as nano-scale sensing [39], and lossless nano-scale waveguides [40]. While bright modes can directly be excited by illuminating the nanoribbon with linearly polarized plane waves (as we later will show), dark plasmonic modes require more intricate excitation schemes [37] such as spatially nonuniform fields [41], evanescent excitation [42], oblique illumination [43], or spatial phase shaping [44].

Particle simulation of plasmons in an Au nanoribbon. (A) Schematic illustrations of dark and bright plasmon modes excited on a nanoribbon. Dark and bright plasmon modes have zero and nonzero net dipole moment respectively across the nanoribbon. (B) Calculated spectral intensity of magnetic fields Bx around the nanoribbon for two different excitation scenarios (schematically shown in inset): edge (indicated by a, red) and center (indicated by d, blue). (C) Various plasmon modes (electric fields Ex in y–z plane) extracted at different resonant energies from spectra in (B), where the region inside the metal has been blacked out for clarity of the surface plasmon modes. In this particle simulation, the Au nanoribbon is modelled using free electron density n = 5.9 × 1028 m−3 with Maxwell-Boltzmann statistics of initial mean electron energy of 0.0375 eV. The nanoribbon is excited by a single macroparticle of charge −500e travelling at a constant velocity of vx = c/2, where c is the speed of light in free space, 1 nm away from the edge of the nanoribbon.
In our PIC simulation of the EELS setup, the electrons within the Au nanoribbon were initialized with a Maxwell-Boltzmann energy distribution, with mean energy of
Figure 2B shows the computed spectra of the plasmon magnetic field Bx, calculated by taking the Fourier transform of the time-domain intensity profile over the entire simulation duration of 100 fs with a time interval of 0.15 fs for excitation locations a (red) and d (blue) using our PIC simulation. The excitation locations and geometry of the nanoribbon are depicted in the inset panel. As expected, we observe that excitation at the edge (location a) results in the formation of both bright and dark modes (indicated by yellow and gray strips respectively) due to the asymmetry of the excitation location with respect to the nanoribbon. Conversely, the excitation at the center (location d) only results in dark modes being formed due to the inherently symmetrical excitation geometry. By plotting the spatial intensity of the plasmon modes, we are able to identify whether each peak corresponds to a dark or bright mode by checking whether the number of high-intensity regions across the Au nanoribbon is odd (dark) or even (bright). The color maps in Figure 2C show the intensity as a function of space for the dark (left panels, location d) and bright (right panels, location a) plasmon modes corresponding to several of the spectral peaks shown in Figure 2B. The cyan arrows in Figure 2C indicate the location of excitation.
To further analyze the electron dynamics within the Au nanoribbon during the plasmonic oscillations, we irradiate the Au nanoribbon with a continuous electromagnetic plane wave centered at 0.78 eV excitation energy (shown in Figure 3A), which corresponds to one of the bright plasmon modes (see Figure 2C). The resulting plasmon field Ex, as in Figure 3B, shows the gradual build up of the mode in about 15 fs (i. e., from 73 to 88 fs), until a steady-state oscillation is reached. Here, the time to establish the mode is about 3 times the plasmon mode cycle.

Diagnostics of electron dynamics in an Au nanoribbon. Time evolution of (A) driving field Ey with 0.78 eV central laser energy and 5 × 108 V/m peak amplitude, and (B) the resulting plasmon field Ex. We plot (C) the spatial distribution and (D) the normalized vertical velocities vy/c of free electrons in the Au nanoribbon at five time points (i)–(v) within a single plasmon cycle of 5.3 fs once steady-state has been attained. The red (blue) dots represent electrons moving upwards (downwards). Unless otherwise stated, the electrons within the Au nanoribbon were initialized with the same parameters as in Figure 2. The nanoribbon is excited by a continuous plane wave linearly polarized along y and propagating in the x direction.
One benefit of PIC simulations of plasmons is the ability to track individual particles; this readily enables the study of electron dynamics of plasmons. We show this in Figure 3C (scatter-plot) and 3D (phase-plot), which depict the steady-state dynamics of the electrons in the Au nanoribbon at five typical times within a single plasmon cycle (cyan dashes in Figure 3A; white dashes in Figure 3B). The temporary local electron densities (moving upwards and downwards) vividly exhibit that electron bunching occurs to establish the plasmon mode, at a frequency of every half cycle. This bunching is strongest while the plasmon field is weakest, as seen at two instances (ii) and (iv) in the figure, suggesting a quarter-cycle time-delay between electron bunching and the plasmon field. It should be noted that the plasmon field data plotted in Figure 3B was collected 5 nm away from the edge of the nanoribbon (at the side of x < 0). For visibility, only electrons with kinetic energy greater than four times the mean value were plotted in Figure 3C, D. The electrons in the Au nanoribbon were initialized with the same parameters as in Figure 2.
The results presented in this section indicate that with regards to plasmon simulation, the PIC method enables us to not only to reproduce results of conventional field simulations, but also allows us a direct view of the actual electron oscillations which underlie the formation of plasmons. In addition to modelling single-electron excitations in EELS experiments, particle simulation could be potentially used to model other electronic excitations, for example, tunneling-electrons excitation of plasmons [45].
3.2 Revealing electron spill-out effects
In the second case study, we leverage on the soft boundary conditions of particle simulation to study the spill-out effect of surface electrons in mesoscopic-scale metal nanoparticles. The term “mesoscopic” refers to deep nanoscale, which lies between the granular microscopic (electronic-scale) and continuous macroscopic (wavelength-scale) regimes [46]. Solving mesoscopic electromagnetic problems has long relied on the well-developed theoretical approaches in either macroscopic or microscopic regimes. However, electromagnetic field-simulation based macroscopic approaches often deal with individual aspects of the omissions, such as nonlocality [34] and free-electron spill-out [24]. Microscopic first-principles approaches, e. g., time-dependent density functional theory [47], are severely constrained by computational demands, thus impractical for multiscale problems. Recently, a general and unified theoretical framework has been developed by introducing a set of mesoscopic boundary conditions to the conventional macroscopic Maxwell’s equations [46]. Motivated by the same objective, our PIC method is an alternative approach starting from a microscopic viewpoint, but is much more computationally efficient than other first-principles approaches by emphasizing only the motions of free electrons.
In an attempt to use a plasma model to simulate plasmons in a solid-state material, there are a few things to note. In plasma simulations, the electrons are separated from the ions and totally free in space, and they are considered all equivalent following classical Maxwell-Boltzmann statistics. But in solid-state physics, the electrons follow Fermi–Dirac statistics, in which each quantum state can only be occupied by a single electron. More importantly, these Fermi electrons are only free to move within the material, as they are still loosely bound to nuclei, therefore energy is required for them to escape from the material. In PIC simulation, an external potential-wall can be set up at the material boundary to represent such an energy barrier. Depending on the height of this potential-wall, the electrons can be partially or totally confined within the material. With the freedom to set the potential-wall or leave it open, we designate such boundary as a soft boundary.
As an example, we apply the PIC method to study the electron spill-out in a sodium (Na) nanowire, which is commonly used to differentiate the standard hydrodynamic model with a hard-wall boundary condition (HW-HDM) [48] and the self-consistent hydrodynamic method (SC-HDM) [24]. The major difference between these two models is that the electrons in HW-HDM are strictly confined in the metallic structure, while SC-HDM generalizes standard hydrodynamic theory to include free-electron spill-out at metal–dielectric interfaces. As a result, SC-HDM predicts red-shifted resonances of infinitely long Na nanowires with 1/R scaling (R: nanowire radius), whereas HW-HDM incorrectly predicts a blue-shift [24]. To elucidate the effect of electron spill-out, we implement either a potential-wall or an open boundary condition in the same PIC framework, corresponding to HW-HDM and SC-HDM, respectively. The potential wall is created by applying an external static electric-field of 117 V/nm over a 2 nm shell surrounding the nanowire. This number is chosen merely to represent a very large potential wall, which is equivalent to 100 times the work function of Na over 2 nm thickness.
As evidently indicated in Figure 4A, on the same basis of Fermi–Dirac statistics, the open boundary PIC simulations predict a red-shift on the 1/R scaling (red line), whereas the potential-wall PIC simulations have an opposite blue-shift prediction (black line). Clearly, it is the spill-out of electrons that has caused the red-shift, but only for Na nanowires of radius R smaller than a critical value
![Figure 4: Revealing electron spill-out effects. (A) Simulated surface plasmon resonant energies of infinitely long Na nanowires as a function of the inverse of their radius (1/R), employing either open or potential-wall boundary conditions in the PIC method, benchmarked to the self-consistent hydrodynamic method (SC-HDM) [24]. (Inset) The Na nanowires are assumed infinitely long in z-direction and 2D PIC simulations were performed. (B) For Na nanowire with radius R = 2 nm excited at its plasmon resonance, the steady-state space distributions of the free electrons density and induced charge density (inset) near the edge of the nanowire, for open (red) and potential-wall (black) boundary conditions. (C) Snapshots of the free electron movements near the upper edge of the Na nanowire with radius R = 2 nm to show the soft boundary in PIC, which can be set as potential-wall (upper panel) or open (lower panel) to restrict or allow surface electrons to spill out. The red (blue) dots represent electrons moving upwards (downwards). In this 2D particle simulation, the Na nanowire is initialized using free electron density n = 2.5 × 1028 m−3 with either Maxwell–Boltzmann of initial mean electron energy of 0.0375 eV or Fermi-dirac statistics with EF = 3.24 eV at T = 300 K. The nanowire is excited by a single negatively charged macroparticle travelling at a constant velocity of vx = c/2, where c is the speed of light in free space, and 0.1 nm away from the edge of the nanowire in (A). To illustrate electron dynamics at plasmon resonance, the Na nanowire is excited by a continuous plane wave with central laser frequency equal to the plasmon resonant frequency, linearly polarized along y and propagating in the x direction in (B) and (C).](/document/doi/10.1515/nanoph-2020-0067/asset/graphic/j_nanoph-2020-0067_fig_004.png)
Revealing electron spill-out effects. (A) Simulated surface plasmon resonant energies of infinitely long Na nanowires as a function of the inverse of their radius (1/R), employing either open or potential-wall boundary conditions in the PIC method, benchmarked to the self-consistent hydrodynamic method (SC-HDM) [24]. (Inset) The Na nanowires are assumed infinitely long in z-direction and 2D PIC simulations were performed. (B) For Na nanowire with radius R = 2 nm excited at its plasmon resonance, the steady-state space distributions of the free electrons density and induced charge density (inset) near the edge of the nanowire, for open (red) and potential-wall (black) boundary conditions. (C) Snapshots of the free electron movements near the upper edge of the Na nanowire with radius R = 2 nm to show the soft boundary in PIC, which can be set as potential-wall (upper panel) or open (lower panel) to restrict or allow surface electrons to spill out. The red (blue) dots represent electrons moving upwards (downwards). In this 2D particle simulation, the Na nanowire is initialized using free electron density n = 2.5 × 1028 m−3 with either Maxwell–Boltzmann of initial mean electron energy of 0.0375 eV or Fermi-dirac statistics with EF = 3.24 eV at T = 300 K. The nanowire is excited by a single negatively charged macroparticle travelling at a constant velocity of vx = c/2, where c is the speed of light in free space, and 0.1 nm away from the edge of the nanowire in (A). To illustrate electron dynamics at plasmon resonance, the Na nanowire is excited by a continuous plane wave with central laser frequency equal to the plasmon resonant frequency, linearly polarized along y and propagating in the x direction in (B) and (C).
Besides predicting the resonant surface plasmon energies, we are also able to compute directly from the PIC model the space distribution of free electron densities, similar to other microscopic theoretical approaches. Taking a particular Na nanowire with R = 2 nm as an example, in Figure 4B, we plot the steady-state free electron density around the edge of the nanowire, with open and potential-wall boundary conditions, after 10 fs from the start of the simulation (note: this steady state is self-evolved and developed from an initialized uniform distribution of electrons). From Figure 4B, the spill-out of the free electrons is more prominently present in the simulation with open boundary. Nevertheless, it is questionable that quite a significant number of electrons are still observable outside the material boundary for the potential-wall boundary condition. This is due to the interpolation algorithm in PIC simulation. To better distinguish the electron spill-out for the two boundary conditions, we also plot the induced charge density (i. e., the difference between ion and electron density) [46] in the inset of Figure 4B. For open boundary (red curve), “+” (or “−”) charge accumulation appears just inside (or outside) the edge of the Na nanowire, indicating the electron spill-out. In contrast, “−” charge accumulation occurs only within the nanowire for the potential-wall boundary (black curve), representing electron spill-in.
More vividly, for the same nanowire, we show the real-time electron movements in the x–y plane near the upper edge of the nanowire at four typical times within a single plasmon cycle T0 in Figure 4C. With red (blue) dots representing electrons moving upwards (downwards), it is pictorially evident that electron bunching occurs periodically, forming the plasmon dipole mode (see t = 0 with most red dots and t = 0.5T0 with most blue dots), for both potential-wall and open boundary conditions. However, the potential-wall boundary restricts the spill-out of the electrons in the vicinity of the Na surface, whereas the open boundary indeed allows some electrons to get out, because their movements are purely determined by the local electromagnetic fields.
The results presented in this section give a demonstration of an important feature of particle simulation, i. e., the soft boundary condition, where we have the flexibility to set the boundary to be open or a potential-wall. This allows us to capture the microscopic dynamics of electrons at the material interface (in this example, the electrons are spilling out of Na into the vacuum) so as to study the effects of such dynamics. We could further modify the settings of our potential-wall boundary condition to include more physics or even some quantum effects. By doing so, the potential of particle simulation is unlimited, e. g., an extension to Ag nanowires by including inter-band transitions (see below), bimetallic nanoplasmonic structures, electron spill-out and charge transport between nearly touching nanoparticles.
4 Limitations and extensions
4.1 Inter-band transitions
Over a wide frequency range, the optical properties of metals can be explained by the free electron model [25]. To be more precise, the free-electron description adequately describes the optical response of metals only for photon energies below the threshold of transitions between electron bands. It breaks down when inter-band transitions occur. For alkali metals, this range extends up to the ultraviolet, while for noble metals inter-band transitions occur at visible and near-infrared frequencies, limiting the validity of this approach. Above their respective band edge thresholds, photons are very efficient in inducing inter-band transitions, where electrons from the filled band below the Fermi surface are excited to higher bands. The main consequence of these processes concerning plasmons is an increased damping and competition between the two excitations at visible frequencies. At the moment, the contribution from valence electrons is neglected in PIC simulations, while the ions (composed of the nucleus, core and valence electrons) only keep the charge neutrality in the metal. To take into account the screening effect of the valence electrons governed by the inter-band transitions, another species of particles with larger effective mass can be introduced into our PIC simulation to represent valence electrons. To simplify, as these heavier valence electrons are moving closely around the nucleus, they effectively cause a highly polarized environment, which can be modelled by a background
4.2 Damping via collisions
The electrons oscillate in response to the applied electromagnetic field, and their motion is damped via collisions occurring with a characteristic collision frequency
4.3 Fermi-electrons in metals
To model solid metals using plasma-physics based PIC, that is normally used to model quantum plasmas [49], the quantum effects become important when the temperature T is lower than the Fermi temperature TF, hence the relevant statistical distribution changes from Maxwell–Boltzmann to Fermi–Dirac. Besides, energy (i. e., work function) is required to free the Fermi-electrons from the metal. It is possible to implement such an energy barrier by setting up a potential-wall at the surface of the metal, as illustrated in our electron spill-out study. So far, we have only tried a very basic implementation of a sharp potential-wall by adding an external static electric-field over an extremely short distance. This electric field will simply push away those electrons approaching the wall. To further improve the potential-wall implementation, we could consider a more accurate profile of the potential barrier and the quantum tunneling probability for each electron [6], [31]. Depending on their kinetic energies, some electrons will be reflected back to the metal, while some will be allowed to transmit through the wall. This idea of a quantum potential-wall could be implemented in future.
4.4 Shape function
Non-locality in plasmonic structures is the result of electron-electron repulsion and consequently a wave-vector dependent response of the system. In PIC simulation, both factors are implicitly taken into account by charge deposition and field interpolation at mesh nodes through the shape function defined in Eq. (4) [10]. This opens up an opportunity to use shape function as a fitting parameter for PIC simulation of non-local effects in plasmons. For example, different orders of shape function could be implemented to account for the range of electron-electron repulsion, which could affect the slope of the size-dependent curve in Figure 4A in a similar manner as the β factor of the hydrodynamic model [34]. The shape function used in our current work is based on second order interpolation, where the function extends over a distance equal to three times the PIC cell size [30], which is set to be 0.1 nm. Higher orders will be more accurate but computationally more expensive [28]. This topic can be further explored in future work to thoroughly interpret the PIC results, together with benchmarked data such as experiments or SC-HDM results.
5 Conclusions and outlook
We have demonstrated a new theoretical framework to model plasmons – particle simulation, which not only solves for electromagnetic fields, but also traces the motions of electrons. In doing so, our numerical experiment by particle simulation is able to provide a microscopic (electronic-scale) view of plasmons. In the first example of studying plasmons inside an Au nanoribbon, we have employed a single-electron source, mimicking realistic EELS experiments, to efficiently excite various plasmon modes in one go. Different excitation scenarios can be handled straightforwardly by our particle simulations. To analyze the electron dynamics for a particular plasmon mode, we have the flexibility to switch to plane-wave excitation at a central frequency and extract information of interest. For example, we have vividly illustrated the spatial bunching of free electrons with their corresponding velocities during a cycle of a plasmon oscillation, noticing that the plasmon field is a quarter-cycle delayed with respect to the electron bunching. Other diagnosis techniques could be applied to reveal different aspects of plasmons.
The other uniqueness of particle simulation lies in the soft boundary conditions. We have made use of this property to demonstrate the possibility of modelling electron spill-out effects in Na nanowires, benchmarked to the well-developed hydrodynamic models. It should be emphasized that particle simulation is able to differentiate between free and valence electrons, and can turn on and off an external potential-wall at the material boundary. This operational flexibility allows us to differentiate alkali and noble metals, and to identify regimes where electron spill-out or retardation effects dominate. Moving forward, our PIC model could be further developed by taking into consideration the inter-band transitions, the electron-electron collisions, and a more comprehensive description of Fermi-electrons in metals and beyond. This is still a largely unexplored field, and is thus likely to become an active area of future research.
Funding source: Science and Engineering Research Council
Award Identifier / Grant number: Young Individual Research Grants (YIRG No. A1784c0
Funding source: National Research Foundation Singapore
Award Identifier / Grant number: NRF2016-NRF-ANR002NRF2017-NRFNSFC002-015
Funding source: MOE Singapore
Award Identifier / Grant number: T2 grant (2018-T2-1-007) and MOE PhD RSSTier 1 (R-284-000-179-133).
Funding source: Agency for Science, Technology and Research
Funding source: National Research Foundation Singapore
Author contribution: All the authors have accepted responsibility for the entire content of this submitted manuscript and approved submission.
Employment or leadership: None declared.
Honorarium: None declared.
Conflict of interest statement: The authors declare no conflicts of interest regarding this article.
Acknowledgments
The Institute of High Performance Computing (IHPC) acknowledges financial support from A*STAR SERC Young Individual Research Grants (YIRG No. A1784c0020), the National Research Foundation Singapore (NRF2017-NRF-NSFC002-015 and NRF2016-NRF-ANR002). J. Z. J. L. and L. K. A. would like to acknowledge the support of the MOE T2 grant (2018-T2-1-007), USA ONRG grant (N62909-19-1-2047) and MOE PhD RSS. M. B. acknowledges support from the Singapore Ministry of Educationʼs Academic Research Fund Tier 1 (R-284-000-179-133). The computational work for this article was partially performed on resources of the National Supercomputing Centre, Singapore (https://www.nscc.sg).
References
[1] M. S. Tame, K. R. McEnery, Ş. K. Özdemir, J. Lee, S. A. Maier, and M. S. Kim, “Quantum plasmonics,” Nat. Phys., vol. 9, no. 6, pp. 329–340, 2013, https://doi.org/10.1038/nphys2615.Search in Google Scholar
[2] D. Xu, X. Xiong, L. Wu, et al., “Quantum plasmonics: new opportunity in fundamental and applied photonics,” Adv. Opt. Photonics, vol. 10, no. 4, pp. 703–756, 2018, https://doi.org/10.1364/AOP.10.000703.Search in Google Scholar
[3] J. M. McMahon, S. K. Gray, and G. C. Schatz, “Nonlocal optical response of metal nanostructures with arbitrary shape,” Phys. Rev. Lett., vol. 103, no. 9, 2009, Art no. 97403, https://doi.org/10.1103/PhysRevLett.103.097403.Search in Google Scholar
[4] A. I. Fernández-Domínguez, A. Wiener, F. J. García-Vidal, S. A. Maier, and J. B. Pendry, “Transformation-optics description of nonlocal effects in plasmonic nanostructures,” Phys. Rev. Lett., vol. 108, no. 10, 2012, Art no. 106802, https://doi.org/10.1103/PhysRevLett.108.106802.Search in Google Scholar
[5] R. Esteban, A. G. Borisov, P. Nordlander, and J. Aizpurua, “Bridging quantum and classical plasmonics with a quantum-corrected model,” Nat. Commun., vol. 3, 2012, Art no. 825, https://doi.org/10.1038/ncomms1806.Search in Google Scholar
[6] L. Wu, S. F. Tan, M. Bosman, et al., “Charge transfser plasmon resonances across silver-molecule-silver junctions: estimating the terahertz conductance of molecules at near-infrared frequencies,” RSC Adv., vol. 6, no. 75, pp. 70884–70894, 2016, https://doi.org/10.1039/C6RA16826D.Search in Google Scholar
[7] P. Das, T. Kumar Chini, and P. James, “Probing higher order surface plasmon modes on individual truncated tetrahedral gold nanoparticle using cathodoluminescence imaging and spectroscopy combined with fdtd simulations,” J. Phys. Chem. C, vol. 116, no. 29, pp. 15610–15619, 2012, https://doi.org/10.1021/jp3047533.Search in Google Scholar
[8] F. J. García de Abajo, “Optical excitations in electron microscopy,” Rev. Modern Phys., vol. 82, no. 1, pp. 209–275, 2010, https://doi.org/10.1103/RevModPhys.82.209.Search in Google Scholar
[9] Z. Fan, M. Bosman, X. Huang, et al., “Stabilization of 4H hexagonal phase in gold nanoribbons,” Nat. Commun., vol. 6, 2015, Art no. 7684, https://doi.org/10.1038/ncomms8684.Search in Google Scholar
[10] J. M. Dawson, “Particle simulation of plasmas,” Rev. Modern Phys., vol. 55, no. 2, 1983, Art no. 403, https://doi.org/10.1103/RevModPhys.55.403.10.1103/RevModPhys.55.403Search in Google Scholar
[11] C. K. Birdsall, and A. B. London, Plasma Physics via Computer Simulation, Boca Raton, Taylor and Francis, 2004.Search in Google Scholar
[12] J. Dawson, “One-dimensional plasma model,” Phys. Fluids, 5, no. 4, pp. 445–459, 1962, https://doi.org/10.1063/1.1706638.Search in Google Scholar
[13] A. B. Langdon, and K. B. Charles, “Theory of plasma simulation using finite-size particles,” Phys. Fluids, vol. 13, no. 8, pp. 2115–2122, 1970, https://doi.org/10.1063/1.1693209.Search in Google Scholar
[14] S. Mondal, V. Narayanan, W. J. Ding, et al., “Direct observation of turbulent magnetic fields in hot, dense laser produced plasmas,” Proc. Natl. Acad. Sci. U.S.A., vol. 109, no. 21, pp. 8011–8015, 2012, https://doi.org/10.1073/pnas.1200753109.Search in Google Scholar
[15] E. Esarey, C. B. Schroeder, and W. P. Leemans, “Physics of laser-driven plasma-based electron accelerators” Rev. Modern Phys., vol. 81:1229–1285, 2009, https://doi.org/10.1103/RevModPhys.81.1229.Search in Google Scholar
[16] A. Macchi, M. Borghesi, and M. Passoni, “Ion acceleration by superintense laser-plasma interaction,” Rev. Modern Phys., vol. 85, pp. 751–793, 2013, https://doi.org/10.1103/RevModPhys.85.751.Search in Google Scholar
[17] X. Q. Yan, H. C. Wu, Z. M. Sheng, J. E. Chen, and J. Meyer-ter Vehn, “Self-organizing GeV, nanocoulomb, collimated proton beam from laser foil interaction at 7 × 1021 W/cm2,” Phys. Rev. Lett., vol. 103, 2009, Art no. 135001, https://doi.org/10.1103/PhysRevLett.103.135001.Search in Google Scholar
[18] K. Mima, A. E. Dangor, R. G. Evans, et al., “Fast heating of ultrahigh-density plasma as a step towards laser fusion ignition,” Nature, vol. 412, pp. 798–802, 2001, https://doi.org/10.1038/35090525.Search in Google Scholar
[19] C. Thaury, and F. Quéré, “High-order harmonic and attosecond pulse generation on plasma mirrors: basic mechanisms,” J. Phys. B, vol. 43, 2010, Art no. 213001, https://doi.org/10.1088/0953-4075/43/21/213001.Search in Google Scholar
[20] G. Liao, Y. Li, H. Liu, et al., “Multimillijoule coherent terahertz bursts from picosecond laser-irradiated metal foils,” Proc. Natl. Acad. Sci., vol. 116, no. 10, pp. 3994–3999, 2019, https://doi.org/10.1073/pnas.1815256116.Search in Google Scholar
[21] W. J. Ding, Z. M. Sheng, and W. S. Koh, “High-field half-cycle terahertz radiation from relativistic laser interaction with thin solid targets,” Appl. Phys. Lett., 103, no. 20, 2013, Art no. 204107, https://doi.org/10.1063/1.4831684.Search in Google Scholar
[22] Z.-M. Sheng, K. Mima, J. Zhang, and H. Sanuki, “Emission of electromagnetic pulses from laser wakefields through linear mode conversion,” Phys. Rev. Lett., vol. 94, 2005, Art no. 095003, https://doi.org/10.1103/PhysRevLett.94.095003.Search in Google Scholar
[23] S. Mondal, Q. Wei, W. J. Ding, et al., “Aligned copper nanorod arrays for highly efficient generation of intense ultra-broadband thz pulses,” Sci. Rep., vol. 7, 2017, Art no. 40058, https://doi.org/10.1038/srep40058.Search in Google Scholar
[24] G. Toscano, J. Straubel, K. Alexander, et al., “Resonance shifts and spill-out effects in self-consistent hydrodynamic nanoplasmonics,” Nat. Commun., vol. 6, 2015, Art no. 7132, https://doi.org/10.1038/ncomms8132.Search in Google Scholar
[25] S. A. Maier, Plasmonics: Fundamentals and Applications, New York, NY, Springer Science & Business Media, 2007.10.1007/0-387-37825-1Search in Google Scholar
[26] T. Umeda, Y. Omura, T. Tominaga, and H. Matsumoto, “A new charge conservation method in electromagnetic particle-in-cell simulations,” Comput. Phys. Commun., vol. 156, no. 1, pp. 73–85, 2003, https://doi.org/10.1016/S0010-4655(03)00437-5.Search in Google Scholar
[27] T.Z. Esirkepov, “Exact charge conservation scheme for particle-in-cell simulation with an arbitrary form-factor,” Comput. Phys. Commun., vol. 135, no. 2, pp. 144–153, 2001, https://doi.org/10.1016/S0010-4655(00)00228-9.Search in Google Scholar
[28] E. Cormier-Michel, B. A. Shadwick, C. G. R. Geddes, E. Esarey, C. B. Schroeder, and W. P. Leemans, “Unphysical kinetic effects in particle-in-cell modeling of laser wakefield accelerators,” Phys. Rev. E, vol. 78, 2008, Art no. 016404, https://doi.org/10.1103/PhysRevE.78.016404.Search in Google Scholar
[29] T. D. Arber, K. Bennett, C. S. Brady, et al., “Contemporary particle-in-cell approach to laser-plasma modelling,” Plasma Phys. Control. Fusion, vol. 57, no. 11, 2015, Art no. 113001, https://doi.org/10.1088/0741-3335/57/11/113001.Search in Google Scholar
[30] J. Derouillat, A. Beck, F. Pérez, et al., “Smilei: a collaborative, opensource, multi-purpose particle-in-cell code for plasma simulation,” Comput. Phys. Commun., vol. 222, pp. 351–373, 2018, https://doi.org/10.1016/j.cpc.2017.09.024.Search in Google Scholar
[31] L. Wu, H. Duan, P. Bai, M. Bosman, K. W. Y. Joel, and E. Li, “Fowler-Nordheim tunneling induced charge transfer plasmons between nearly touching nanoparticles,” ACS Nano, vol. 7, no. 1, pp. 707–716, 2012, https://doi.org/10.1021/nn304970v.Search in Google Scholar
[32] Z. Zhu, S. Joshi, S. Grover, and G. Moddel, “Graphene geometric diodes for terahertz rectennas,” J. Phys. D Appl. Phys., vol. 46, no. 18, 2013, Art no. 185101, https://doi.org/10.1088/0022-3727/46/18/185101.Search in Google Scholar
[33] S. F. Tan, L. Wu, J. K. W. Yang, P. Bai, M. Bosman, and C. A. Nijhuis, “Quantum plasmon resonances controlled by molecular tunnel junctions,” Science, vol. 343, no. 6178, pp. 1496–1499, 2014, https://doi.org/10.1126/science.1248797.Search in Google Scholar
[34] S. Raza, I. B. Sergey, M. Wubs, and N. A. Mortensen, “Nonlocal optical response in metallic nanostructures,” J. Phys. Condensed Matter, vol. 27, no. 18, 2015, Art no. 183204, https://doi.org/10.1088/0953-8984/27/18/183204.Search in Google Scholar
[35] J. A. Scholl, A. L. Koh, and J. A. Dionne, “Quantum plasmon resonances of individual metallic nanoparticles,” Nature, vol. 483, no. 7390, pp. 421–427, 2012, https://doi.org/10.1038/nature10904.Search in Google Scholar
[36] M. Bosman, E. Ye, S. F. Tan, et al., “Surface plasmon damping quantified with an electron nanoprobe,” Sci. Rep., vol. 3, 2013, Art no. 1312, https://doi.org/10.1038/srep01312.Search in Google Scholar
[37] N. J. Halas, S. Lal, W.-S. Chang, S. Link, and P. Nordlander, “Plasmons in strongly coupled metallic nanostructures,” Chem. Rev., vol. 111, no. 6, pp. 3913–3961, 2011, https://doi.org/10.1021/cr200061k.Search in Google Scholar
[38] D. E. Gómez, Z. Q. Teo, M. Altissimo, T. J. Davis, S. Earl, and A. Roberts, “The dark side of plasmonics,” Nano Lett., vol. 13, no. 8, pp. 3722–3728, 2013, https://doi.org/10.1021/nl401656e.Search in Google Scholar
[39] S. Zhang, and H. Xu, “Tunable dark plasmons in a metallic nanocube dimer: toward ultimate sensitivity nanoplasmonic sensors,” Nanoscale, vol. 8, no. 28, pp. 13722–13729, 2016, https://doi.org/10.1039/C6NR03806A.Search in Google Scholar
[40] D. SolisJr, B. Willingham, L. N. Scott, et al., “Electromagnetic energy transport in nanoparticle chains via dark plasmon modes,” Nano Lett., vol. 12, no. 3, pp. 1349–1353, 2012, https://doi.org/10.1021/nl2039327.Search in Google Scholar
[41] J.-S. Huang, J. Kern, P. Geisler, et al., “Mode imaging and selection in strongly coupled nanoantennas,” Nano Lett., vol. 10, no. 6, pp. 2105–2110, 2010, https://doi.org/10.1021/nl100614p.Search in Google Scholar
[42] S.-C. Yang, H. Kobori, C.-L. He, et al., “Plasmon hybridization in individual gold nanocrystal dimers: direct observation of bright and dark modes,” Nano Lett., vol. 10, no. 2, pp. 632–637, 2010, https://doi.org/10.1021/nl903693v.Search in Google Scholar
[43] W. Zhou, and T. W. Odom, “Tunable subradiant lattice plasmons by out-ofplane dipolar interactions,” Nat. Nanotechnol., vol. 6, no. 7, pp. 423–427, 2011, https://doi.org/10.1038/nnano.2011.72.Search in Google Scholar
[44] G. Volpe, S. Cherukulappurath, R. J. Parramon, G. Molina-Terriza, and R. Quidant, “Controlling the optical near field of nanoantennas with spatial phase-shaped beams,” Nano Lett., vol. 9, no. 10, pp. 3608–3611, 2009, https://doi.org/10.1021/nl901821s.Search in Google Scholar
[45] W. Du, T. Wang, H.-S. Chu, et al., “On-chip molecular electronic plasmon sources based on self-assembled monolayer tunnel junctions,” Nat. Photonics, vol. 10, no. 4, pp. 274–280, 2016, https://doi.org/10.1038/nphoton.2016.43.Search in Google Scholar
[46] Y. Yang, D. Zhu, W. Yan, et al., “A general theoretical and experimental framework for nanoscale electromagnetism,” Nature, vol. 576:248–252, 2019, https://doi.org/10.1038/s41586-019-1803-1.Search in Google Scholar
[47] J. Zuloaga, E. Prodan, and P. Nordlander, “Quantum description of the plasmon resonances of a nanoparticle dimer,” Nano Lett., vol. 9, no. 2, pp. 887–891, 2009, https://doi.org/10.1021/nl803811g.Search in Google Scholar
[48] V. T. Tatiana, P. Nordlander, J. Aizpurua, and A. G. Borisov, “Quantum effects and nonlocality in strongly coupled plasmonic nanowire dimers,” Opt. Express, vol. 21, no. 22, pp. 27306–27325, 2013, https://doi.org/10.1364/OE.21.027306.Search in Google Scholar
[49] G. Manfredi, “How to model quantum plasmas,” Fields Inst. Commun., vol. 46, pp. 263–287, 2005.10.1090/fic/046/10Search in Google Scholar
© 2020 Wen Jun Ding et al., published by De Gruyter, Berlin/Boston
This work is licensed under the Creative Commons Attribution 4.0 International License.