Simple model of saturable localised surface plasmon

Localised surface plasmons (LSPs) are now applied to various fields, such as bio-sensing, solar cell, molecular fluorescence enhancement and quantum-controlled devices at nanometre scale. Recent experiments show that LSPs are optically saturated by high-intensity light. Absorption saturation arises as a result of strong optical nonlinearity and cannot be explained by the conventional boson model of LSPs. Here, we propose a simple model of saturable LSPs using an effective dipole approximation. The strategy is to directly compare the classical linear optical response of an LSP with that obtained from a saturable quantum two-level system in the limit of weak excitation. The second quantization can then be performed by replacing a classical polarizability with a quantum dipole operator. Taking an ellipsoidal nanometal as an example, we analyse in detail the optical response of a single ellipsoidal nanometal to validate our model. Our numerical results show that the plasmon resonance frequency and spectral linewidth decrease as the aspect ratio of the ellipsoid increases, which is similar to the size dependence observed in early experiments.

Localised surface plasmons (LSPs) are quantized plasma oscillations formed near the surfaces of metal nanoparticles. In contrast to propagating surface plasmons, which often require carefully-constructed optical arrangement for phase matching, LSPs can be excited easily by direct irradiation of light. It is well known that LSPs can focus light to nanometre scale and strongly enhance the electric field near nanometals, which is called the antenna effect. The research field dealing with such plasmon characteristics is called "plasmonics", and has been extensively investigated. In fact, many applications of LSPs, such as bio-sensing 1,2 , solar cell [3][4][5] and molecular fluorescence enhancement [6][7][8] , have been reported.
For several years, the studies focused on the quantum nature of plasmons are now rapidly evolving, in which the light field interacts with matter at nanometre scale and plasmons play a role in controlling the light-matter interactions at the quantum level. This new research field is called "quantum plasmonics 9,10 ". Quantum plasmonics has opened up a new frontier in the study of the fundamental physics of LSPs [11][12][13][14][15] , namely vacuum Rabi splitting and the Fano effect, and now provide new potential applications of plasmons, such as single-photon sources 16,17 , entangled-photon sources 18,19 and quantum-controlled devices at the nanometre scale [20][21][22][23][24] .
As is well known, however, the optical responses of LSPs are generally analysed using classical electromagnetism, in spite of the fact that LSPs are quantized plasma oscillations. This is because that the LSP antenna effect can be understood without the quantum nature of plasmon and can be explained simply by Maxwell's equations. However, quantum plasmonics requires a quantum-mechanical treatment of plasmons, in other words, the second quantization of plasmons 25,26 . In particular, the understanding of the quantum properties of a single LSP, namely dipole moment and its relaxation rate, is required to understand the vacuum Rabi splitting in LSP-matter interaction, because the strong coupling is realized when the coupling constant rate g of LSP-matter interaction is lager than the relaxation rates of the LSP and the matter.
Conventionally, the second quantization of surface plasmons is based on the framework of cavity QED theory and is described by the boson model. In general, the boson model is justified by the fact that the Pauli exclusion principle hardly affects collective excitation. For example, in a lattice model of condensed-matter physics, the prohibition of two excitations at the same lattice corresponds directly to the Pauli exclusion principle. This prohibition component constitutes only N −1 of a collective excitation mode, where N is the number of lattice, and therefore can be ignored for N 1  . This boson approximation is valid for collective excitations of general quasiparticles. Recently, however, it has been experimentally reported that LSPs are optically saturated by high-intensity light 27 . Absorption saturation arises as a result of strong optical nonlinearity and cannot be explained by the above boson model. In contrast to propagating surface plasmons formed at the surface of planar metallic films, LSPs are localised in literally at nanometre scale. Therefore, the Pauli exclusion principle might be non-negligible, especially for small metal nanoparticles.
In this study, we propose a simple model of a saturable LSP for small metal nanoparticles using an effective dipole approximation. Our strategy is to directly compare the classical linear optical response of an LSP with that obtained from a saturable quantum two-level system. The second quantization is then performed by replacing a classical polarizability of the LSP with a quantum dipole operator described by the two-level system. Taking an ellipsoidal nanometal as an example, we introduce an optical response function of an LSP, based on the optical Bloch equations, and validate our method by analysing in detail the size dependence of the plasmon resonance frequency and the relaxation decay rate of a single LSP. Our numerical results show that the plasmon resonance frequency and spectral linewidth decrease as the aspect ratio of the ellipsoidal nanometal increases, which is similar to the size dependence observed in early experiments.

Results
Second quantization of LSPs using effective dipole approximation. We start by considering a small metal ellipsoid with semiaxes a x , a y and a z interacting with incident light E in with a wavelength of λ, as depicted in Fig. 1. Assuming  λ a , we restrict ourselves to dipole radiation from the metal ellipsoid, because the quadrupole and magnetic dipole radiations are both negligibly small. When we focus on the quasi-static approximation 28 and plasmon excitation at low-intensity light, dipole moment p i along the principal axes (i = x, y, z) can be described as a linear response of E in,i , given by where ε m is the dielectric constant of the background. The polarizability along the principal axes, α i , depends on the complex dielectric function ε(ω) of the ellipsoid: In order to quantize the LSP of the ellipsoidal nanometal, we follow the experimental facts 29 : The LSP absorbs a photon with a specific energy ћω R and the LSP spectrum can be approximated by Lorentzian function. First, we derive the Lorentzian function from Eq. (2). By assuming the incident light with a frequency close to ω R , we approximate Eq. (2) by taking the Taylor expansion of ε(ω) around ω R . The plasmon resonance enhancement takes place under the condition that |ε m + L i (ε(ω) − ε m )| is a minimum, which for the case of small or slowly-varying Im[ε(ω)] around ω R can simplify to Using this condition, the complex dielectric function ε(ω) can be approximated as 14 (2), we obtain  , the first term in parenthesis in Eq. (6) can be ignored, and we can rewrite Eq. (1) as the shape of the Lorentzian function, Next, we follow another experimental fact that the LSPs are saturable 27 , and derive the dipole moment quantum-mechanically by considering a saturable quantum dipole operator. The simplest way to treat the saturation is to adopt the dipole operator of a two-level system, defined as p d e g g e ( ) i î = | 〉〈 | + | 〉〈 | , where |g〉 and |e〉 are the ground and excited states of a quantum two-level system and d i is the dipole moment along the i axis. In the limit of weak excitation, the optical response of arbitrary quantum two-level systems can be approximated by where Γ is the dipole relaxation rate and ω 0 is the energy difference between |g〉 and |e〉.
Finally, we compare the classical p i of Eq. (9) are satisfied, equation (9) becomes the same as Eq. (10). The LSP of the ellipsoidal nanometal can thus be reduced to an effective quantum two-level system with plasmon resonance frequency ω R and dipole relaxation rate γ, as depicted in Fig. 2. The Hamiltonian of the LSP is given by | 〉〈 |. It should be noted that, for the case of LSP excitation by low-intensity light, the nonlinearity of the LSP is negligibly small. In this case, the LSP can be approximated by using a bosonic operator cˆ instead of p i in Eq. (10), as can be seen in the conventional quantization methods of surface plasmons.

Plasmonic Bloch equations.
Here, we extend p î 〈 〉 to the case for no limit of weak excitation, presuming that the light polarization is parallel to a semiaxis a i of the ellipsoid. In general, fully quantum-mechanical dynamics of a two-level system can be described by the quantum Langevin equations. Introducing the dipole operator is an input filed describing quantum fluctuation. With the input-output relation, , Eqs (14) and (15) can be directly applied to the analysis of the quantum properties of photons emitted from a single LSP, such as single-photon and entangled-photon generations.
In this study, however, we focus on the optical response of LSPs driven by classical light in order to compare to conventional experimental results. In this case, the quantum Langevin equations can be reduced to the optical Bloch equations, by calculating the expectation values of operators, described as As in the second quantization of the LSP, we consider a cw incident light given by E(t) = E in,i exp(−iωt). In the case of the cw incident light, it is possible to describe the response of the LSP in terms of the steady state dipole expectation value ŝ σ 〈 〉 − oscillating at the same frequency as E(t). σ 〈 〉 ẑ then oscillates slowly as compared to the frequency ω, and hence z s σ 〈 〉 can be considered as a constant. The steady state dipole can then be obtained from Eq. (16) as is proportional to the steady state value of the inversion ˆz s σ 〈 〉. This steady state value can be obtained from Eq. (17), using the relation between dipole and inversion given by Eq. (18), and is given by where Ω = 2d i E in,i ћ −1 is the Rabi frequency. Using Eq. (19), σ 〈 〉 − ŝ can be rewritten as Optical response function of a single LSP. In order to compare the optical response of the plasmonic Bloch equations to the experimental results, we here introduce an optical response function χ of a single quantized LSP. The imaginary part of χ corresponds to the spectrum of a single LSP. From Eq. (21), χ is simply defined asˆε Though ω R and η are formally derived from Eqs (4) and (7), respectively, their values cannot be determined without an explicit expression of ε(ω). In what follows, for simplicity, we consider a noble nanometal and use the Drude model, where ω P is the plasma frequency and ε ∞ is the dielectric constant of the nanometal. From Eq. (4), ω R can be derived as The remaining unknown parameter is γ, however it cannot be derived directly from our method. Since γ corresponds to radiative relaxation of a single LSP, we evaluate it from the spontaneous emission rate, given by Using Eqs (11) and (25), |d i | 2 can be rewritten as Thus, the plasmon resonance frequency ω R , plasmon radiative decay rate γ, and optical response function χ for a single LSP can be expressed in terms of only the plasmon and geometrical parameters of nanometals, namely, ω P , a i , L i , ε m and ε ∞ .
Optical spectra of a single LSP: Dependence on metal particle size. Here, we numerically analyse the particle-size dependence of optical spectra obtained from a single LSP using Im [χ]. The details of calculation parameters are given in the section of Methods. Figure 3 shows Im[χ] for a x = 12.5, 25.0, 37.5 and 50.0 nm, normalized by its peak maximum at ω R . Both the central frequency ω R and FWHM 2γ decrease as a x increases. Details of the change in ω R and γ are shown in Figs 4 and 5, respectively, as a function of a x . One can find that ω R decreases monotonically as a x increases, whereas γ decreases gradually, in particular in the range of small aspect ratio a x /a < 2: From a x = 12.5 to 30 nm ω R decreases from 3.2 to 2.5 eV, whereas γ is nearly unchanged from 25.5 to 25 meV.
A similar tendency for γ to decrease as the aspect ratio increases can be found experimentally for gold nanorods 30 . For comparison with the experiment, the quality factors Q = ω R /2γ are plotted in Fig. 6. The Q values  observed in the experiments are, however, smaller, being around a third of that obtained from our method. In addition, the size dependence exhibits the opposite trend: Q obtained from our method gradually decreases with increase in a x , whereas Q in the experiment increases with increase in a x and saturates for large a x . This is due to the inhomogeneity of the nanometals in the actual experiment. Our result for size dependency of ω R indicates that, for small a x , the small difference in a x , Δa x , causes large inhomogeneous linewidth: For example, for a x < 20 nm, a size inhomogeneity of Δa = 0.5 nm leads to an energy fluctuation of Δω R ≈ 25 meV, comparable to γ. Therefore, for small a x , size inhomogeneities in the nanometals directly lead to overestimation of the FWHM of the spectra. On the other hand, for large a x > 50 nm, errors of Δa = 0.5 nm lead to Δω R < 10 meV, which is smaller than γ, and have little influence on the FWHM. Consequently, in the actual experiment, Q values for small a x can become smaller than theoretically predicted values owing to size inhomogeneity of nanometals, and exhibit a tendency to increase as a x increases. The size dependence of Q observed in the experiment can thus be explained qualitatively using ω R and γ for a single LSP, which can be obtained from our method.

Summary and Discussion
In summary, we have proposed a simple model of a saturable LSP using an effective dipole approximation. Taking an ellipsoidal nanometal as an example, we also have derived an optical response function of a single LSP from plasmonic Bloch equations, based on the optical Bloch equations, and have analysed in detail the size dependence of ω R and γ of a single LSP. We have shown that ω R and γ decrease as the aspect ratio of an ellipsoid increases, which is similar to the size dependence observed in the early experiments. We also have shown that, for small nanometals, ω R is very sensitive to the inhomogeneity of particle size, leading to broadening of the FWHM of spectra obtained from a single LSP.
Although in this study we adopt a simple two-level model to describe the saturation of LSPs, a more rigorous approach, e.g., using lattice model for finite systems, could also be developed. However, we would like to emphasize that the proposed model can be simply applied to various shapes of LSP by varying the semiaxes of a x , a y and a z , for example  a a a x y z = for a disk-like nanometal or =  a a a x y z for a rod-like nanometal, without explicitly considering the system size N of lattice model. In fact, the optical response of a single LSP is sensitive to the shape of nanometal, and hence our ellipsoidal model is quite useful. In addition, if we use the Drude-Lorentz model for the complex dielectric function ε(ω), the range of application of our method can be easily extended. Thus, our proposed method has advantages in simplicity and usability.
Finally, we should refer to the saturation of ω R for large a, strictly speaking, less dependence of ω R on a. This is because that ω R approaches to the bulk limit, in which boson model is approximately valid. In fact, for large semiaxis a x , ω R becomes almost constant. However, in this study, we presume a small nanometal to analyse the saturation effect of LSP and ignore higher-order multipole radiations. Therefore, if we would like to analyse the connection to the bulk limit rigorously, we have to consider the Mie scattering theory, which is beyond our model. Since the lowest-order approximation of scattering problem, used in this study, is adequate for nanoparticles of dimensions below 100 nm 28 , our proposed method is applicable for small nanoparticles. We hope that our results will facilitate the spread of the research field of quantum plasmonics.

Methods
Plasmonic parameters. In the calculation of optical spectra of a single LSP, we use the plasmon parameters of ω P = 11.586 eV and ε ∞ = 8.926 31 . For the background and semiaxes of the ellipsoidal nanometal we refer to ref. 30 : ε m = 2.25, a y = a z = a = 12.5 nm (fixed) and a x varies from 12.5 to 60 nm.  Randomly oriented dipoles. |d i | and γ in the main text are obtained for a single LSP oriented in the same direction as the polarization of incident light. In realistic systems, however, the LSPs might be randomly oriented. In order to more accurately evaluate d i and γ for comparison with the ensemble average of experimental results, we have to replace |d i | 2 → ∑ i |d i | 2 3 −1 in Eq. (26). In the actual calculation, randomly oriented dipoles are adopted.