Spasers with retardation and gain saturation: electrodynamic description of fields and optical cross-sections

Realistic spasers are numerically modeled within classical electrodynamics scattering framework using intensity-dependent Lorentzian dielectric function. Quantum mechanical effects are accounted for via saturation broadening. Spasers based on silver nano-shells and nanorods with strong field inhomogeneity and retardation are studied in detail. Fields and optical cross-sections are exhaustively analyzed upon variation of three control parameters: the amplitude of the gain Lorentzian, the detuning of the driving frequency from the spaser generation frequency, and the strength of the external E-field. An externally driven spaser demonstrates bistability for E-fields and optical cross-sections, while a freely generating spaser corresponds to the limiting case of vanishing external field. Gain saturation removes singularities and unphysical post-threshold behavior frequently reported with linear simulations. A small shift of the spaser generation frequency with increasing available gain level is observed. ©2015 Optical Society of America OCIS codes: (250.5403) Plasmonics; (240.6680) Surface plasmons; (350.4238) Nanophotonics and photonic crystals; (140.3430) Laser theory; (290.5825) Scattering theory; (190.1450) Bistability. References and links 1. D. J. Bergman and M. I. Stockman, “Surface plasmon amplification by stimulated emission of radiation: quantum generation of coherent surface plasmons in nanosystems,” Phys. Rev. Lett. 90(2), 027402 (2003). 2. X. G. Meng, U. Guler, A. V. Kildishev, K. Fujita, K. Tanaka, and V. M. Shalaev, “Unidirectional spaser in symmetry-broken plasmonic core-shell nanocavity,” Sci. Rep. 3, 1241 (2013). 3. X. L. Zhong and Z. Y. Li, “All-analytical semiclassical theory of spaser performance in a plasmonic nanocavity,” Phys. Rev. B 88(8), 085101 (2013). 4. N. Arnold, B. Ding, C. Hrelescu, and T. A. Klar, “Dye-doped spheres with plasmonic semi-shells: Lasing modes and scattering at realistic gain levels,” Beilstein J. Nanotechnol. 4, 974–987 (2013). 5. J. B. Khurgin and G. Sun, “Comparative analysis of spasers, vertical-cavity surface-emitting lasers and surfaceplasmon-emitting diodes,” Nat. Photonics 8(6), 468–473 (2014). 6. M. A. Noginov, G. Zhu, A. M. Belgrave, R. Bakker, V. M. Shalaev, E. E. Narimanov, S. Stout, E. Herz, T. Suteewong, and U. Wiesner, “Demonstration of a spaser-based nanolaser,” Nature 460(7259), 1110–1112 (2009). 7. X. Meng, A. V. Kildishev, K. Fujita, K. Tanaka, and V. M. Shalaev, “Wavelength-tunable spasing in the visible,” Nano Lett. 13(9), 4106–4112 (2013). 8. N. M. Lawandy, “Localized surface plasmon singularities in amplifying media,” Appl. Phys. Lett. 85(21), 5040– 5042 (2004). 9. I. E. Protsenko, A. V. Uskov, O. A. Zaimidoroga, V. N. Samoilov, and E. P. O’Reilly, “Dipole nanolaser,” Phys. Rev. A 71(6), 063812 (2005). 10. A. Y. Smuk and N. M. Lawandy, “Spheroidal particle plasmons in amplifying media,” Appl. Phys. B 84(1-2), 125–129 (2006). 11. J. A. Gordon and R. W. Ziolkowski, “The design and simulated performance of a coated nano-particle laser,” Opt. Express 15(5), 2622–2653 (2007). 12. I. E. Protsenko, “Quantum theory of dipole nanolasers,” J. Russ. Laser Res. 33(6), 559–577 (2012). #248006 Received 17 Aug 2015; revised 23 Sep 2015; accepted 23 Sep 2015; published 15 Oct 2015 © 2015 OSA 1 Nov 2015 | Vol. 5, No. 11 | DOI:10.1364/OME.5.002546 | OPTICAL MATERIALS EXPRESS 2546 13. T. A. Klar, A. V. Kildishev, V. P. Drachev, and V. M. Shalaev, “Negative-index metamaterials: Going optical,” IEEE J. Sel. Top. Quantum Electron. 12(6), 1106–1115 (2006). 14. M. Wegener, J. L. García-Pomar, C. M. Soukoulis, N. Meinzer, M. Ruther, and S. Linden, “Toy model for plasmonic metamaterial resonances coupled to two-level system gain,” Opt. Express 16(24), 19785–19798 (2008). 15. S. Xiao, V. P. Drachev, A. V. Kildishev, X. Ni, U. K. Chettiar, H. K. Yuan, and V. M. Shalaev, “Loss-free and active optical negative-index metamaterials,” Nature 466(7307), 735–738 (2010). 16. M. I. Stockman, “Spaser action, loss compensation, and stability in plasmonic systems with gain,” Phys. Rev. Lett. 106(15), 156802 (2011). 17. N. Arnold, C. Hrelescu, and T. A. Klar, Institute of Applied Physics, Johannes Kepler University, Altenberger Straße 69, 4040, Linz, Austria, are preparing a manuscript to be called “Minimal Spaser Threshold within Electrodynamic Framework: Shape, Size and Modes.” 18. M. I. Stockman, “Nanoplasmonics: past, present, and glimpse into future,” Opt. Express 19(22), 22029–22106 (2011). 19. A. E. Siegman, Lasers (University Science Books, Mill Valley, Calif., 1986). 20. D. G. Baranov, E. S. Andrianov, A. P. Vinogradov, and A. A. Lisyansky, “Exactly solvable toy model for surface plasmon amplification by stimulated emission of radiation,” Opt. Express 21(9), 10779–10791 (2013). 21. T. Klar, M. Perner, S. Grosse, G. von Plessen, W. Spirkl, and J. Feldmann, “Surface-plasmon resonances in single metallic nanoparticles,” Phys. Rev. Lett. 80(19), 4249–4252 (1998). 22. M. I. Stockman, “The spaser as a nanoscale quantum generator and ultrafast amplifier,” J. Opt. 12(2), 024004 (2010). 23. V. M. Parfenyev and S. S. Vergeles, “Quantum theory of a spaser-based nanolaser,” Opt. Express 22(11), 13671–13679 (2014). 24. A. P. Vinogradov, E. S. Andrianov, A. A. Pukhov, A. V. Dorofeenko, and A. A. Lisyansky, “Quantum plasmonics of metamaterials: loss compensation using spasers,” Phys-Usp 55(10), 1046–1053 (2012). 25. E. S. Andrianov, D. G. Baranov, A. A. Pukhov, A. V. Dorofeenko, A. P. Vinogradov, and A. A. Lisyansky, “Loss compensation by spasers in plasmonic systems,” Opt. Express 21(11), 13467–13478 (2013). 26. S. Wuestner, A. Pusch, K. L. Tsakmakidis, J. M. Hamm, and O. Hess, “Overcoming losses with gain in a negative refractive index metamaterial,” Phys. Rev. Lett. 105(12), 127401 (2010). 27. A. Pusch, S. Wuestner, J. M. Hamm, K. L. Tsakmakidis, and O. Hess, “Coherent amplification and noise in gainenhanced nanoplasmonic metamaterials: a Maxwell-Bloch Langevin approach,” ACS Nano 6(3), 2420–2431 (2012). 28. N. Arnold, L. J. Prokopeva, and A. V. Kildishev, “Modeling the local response of gain media in time-domain,” Annual Review of Progress in Applied Computational Electromagnetics 29, 771–776 (2013). 29. J. Pan, Z. Chen, J. Chen, P. Zhan, C. J. Tang, and Z. L. Wang, “Low-threshold plasmonic lasing based on high-Q dipole void mode in a metallic nanoshell,” Opt. Lett. 37(7), 1181–1183 (2012). 30. R. W. Boyd, Nonlinear optics (Academic Press, 2003). 31. L. Novotny and B. Hecht, Principles of Nano-optics (Cambridge University Press, 2012). 32. C. F. Bohren and D. R. Huffman, Absorption and Scattering of Light by Small Particles (Wiley, 1983). 33. J. A. Wiser, “Harmonic resonance dynamics of the periodically forced hopf oscillator,” in Mathematics (Ohio State Unversity, 2013), p. 178. 34. E. S. Andrianov, A. A. Pukhov, A. V. Dorofeenko, A. P. Vinogradov, and A. A. Lisyansky, “Forced synchronization of spaser by an external optical wave,” Opt. Express 19(25), 24849–24857 (2011). 35. S. Y. Liu, J. Li, F. Zhou, L. Gan, and Z. Y. Li, “Efficient surface plasmon amplification from gain-assisted gold nanorods,” Opt. Lett. 36(7), 1296–1298 (2011). 36. H. P. Zhang, J. Zhou, W. B. Zou, and M. He, “Surface plasmon amplification characteristics of an active threelayer nanoshell-based spaser,” J. Appl. Phys. 112, 074309 (2012). 37. P. Ding, J. N. He, J. Q. Wang, C. Z. Fan, G. W. Cai, and E. J. Liang, “Low-threshold surface plasmon amplification from a gain-assisted core-shell nanoparticle with broken symmetry,” J. Opt. 15(10), 105001 (2013). 38. W. R. Zhu, M. Premaratne, S. D. Gunapala, G. P. Agrawal, and M. I. Stockman, “Quasi-static analysis of controllable optical cross-sections of a layered nanoparticle with a sandwiched gain layer,” J. Opt. 16(7), 075003 (2014). 39. J. Song, Y. L. Tian, S. Ye, L. C. Chen, X. Peng, and J. L. Qu, “Characteristic analysis of low-threshold plasmonic lasers using Ag nanoparticles with various shapes using photochemical synthesis,” J. Lightwave Technol. 33(15), 3215–3223 (2015). 40. X. Fan, Z. Shen, and B. Luk’yanchuk, “Huge light scattering from active anisotropic spherical particles,” Opt. Express 18(24), 24868–24880 (2010). 41. A. A. Lisyansky, E. S. Andrianov, A. V. Dorofeenko, A. A. Pukhov, and A. P. Vinogradov, “Forced spaser oscillations,” Proc. SPIE 8457, 84570X (2012). 42. U. Kreibig and M. Vollmer, Optical Properties of Metal Clusters (Springer, 1995). 43. G. Raschke, “Molekulare Erkennung mit einzelnen Gold–Nanopartikeln,” in Fakultät für Physik (Ludwig– Maximilians–Universität, 2005), p. 146. 44. M. Quinten, Optical Properties of Nanoparticle Systems: Mie and Beyond (Wiley-VCH, Weinheim, Germany, 2011). #248006 Received 17 Aug 2015; revised 23 Sep 2015; accepted 23 Sep 2015; published 15 Oct 2015 © 2015 OSA 1 Nov 2015 | Vol. 5, No. 11 | DOI:10.1364/OME.5.002546 | OPTICAL MATERIALS EXPRESS 2547 45. P. B. Johnson and R. W. Christy, “Optical constants of noble metals,” Phys. Rev. B 6(12), 4370–4379 (1972). 46. M. Born and E. Wolf, Principles of Optics: Electromagnetic Theory of Propagation, Interference and Diffraction of Light (Cambridge University Press, 1999). 47. W. Liu, A. E. Miroshnichenko, R. F. Oulton, D. N. Neshev, O. Hess, and Y. S. Kivshar, “Scattering of core-shell nanowires with the interference of electric and magnetic resonances,” Opt. Lett. 38(14), 2621–2624 (2013). 48. B. S. Luk’yanchuk, N. V. Voshchinnikov, R. Paniagua-Dominguez, and A. I. Kuznetsov, “Optimum forward light scattering by spherical and spheroidal dielectric nanoparticles with high refractive index,” ACS Photonics 2(7), 993–999 (2015). 49. H. Kuwata, H. Tamaru, K. Esumi, and K. Miyano, “Resonant light scattering from metal nanoparticles: Practical analysis beyond Rayleigh approximation,” Appl. Phys. Lett. 83(22), 4625–4627 (2003). 50. V. M. Parfenyev and S. S. Vergeles, “Intensity-dependent frequency shift in surface plasmon amplification by stimulated emission of radiation,” Phys. 


Introduction
Surface plasmon amplification by stimulated emission (spaser) has been extensively discussed in the literature over the last decade theoretically [1][2][3][4][5] and experimentally [6,7].Similar concepts were discussed for nano-lasers [8][9][10][11][12] and loss compensation in metamaterials [13][14][15].Detailed analysis often relies on a quantum-mechanical description of the gain medium as a collection of individual chromophores [16].In practical realizations, it is the concentration of active agents, rather than their individual quantum mechanical properties, what limits the feasibility of a spaser.Hundreds to thousands of lasing molecules are required for generation in a realistic finite-sized spaser, especially taking into account that losses in metallic nanostructures increase with decreasing size due to boundary scattering and other parasitic effects.This justifies the description of a spaser based on macroscopic electrodynamics.The threshold conditions obtained in this way from the linear theory [17] are fully consistent with the more elaborate quantum mechanical results discussed in [18].However, linear theory becomes inadequate above threshold.For example, arbitrarily high field values are obtained in numerical studies near the threshold singularity, and fields and the scattering cross-section start to decrease, as available gain exceeds threshold value.From the physical point of view, the post-threshold operation of a spaser (or laser) is governed by a nonlinear gain saturation.A standard way to model the saturation in continuous wave (CW) lasers [19] is using an intensity-dependent dielectric function, where a Lorentzian gain includes saturation broadening.This theoretical approach was applied to the quasi-static (QS) core-shell (CS) dipolar spaser by Baranov et al. [20].In this system, the field in the gain core encapsulated within a metal shell is constant even in the saturated regime, thus allowing for an analytical treatment.
In the present work, we extend this approach to arbitrary finite-sized structures and multipolar modes with spatially inhomogeneous saturation, retardation and radiative losses.In doing so, we present our findings within the conventional framework based on optical crosssections, convenient for comparison with experiments.Saturation introduces non-linearity into Maxwell's equations written for a fixed frequency, which in the general case are solved numerically, using the finite-element frequency-domain (FEFD) nonlinear COMSOL Multiphysics solver.This quasi-static approach requires pulse durations longer than times required to establish a regime of stationary spasing, which are of the order of few picoseconds [21,22].For shorter pulses, time dependent analysis, based either on density matrix formalism [22,23], or similar equations of quantum-mechanical dynamics [24,25], or multilevel rate equations coupled with the quasi-classical description of local polarizations [26][27][28] should be used.Our framework is equivalent to a stationary cases of these approaches, and possesses the usual advantages of the FEFD methods, such as weaker stability constrains with direct solvers, conformal tetrahedral meshes, high-order basis functions, all leading to a possibility of quick and efficient analysis of much larger structures and parameter ranges than those which are accessible with time domain simulations.
The paper is organized as follows.Sec. 2 deals with the quasi-static core-shell dipolar spaser with gain saturation.The analytical results are discussed in Sec 2.1.In, Sec.2.1.1,we study the dependences on frequency and external field, in Sec.2.1.2,we study the dependences on frequency and available gain, and compare the behavior of optical crosssections in the linear and saturated cases.In Sec.2.2, we compare analytical and numerical results, thus validating our numerical framework.In Sec. 3, we apply the numerical framework to a spaser based on a quadrupolar mode of a large silver shell, sandwiched between a gain core and an outer gain shell.Such a spaser has strong spatial field inhomogeneity combined with retardation effects, which makes an analytical treatment impossible.In Sec. 4, we numerically investigate spasing of the dipolar mode of a small prolate silver spheroid immersed in gain material.In Sec.4.1, we study the influence of the size and geometry of a gain capsule enclosing the spheroid.In Sec.4.2, we discuss the changes of the spaser generation frequency with available gain.Finally, the conclusions summarize the results from the different sections.To facilitate the readability, the necessary derivations and auxiliary materials are put into 4 appendixes.Appendix A relates parameters used in the calculations to the microscopic properties of the gain molecules.Appendix B contains the reference formulas for the quasi-static core-shell spaser, while in Appendix C, we derive parameters of such a spaser near the generation threshold with and without an external field.Appendix D includes additional Figs., which further elucidate the interrelation between multiple parameters, fields, and optical cross-sections and briefly discusses the usage of spasers for loss compensation in metamaterials.

Analytical results
2.1.1Core field and scattering cross-sections as a function of frequency and external field Spasers, as well as nano-and micro-lasers are often studied within the linear electrodynamic scattering framework [2,4,29].Such an analysis correctly predicts the threshold gain level and the generation wavelength, but it becomes unphysical in the post-threshold operation regime [17,18], which, similarly to conventional lasers, is governed largely by gain saturation.
Adhering to macroscopic electrodynamics, we describe the gain saturation of the lasing transition by an intensity-dependent Lorentzian dielectric function G adapted from [30,31].For simplicity, we use the resonant approximation here, but the full version of Lorentzian, or more complex intensity-dependent line profiles can be straightforwardly used in the same framework.
Here h ε is the background "host" value; the Lorentzian is centered at L ω and has a native full where sat E is the saturation field.Qualitatively, the weakening is more important than the broadening, as it limits the intensity of a spaser; but the behavior of the spaser driven by an external harmonic wave is quite sensitive to frequency detuning from the spaser generation frequency [20], which makes the broadening crucial for the quantitative analysis.The relation of sat E and other parameters to the atomic characteristics and the pumping rate is discussed in Appendix A.
The dielectric function (1) makes Maxwell's equations nonlinear.If the gain saturation is spatially inhomogeneous, their analytical solution is problematic.If the field is homogeneous within the gain material, the non-linearity remains purely algebraic and the problem can be tackled analytically.This was done for the small core-shell spaser with a gain core and a metallic shell [20], utilizing the fact that the field in the core is constant in the quasi-static approximation.This holds also for (concentric and co-oriented) multi-shell ellipsoids, as can be deduced from the Ref. [32], p. 149 top, also pp.142, 148 and Eq.(5.29).
In this section, we recapitulate the analytical results from [20], recast them into the scattering framework convenient for further use, derive several functional dependences for the free and driven quasi-static spaser, and compare these predictions with numerical results for a small finite-sized core-shell spaser with retardation.
The core-shell geometry is shown in the Fig. 1, in the left panel, and the typical field distribution in the right panel.The constant field 1 e inside the core can be derived from electrostatic equations (Appendix B), and is related to the external field 3 e in the following way: Fig. 1.Geometry (left) and typical field distribution (right) for the (almost) quasi-static dipolar core-shell spaser.Gain core radius a 1 = 10 nm, Ag shell thickness h = 2.7 nm, background material is air.The right panel shows the distribution of the normalized field amplitude ⏐(E)⏐/E sat in the plane of incidence defined by the k-vector and external field vector e 3 , for a post-threshold gain, ε L = 1.5ε thr and an external field amplitude e 3 = 0.006 E sat .The shape of the field distribution is similar for other parameter values and (in the vicinity of the structure) is similar to that in a freely generating spaser.The field e 1 inside the core is practically constant.
The denominator D is the determinant of the coefficients matrix for the linear equations originating from the boundary conditions at two interfaces.With a saturating 1 ω ε of the complex-valued equation: Here, the r.h.s., and, if relevant, h ε , should be treated as functions of thr ω (see Eqs. ( 30)- (31) in Appendix C.1).
The following parameters were used in the calculations: Shell material is silver with , where p 4.9, 9.5eV, =0.05eV These values were taken from [20], to allow qualitative comparison with their results.For these values, spasing threshold, has the lowest values, corresponding to the minimum  , both of these dependences can be approximated by a single cubic equation with a complex constant C (for derivation see Eqs. ( 37)- (39) in Appendix C.2): C E e = follows from (5) and is indicated by the arrow in Fig. 2(a).The freely generating spaser without an external field, 3 0 e = , corresponds to the tip of the conical surface in Fig. 2(c ε above threshold, in agreement with findings in [20,22,23]: ε ε = used here, this region is rather narrow, within ± 0.2% of the generation frequency, similarly to the observations made in [20].This necessitates accurate accounting for frequency in numerical calculations, especially near the threshold (the bistability region increases with the available gain L ε ).Of course, in real systems, the experimental inhomogeneities and variance of the particle sizes and shapes will lead to a variation of individual spaser frequencies.This will broaden all fine frequency dependences in an ensemble of particles.However, the spectral predictions remain valid for each individual particle.The projection of the bistable solution surface onto the plane  A full-scale stability analysis for the solutions shown in Fig. 2 should rely on dynamic equations, which is beyond the scope of the present article.Nevertheless, we would like to briefly comment on it.A multi-valued solution surface shown in Fig. 2(c) contains stable and unstable regions.From the mathematical analysis of forced Hopf oscillators we expect that in the bistable region the upper branch is stable, while the middle solution is saddle-like, and the lower branch is unstable [33].This qualitatively agrees with the numerical stability analysis carried out in [25,34] on the basis of quantum-mechanical dynamics for a similar system, and in [20] using linear perturbation analysis based on the spectral properties of the polarizability. .This corresponds to the zero of the complex resonant denominator D in (2), or its more intricate versions with retardation studied in [17].Within the linear theory, arbitrarily high values can be obtained numerically near the singularity, and all the quantities start to decrease again as available gain exceeds threshold value.This unphysical artefact was often reported in numerical calculations performed within the electromagnetic scattering framework [2,4,[35][36][37][38][39]. A similar effect was observed for particles where a gain in the form of 0 ε ′′ < was introduced into the dielectric components of an anisotropic hyperbolic permittivity [40].In contrast, saturation  The behavior of other cross-sections is more complex.The scattering cross-section in Fig. 4(b) monotonously decreases with the external field everywhere; the absorption cross-section in Fig. 4(c) inverts from negative to positive values, starting from the center of the line, reaches a certain maximum value, and then gradually decreases with the external field, while its values at the wings become positive everywhere.The behavior of the extinction crosssection in Fig. 4(a) is influenced by the fact that in the linear problem the singularity is asymmetric (solid curve in Fig. 4(a)).With gain saturation, this singularity is smoothened as the external field 3 e increases, and the extinction remains positive everywhere to the red side of the resonance.At the blue side, ext 0 σ ≤ in the region indicated by the dashed ellipse; this region depends on the external field 3 e .For small fields 3 e , weakly saturated gain can overcompensate the losses, while for larger 3 e , the gain saturates, and ext 0 σ > everywhere (dash-dotted curve in Fig. 4 , but because all optical cross-sections depend continuously on parameters, similar behavior can be observed within a certain gain interval both below and above the threshold.More details are given in Appendix D, Figs.12-14.

Comparison between the analytical and numerical results
To validate our computational framework, we compare the analytical predictions for the quasi-static dipolar core shell spaser from the previous section with the numerical solution of nonlinear Maxwell's equations.The latter is obtained using the Electromagnetic Waves, Frequency Domain Solver (COMSOL Multiphysics, Wave Optics Module).
Let us first discuss the threshold values.The quasi-static threshold condition (4) does not include radiative losses and inhomogeneous depolarization within the core.These retardation effects are included into the dipolar scattering coefficient in the rigorous multi-shell Mie theory [42][43][44].Equating the corresponding denominator to zero, one obtains the threshold values with retardation: thr thr thr Mie Mie Mie 2.4840 eV ( 499.13 nm), 0.1683 The generation frequency is somewhat red-shifted as compared to the quasi-static values (4), and the threshold gain level is noticeably larger.Thus, even for the rather small structure studied here (shell outer radius 2 12.7 nm a = ), the retardation is not negligible.In our numerical simulations, we find: thr thr thr Num Num Num 2.4839 eV ( 499.15 nm), 0.1684 The numerical threshold is basically identical with the Mie values, and was used for normalization of numerical results.With retardation, the field 1 e and consequently, the saturation s vary inside the core, and analytical formulas describing them are not available.However, we expect that the behavior of the retarded system with respect to available gain and frequency normalized to the threshold values L t h r Mie / ε ε and thr Mie / ω ω will be similar to the non-retarded case with respect to L t h r t h r QS QS / , / ε ε ω ω .This comparison is shown in Fig. 5, where solid curves are obtained from the analytical quasi-static solution (2), and full symbols are obtained from FEFD modeling.
Both results demonstrate good agreement, when the normalization to their respective threshold values is used.(where only the upper curve in Fig. 5(a) exists), the upper branch of solution can be numerically explored.
We find generally a close match between the analytical and numerical results.Slight discrepancy for larger fields 3 e , especially near the disappearance of the lower solution branch are related to the fact, that the analytic results are obtained with quasi-static formulas, while our numerical results are fully retarded.

Large shell quadrupolar spaser
For the system studied above, the field inside the gain material is almost constant and the existing analytical solution provides guidelines for numerical studies.Let us now numerically model larger, inherently inhomogeneous CS structure with retardation, for which up to now no analytical solution has been reported.As a first example we will discuss silver shells immersed in gain material.Optical properties of silver are taken from [45].We choose a structure with a quadrupolar resonance around 756 nm, where an essential (though not global) optimum for spasing exists.To make the system finite, we use a finite outer gain layer (see inset in Fig. 6(b)).This arrangement utilizes gain on both sides of the metal, which decreases the threshold and makes it close to the case with the infinite gain background.The quadrupolar mode has lower radiative losses than the dipolar one, and hence the threshold value remains close to its minimal quasi-static limit even for the relatively large structure used below [17].
We choose the sizes of the core and both shells, as well as the parameters of the gain Lorentzian (1) such, that the metal is not excessively thin, and that the quadrupolar resonance is close to the optimal wavelength, namely: core radius 1 50 nm a = , Ag shell thickness For the sake of generality, we assume a spatially homogeneous available gain, i.e., L ( ) const ε = r in the absence of lasing saturation.This holds e.g., in the case of full pumping saturation, when all available chromophores are in the upper lasing state.The threshold, obtained theoretically from the zero of a denominator in the quadrupolar Mie scattering coefficient for the multi-shell structure [42][43][44] However, due to discretization, the numerically found threshold has a slightly lower frequency thr thr thr Num Num Num 1.6368 eV ( 757.49 nm), 0.03084 The numerical threshold was used for normalization.The quadrupolar mode has an inherently inhomogeneous field distribution even inside the core, and hence the analytical approach similar to that used in [20] cannot be applied.The fields are shown in the inset to Fig. 6 The overall behavior of the system is similar to that of the small CS structure in Fig. 5.The field is now strongly inhomogeneous (see inset in Fig. 6(a)), and the plotted maximum field in the gain material, right outside the metallic shell is higher for the same over-pumping L t h r 1.5 The far field behavior of the scattered field closely follows a standard quadrupolar pattern [32,46].This is related to two factors.First, especially near the threshold, the gain compensates better the damping of the resonant quadrupolar mode, as compared to non-resonant modes, making this mode even more dominant [4].Second, for the numbers used here with low-n materials and relatively small structures, an interference between the electric and magnetic modes, which can lead to pronounced changes in far field directivity [47,48], is not observed.This holds also for the dipolar structures studied in Secs.2.2 and 4. , which is about trice higher than for the small dipolar CS spaser in Fig. 5(a).The upper branch does not extend up to the freely generating regime 4 0 e = , but disappears around 4 s a t 0.00138 e E = .This is due to slight detuning from the (numerical) threshold frequency thr,Num ω from (10), which translates into the relative detuning of , the lower branch is the only existing solution.Seeding of a spaser by a weak, slightly detuned external field in this parameter range results in amplification.The amplification coefficient follows from the Eq. ( 5) (see also Eq. ( 46) Appendix C.3), and is large near threshold: If the (detuned) seeding is done with a larger field, 4 s a t 0.00138 e E > , we expect the upper, spaser-like, solution to be realized in practice., as labeled in the plot.As the frequency is gradually changed, one sees abrupt changes in the gain field at the points where the upper "spaser-like" solution appears or disappears.
Figure 7 shows the dependence of all cross-sections on the external field at three values of available gain: L t h r 0.5, 1, and 1.5 The behavior is complex, especially above the threshold, where the bistability exists.The details of this behavior can be understood by comparing with the quasi-static case shown in Fig. 2(d) for the scattering and Figs.12(c) and 12(d) for the absorption and extinction at small blue detuning from the generation frequency.As the external field decreases (Fig. 7(c)), the solution first follows the "generation" branch with increasing scattering and exceedingly negative absorption, while the extinction first increases and then decreases into the negative range.However, at any finite detuning from thr ω , for small enough external field, the "generation" branch ceases to exist, and the low field "amplification" branch is the only available.Alternatively, with an increasing external field, one starts from the lower branch, and the transition to the higher "generation" state occurs for higher external fields, thus demonstrating a hysteresis, though we do not rigorously analyze here the stability of the branches.In the region 4 s a t 0.00138 e E < (as in Fig. 6(a)), the only existing solution has both abs 0 σ < and ext 0 σ < (see are marked dashed in Figs.7(c) and 7(d)).For these external fields 4 e , the considered CS structure is expected to operate as a stable amplifier.

Influence of the thickness of the gain capsule
Minimal thresholds can be realized by a simpler geometry of a metallic spheroid immersed into gain media.The dipolar modes of such spheroids can be tuned to a wavelength optimal for spasing by changing their aspect ratio.In spite of radiative losses, the thresholds remain close to the quasi-static minimum for long semi-axes up to 20-30 nm [17].First we consider prolate silver spheroids with the optical properties from [45].The length of the longer semiaxis is γ .The geometry of this structure is shown in the left inset in Fig. 8(a).The threshold estimated theoretically from the zero of the dipolar spheroidal denominator with retardation corrections in an infinite gain medium [17,49] However, as the outer material is air with 3 1 ε = , the numerically found threshold has slightly higher frequency and noticeably lower threshold, due to lower radiative losses into air background, as compared to the infinite background with h 2.6 ε = , for which the estimation (12) This numerical threshold was used for normalization.
The upper curve in Fig. 8 The narrow field confinement near the spheroidal structure suggests the following practical question: how little gain material is sufficient, for the threshold to remain close to the small value obtained for the 2 80 nm a = case?We did not investigate this exhaustively, but the estimations can be obtained from the field distribution shown in the upper panels of Fig. 8(a), or from the quasi-static formulas for the coated ellipsoid [32], Eq. (5.35).To account for finite sizes, and stay close to the possible experimental realizations, we numerically found the threshold for a similar Ag spheroid inside a spheroidal gain capsule with the same optical properties as before, but with the thickness 2 10 nm  h = along all semiaxes.For a thin gain shell, the resonance shifts to the blue.In order to keep the resonance in the same spectral range (where the Ag optical properties are almost optimal), a thinner spheroid was used, with the same long semi-axis 20 The lower curve in Fig. 8(a) shows the dependence of the maximum gain field on the external field, for the spaser with thin gain shell.The geometry is shown in the right inset in the Fig. 8(a), while the field distribution can be seen in the right part of the upper panel.The field distributions, as well as the initial increase in the gain field with available gain are comparable.Note, however, that the threshold ( 14) for the thin gain shell was about 20% higher than for the thicker gain (13), so that the absolute levels of gain parameter L ε are respectively higher as well.Nevertheless, the threshold gain in ( 14) remains well within reach of commercially available dye-doped polymers.An interesting observation is the difference in the behavior of the maximum field in two cases.While for the thin gain shell (red curve) max sat in agreement with the Eq. ( 6), for the thicker gain (black curve) this dependence is close to linear.This is a somewhat unexpected result, probably indicating, that as the saturation depletes the population inversion in the immediate vicinity of the spheroid, in the case of a thick gain shell, more remote, non-depleted regions start contributing more to the spasing, increasing the maximal field.

Changes in the generation frequency with the available gain
The quasi-static spaser generation frequency does not depend on the available gain level L ε .
This has been obtained in [22] and [20], and requires certain conditions, see Eqs. ( 32)-( 34) in Appendix C.1.This may sound counterintuitive, as the frequency of a linear oscillator does appreciably change with damping.The physical explanation is that a spaser as a nonlinear generator establishes an amplitude, which exactly nullifies the gain/losses over one oscillation cycle.However, the aforementioned derivations in Refs.[20,22].or Appendix C.1 rely on the absence of retardation and homogeneity of the field in the gain medium.The dependence of the generation frequency on the spaser intensity was discussed in [50].It is based on "spatial hole burning", i.e., essentially on changes in the spatial profile of G ε with saturation, which influence the electrodynamic part of the problem.This effect is also expected for a highly nonlinear and spatially inhomogeneous retarded spheroidal spaser.The spaser generation frequency can be determined numerically in several ways.First, as a singularity of the linear problem, where all the fields diverge (see Fig. 3(a), 3(c)), and second, as a frequency where the bistable cone touches the plane of zero external field ext 0 e = for a given gain level (see Fig. 2(c)).In practice, this can be found as a spectral maximum (or middle) of the small-field response (the innermost curve with 4 s a t / 0 .0 0 1 5 e E << in the Fig. 6(b)).An even better numerical recipe is provided by the Figs.

Conclusions
The spaser is an essentially non-linear system, but its CW operation, similarly to that of a laser, can often be understood within macroscopic electrodynamics, if gain saturation is accounted for.The approach developed in this paper is based on non-linear Maxwell equations with an intensity dependent dielectric function 2 G ( ( ) ) ε E r for the gain material.
The numerical framework uses a finite element, frequency domain nonlinear solver (Wave Optics Module, COMSOL Multiphysics).Here, the spectral gain profile follows Lorentzian dielectric function with saturation broadening.This approach allows realistic modeling of arbitrarily shaped, finite-sized spasers with retardation, radiative losses and spatially inhomogeneous field and saturation.Both the freely generating, and optically driven regimes are investigated.The optically driven spaser is governed by three control parameters: the frequency of the incident wave, the external field, and the available gain.The latter is quantified by the amplitude of the Lorentzian gain dielectric function, which includes both the concentration of the gain molecules and the pumping level.Different spasers demonstrate similar qualitative dependence of the spaser field on these parameters, if they are normalized to the threshold values and the saturation field for the gain material.These functional dependences are found from the analysis of the quasi-static core-shell spaser with a gain core and metallic shell [20].The results are presented in the conventional electrodynamic scattering framework, convenient for the comparison with the experiments.The freely generating spaser is modelled as a driven spaser with the vanishingly small external field.The numerical results for the quasi-static spaser are in good agreement with the theory.However, as the spaser response is very sensitive to the frequency detuning from the exact resonance, one has to include retardation effects into realistic modeling.
In distinction to linear analysis, non-linear optical gain saturation removes the threshold singularities in the fields and optical cross-sections, and makes the internal fields and scattering cross-section at resonant frequency monotonously increasing with the available gain.Simultaneously, the spectral properties of an externally driven spaser above threshold experience discontinuous changes if the driving frequency is scanned.The optical crosssections of the driven spaser depend on the value of the external field.This is especially pronounced for small excitation fields above the spasing threshold.
We modelled two realistic spaser geometries: a) A thin silver shell of finite dimensions (inner radius 50 nm, thickness 4.3 nm) sandwiched between the core and second shell of identical gain material.Here the strongly inhomogeneous quadrupolar mode was chosen for spasing.b) A dipolar spaser, based on prolate silver spheroids (large semi-axis of 20 nm, aspect ratio ~4) immersed into a spherical or spheroidal gain shell.all cases, shape optimization allowed to achieve low thresholds, close to the analytically predicted minimal values for the quasi-static systems [17,18].In the case of spheroids, a gain layer of 10 nm increases the threshold by only about 20% as compared to the infinite gain surrounding.
The results suggest the use of multipolar modes, which have lower radiative losses (in relative terms), and allow larger structure sizes.Larger structures are easier to manufacture, and they have lower surface scattering of the electrons and better surface roughness (which both adversely affect the spasing thresholds).In addition, the fraction of gain material affected by the non-radiative quenching of chromophores in the immediate vicinity of metal surfaces is decreased in larger multipolar spasers, compared to the very small spasers, based on the dipolar modes of the MNP.Finally, larger mode volume may allow larger number of plasmons in a spasing mode.
The quasi-static spaser generation frequency is typically independent of the available gain level.We find that this is not the case for structures with spatially inhomogeneous gain saturation and retardation.For the silver spheroid, a relative red shift of about 4 10 − is obtained numerically for gain levels few times higher than the threshold value.

Appendix A: Parameters of the dielectric function, atomic characteristics and pumping rate
Let us relate the parameters of the Lorentzian (1) to the pumping rate and the microscopic characteristics of the gain material.The gain material shall be a four level system; the levels are numbered from 0 to 3, from low to high energy, with i N being the respective concentrations.Let us further assume fast relaxation of the upper pumping level 3→2, and lower lasing level 1→0.These are realistic approximations for organic dyes and pulse durations longer than several ps.Then the populations Here N γ is the total relaxation rate of the population inversion (which may include both radiative and non-radiative contributions), in this case upper lasing level 2. For dyes, it typically combines the rates of radiative transition 2 The stationary solution of Eq. ( 15) is Here we introduced the unsaturated population inversion due to pumping, in the absence of lasing: This initial inversion increases linearly with p I for small pumping rates p  17) accounts for the saturation of the lasing transition.The lasing saturation intensity sat I is defined as lasing intensity, where the midband gain drops to half of its unsaturated value.There exists a difference between the low and high pumping rates.
The first equality holds for the idealized 4-level system discussed here, while for the 2-level system, often used in the literature, the relaxation rate N γ here and in all subsequent expressions should be divided by 2. For the estimations we used the following typical parameters: .For high pumping levels (which still obey the restriction ( 16)), the lasing saturation intensity is given by the pumping intensity.The turnover between these two regimes occurs at a pumping saturation intensity, similar to the one estimated in (19) for the emission line: Let us now express the saturation field via the (absolute value of) matrix element for the dipole moment 12 μ .We recall that the emission cross-section for the Lorentzian line with the FWHM L γ is [19] Eq. ( 7.5.62),or [51], Eq. ( 2 Here 12 σ refers to the macroscopic, orientation-averaged cross-section, while 12,max σ assumes an ideal alignment between the dipole moment and the incident field.n and ε are the refractive index and dielectric function the of the host matrix, while k is the vacuum wavenumber.For reference purposes, the last equality is written both in CGS and SI units. Substituting the expression ( 20) into (19), we obtain for the saturation field This expression for sat E holds in in both CGS and SI units.For the estimations we assumed .This coincides with the notations in [20], and is 6 smaller than the estimation (21).
The amplitude of the Lorentzian gain profile L ε in the expression ( 1) is proportional to the atomic polarizability α and to the equilibrium population inversion 2 (0) Comparison with the expressions (A.52) from [31] or (6.3.17) from [52] results in the following CGS expressions (they should be divided by 0 4πε in SI units): .From the Eq. ( 18) we see, that half of this value is reached at pumping saturation intensity p,sat sat , estimated after the Eq.(19).Such values of L ε lie above the spasing threshold near optimal conditions (see main text), while such pumping intensities are easily achievable even in the CW regime (e.g., 1 W laser focused into ).The near-stationary inversion is achieved on a time scale estimated from Eq. ( 15) as: Under normal conditions, this time scale is in the ns range.Purcell enhancement of spontaneous emission into the resonant (and non-resonant) plasmonic modes of small structures may drastically increase γ ∝ E r [5].This will shorten inv τ , and significantly modify the estimations ( 19), (21), resulting in the saturation parameters in the range .This will necessitate much higher pumping rates, which may lead to overheating issues [5,53].In addition, even in an ideal metal, surface scattering of the electrons [54], as well as surface roughness, may significantly increase losses and thresholds, and contribute to heating.Although all these issues are beyond the scope of the current paper, we note that heating problems may be circumvented by a pulsed excitation with sub-nanosecond pulses.As noted above, frequency domain calculations will still give the right answers, as long as the pulses are longer than typical vibrational relaxation times of the molecules, which are usually still shorter than the Purcell-enhanced relaxation rates.
In a more accurate Lorentz-Lorenz-type approach, the field acting on the dye molecule is the local field in the virtual spherical cavity.For realistic numbers, the changes in dielectric function due to active molecules (e.g., dyes) are small, which adds a pre-factor the gain Lorentzian and to the saturation term s in the denominator of expression (1).This is equivalent to the following rescaling of L ε and saturation intensity or field indicated by the additional subscript "LL": Finally, to give a link to a quantum mechanical description of a spaser [1,5,16,18,22,23,55], we note, that in the low post-threshold regime, say for ).The corresponding number of plasmons are of the order of (in SI units): Here V is the mode volume, which can be visually estimated from the right panel in Fig. 1, inset in Fig. 6(a), and upper panels in Fig. 8. Using Purcell-enhanced saturation field estimated after the Eq. ( 23) and the modal volume , Eq. ( 25) gives plasmon numbers pl ~1  n near the threshold.This is consistent with the predictions of quantummechanical analysis [5,22,23].
Appendix B: Quasi-static core-shell structure in an external field The geometry of the problem is shown in the right panel in Fig. 1; the subscript 1 refers to the core, 2 to the shell, and 3 to the ambient.The external field 3 E , all dipoles and constant field contributions are oriented in x-direction.Let us search for the quasi-static solution, which is a combination of constant and dipolar potentials/fields in different regions [56].
In these notations 1 e represents the field in the core (where the singular dipolar term must vanish), while These 4 equations can be solved for 4 parameters 1 2 2 3 , , , e e d d , which determine all fields as a function of the external field is 3 e .This leads to the expression (2) for core field in the main text.Other quantities have the same denominator, and can be found similarly.The dipole moment of the whole structure 3 d is: All cross-sections follow from here, or as the first terms of the general multipolar multi-shell Mie expansion, e.g., Eq. (2.8a) with (A.1-A.2) from [43], or Eq.(7.9) from [44]:

C.1. Spaser threshold and generation without an external field
Let us first analyze the freely generating spaser without an external field, 3 0 e = .In this case, the linear system (27), which leads to the expression (2) becomes homogeneous.A nontrivial solution of this system exists only when the corresponding determinant, (denominator D in (2)) equals zero.This becomes a condition for the normalized core intensity .The phase of the core field is of no importance for the generation analysis.The consistency condition 0 D = can be resolved with respect to 1 ε : ) This complex equation defines 2 real quantities: the normalized core intensity s and the generation frequency gen ω , considered as a function of all parameters and available gain L ε .
Let us study the threshold condition and the post-threshold operation separately.Slightly above threshold, the spaser field can be arbitrarily small.This means that the equality (30) should also hold for vanishing core intensity 0 s → .This defines the threshold frequency thr ω and available gain level L t h r Above the threshold, the Eq. ( 30) should hold for 0 s ≠ , so that at the generation frequency gen ω : In several important cases, we can easily find a physical solution of this equation.such, that the complex number on the l.h.s.does not change as well, the Eq. ( 32) is fulfilled automatically.This requires: which defines the intensity s.For Lorentzian gain (1) this results in: This value of s together with gen thr ω ω = provide a solution of the complex equation (two real equations) (32).Thus, above the threshold, the generation frequency does not change, while the intensity s is linear with the available gain L ε .The slope of this linear dependence is defined by (34).Note, that the threshold parameters thr ω and thr ε in this dependence should be found from Eq. ( 31), which depends on the dispersion of all materials and the geometry.This derivation remains valid for arbitrary dependence G L ( , , ) s ε ω ε , as long as the gain material can be described by a single spatially constant field (or saturation parameter s), and if the solution of the Eq. ( 33) results in real-valued s.The first requirement holds e.g., for the ellipsoid gain cores (see section 2.1.1,above Eq.( 2)), or gain material in a form of a point dipole.The second requirement is fulfilled for example for the full (non-resonant) Lorentzian profile used in [20].It does not hold, however, for a double Lorentzian profile, which consists of a "positive" absorption Lorentzian and a "negative, inverted" emission Lorentzian, as used in [4] to describe spectral properties of lasing dye chromophores.This may result in gen L thr thr ( ) for such systems, even in the quasi-static geometries with gain in the ellipsoid core or point dipole.

C.2. Near-threshold behavior of a spaser driven by an external field
In the presence of an external harmonic field 3 0 e ≠ , all quantities remain finite even when the spaser is driven at its native generation frequency e .There exists a phase difference between 3 e and 1 e , or other quantities.As 1 e enters all expressions non-linearly, it is usually more convenient to treat 1 e as real and calculate the phases "backwards".
To obtain general dependences, we rewrite (2) with respect to its denominator D: Exactly at the threshold, the denominator thr 0 D = .We expect all physical quantities to be continuous close to the threshold.Moreover, we expect that in this region 1 3 e e >> , as the spaser can generate even with 3 0 e = .These assumptions are confirmed by the numerical studies and self-consistent results of this Appendix.This makes the r.h.s. in (35) small, and we can Taylor-expand D on the l.h.s.Here, all dielectric functions This is a cubic equation for 1 Note, that C is a complex number, and a phase difference exists between 1 e and 3 e .Exactly at the threshold L 0 δ = , and the core field from (37) becomes: The denominator D at threshold follows from ( 35) with (39), which allows one to find all cross-sections using formulas from the Appendix B. For example, the scattering cross-section from (29) becomes at threshold: The scattering cross-section is not the most natural parameter for a driven spaser, because it can become arbitrarily large for small external fields 3 s a t e E << (as all numerical pre-factors in (40) are finite).
The freely generating spaser is described by the relation (37) with zero external field

C.3. Taylor expansion of the denominator near the threshold
The Taylor expansion of denominator D in (35) with the help of chain rule results in: The 1 st term disappears, while the 2 nd term is a complex, but non-singular expression depending on the dispersion of all materials and structure geometry.The outer derivative in the 3 rd term is also non-singular: The dielectric function 1 ε is defined by the expression (1).Its threshold value is: The partial derivatives in (41) Expansion (41) with the derivatives from the above leads to the expression (36) This is the relation (37) with the notations (38) This holds as long as e δ = at zero external field * 3 0 e = .Further from the threshold all quantities cease to be small, and the results of this appendix may serve as qualitative guidelines only.Figure 11 shows the influence of the external field 3 e on the core field amplitude 1 e and on the scattering cross-section sca σ for the externally driven quasi-static core-shell spaser near the threshold; it should be compared with the Fig. 3  in a freely generating spaser 3 0 e → above threshold (i.e., linear dependence of intensity s, see Eq. ( 6)) can also be better seen in Fig. 11(a) than in Fig. 3(b)., is related to the freely generating spaser.Figure 14 is also helpful for the understanding of numerical results for the quadrupolar and spheroidal cases, in particular for the Fig. 7 from the main text.For the linear problem, the absorption in Fig. 14(a) reveals a negative singular "pillar", preceded by a maximum at lower gain values.This can be explained as follows.At low gain levels, plasmonic resonance has high absorption.At the same time, the influence of gain is designed such, as to provide singularity at threshold gain near this resonance.This singularity is positive for scattering, due to radiative losses.Energetically, this implies that absorption singularity is negative, but the absorption enhanced by the plasmonic resonance still shows up at somewhat lower gain values.Saturation bends the negative absorption pillar in such a way, that it becomes positive in a certain region in Figs.14(c) and 14(e).

D.2. Behavior of absorption and extinction and loss compensation
Linear extinction in Fig. 14(b) has both positive and negative singular "pillars", which touch at the threshold value -this behavior is expected for the real part of a complex function near a simple pole (see Eq. ( 29)).These pillars bend in the same direction and merge, which can be better seen in Fig. 14(f).Similarly to Figs. 12(d) and 13(d), Figs.14(d) and 14(f) also illustrate, that that above threshold ext 0 σ = can be achieved on both sides of the resonance.
The condition ext 0 σ = can be considered as a proxy for loss compensation in metamaterials, and our findings are in agreement with those from [16,20,25,41].

D.3. On loss compensation by spasers in metamaterials
Several further comments are in place here.The condition ext 0 σ = , which is equivalent to a condition 3 0 d ′′ = for the overall dipole moment of a spaser used in [20] (see Eqns. ( 28)-( 29)), are considered above for the vacuum background.In real metamaterials (MM) the background usually contains other absorbing materials.This will shift all the resonances and will require overcompensation of losses.Such an overcompensation condition involves all 3 important parameters: available gain L ε , external field 3 e and frequency ω , and defines a 2D surface in this 3D space.It might be difficult to make these parameters, especially the external field and the pumping spatially homogeneous, which might complicate the practical realization.
In realistic compensating arrangement, many small spasers densely embedded into the material should be used.Each spaser will see not the external field alone, but the selfconsistent mean field, which includes the radiation of other spasers as well as scattering from the background (meta)material.In other words, the spasers can be treated as artificial molecules with the number density N .Alternatively, one can describe the whole material macroscopically with the combination of absorbing and lasing dielectric functions, similarly to time domain approach used in [26], or frequency domain description with the saturated gain dielectric function, suggested in the present work.
gain strength in the line center.Gain saturation weakens and broadens the Lorentzian, as quantified by the parameter .(1), relation (2) is a non-linear equation for the core field 1 e .For small values of available gain L ε it has a single solution, while above a certain threshold value L the spaser generation frequency thr ω .At threshold, 0 D = and threshold values for the dipolar quasi-static core-shell spaser are given by the solutions thr thr ( , )

Figure 2 (
Figure 2(a) shows the amplitude of the core field 1 e obtained from the Eq.(2) as a function of frequency ω and external driving (or seeding) field 3 e at threshold gain strength L ), marked by the arrow.Its core field 1 e and generation frequency gen ω are derived in Appendix C.1, Eq.(34).Under the assumptions discussed there, in the quasi-static case, the generation frequency remains constant, intensity s is linear with the available gain L

Figures 2 (
Figures 2(b) and 2(d) show the normalized scattering cross-section 2 sca 2 / a σ π , calculated from the Eq.(29) in Appendix B for threshold and post-threshold values of L ε .The Figs. 2(c)and 2(d) are topologically equivalent.In both cases, there exists a region of parameters with 3 solutions (bistability).For the over-pumping L t h r / 1.5 Appendix D, Fig. 10(a).Further, the dependence of the core field 1 e on the pair of parameters L 3 ( , ) e ε is illustrated by the Figs. 9 and 10(b).

2. 1 . 2 Fig. 3 .
Fig. 3. Influence of saturation on internal field e 1 and scattering cross-section σ sca for the driven quasi-static dipolar core-shell spaser.Dependence of e 1 and σ sca on frequency ω and available gain ε L .External field e 3 =0.006Esat .Other parameters are as in Fig. 2. (a) and (b) show core field amplitudes e 1 , while (c) and (d) show scattering cross-sections.(a) and (c) refer to the linear problem without saturation, while (b) and (d) include saturation.Let us now study the dependence of core field 1 e on the pair . All quantities become finite (Figs.3(b), 3(d)); the core field exactly at threshold is given after the Eq.(5), and the scattering cross-section by Eq. (40) in Appendix C.2. Saturation is important also below threshold.The singular "pillar" near the threshold becomes tilted; as a result, the core field and scattering cross-section at a resonant frequency thr .Both tilted pillars narrow with increasing L ε , and the upper branch is close to the regime of the freely generating spaser (arrow in Fig.3(b)).At a large enough gain level and fixed external field, the core field 1 e and the scattering cross-section demonstrate abrupt changes as a function of frequency.The transition to the linear problem with 3 0 e → is discussed in Appendix D (Fig.11).

Fig. 4 .
Fig. 4. Influence of the external field e 3 on the saturated optical cross-sections for the driven quasi-static dipolar core-shell spaser at threshold gain level ε L = ε thr .Other parameters are as in Fig. 2. Panels (a), (b) and (c) show extinction, scattering and absorption respectively.The curves are labeled in panel (b).Solid curves refer to the linear problem without saturation, s = 0. Curves with nonlinear saturation, s = |e 1 | 2 /E sat 2 are shown for several external field values e 3 listed in (b).For clarity, not all curves are shown in each panel -(a) and (c) include curves at small e 3 , which better illustrate the transition to the linear case without saturation.The dashed ellipse in (a) marks the region where, for low values of the external field e 3 , the extinction cross-section σ ext ≤ 0, which corresponds to a loss compensation.(d) shows all three optical cross-sections for the external field value e 3 = 0.006E sat without saturation (thick curves) and with saturation (thin curves).

Figure 4
and saturated cases.Unsaturated, linear cross-sections with 0 s = are shown by the solid curves in Figs.4(a), 4(b), and 4(c), while the cross-sections for the nonlinear case with saturation of external field 3 e , labeled in the panel (b).
(a)).The boundary of the region ext 0 σ ≤ is often used as a proxy for loss compensation by spasers in metamaterials, which exactly at threshold L on both sides of the resonance (Fig.12(d) in Appendix D), see also Refs.

Fig. 5 .
Fig.5.Analytical (lines) and numerical (symbols) results for a small dipolar core-shell spaser in air, which is shown in Fig.1.Analytical results are normalized to the quasi-static (QS) values ε thr,QS and ω thr,QS , while numerical ones are normalized to the numerical thresholds ε thr,Num and ω thr,Num .The numerical field value in the core center was used.(a) The dependence of the amplitude of the core field ⏐e 1 ⏐ on the external field e 3 for three values of available gain ε L (labeled in the plot) exactly on resonance ω = ω thr .The solid curve for ε L = ε thr is a cut of the surface shown in Fig.2(a) at ω/ω thr = 1, while the solid bistable curve for ε L = 1.5ε thr is a cut of the Fig.2(c), also at ω/ω thr = 1.The parameters marked by a full circle on the upper curve were used to draw a field distribution in the right panel of Fig. 1.(b) The dependence of the amplitude of the core field ⏐e 1 ⏐ on frequency detuning above threshold (ε L = 1.5ε thr ) for several values of external field labeled in the plot.Solid curves correspond to the cuts of the surface shown in Fig. 2(c) at corresponding e 3 /E sat values (curve for e 3 /E sat = 0.0042 is not shown).

Figure 5 (
Figure 5(a) shows the dependence of the core field amplitude 1 e on the external field 3 e , on resonance,

.
h = , outer gain shell 3 50 nm h = , gain material polystyrene (PS) with the dielectric function h The material outside the spaser is air with 4 1 ε = .
field 4 s a t / 0.0124 e E = , as marked by the small filled circle on the upper curve in Fig. 6(a).As a measure of field in the gain medium we chose its maximum value at the outer side of the Ag-shell at near 45 degrees from the incident wave direction in the polarization plane.The Figs. 6(a) and 6(b) show the dependence of this maximum field on the external field 4 e like in Figs.2(a) and 2(b).

Fig. 6 .
Fig. 6.Large shell quadrupolar spaser driven by the external field, similar to Fig. 5. Symbols represent numerical results, curves are guides for the eye.The inset in (b) shows the geometry of the structure.Core radius a 1 = 50 nm, Ag shell thickness h 2 = 4.3 nm, outer gain shell h 3 = 50 nm; material parameters are described in the text.Results are normalized to the numerical ε thr, Num and ω thr, Num values given in the text.(a) The dependence of the maximum field in the gain material E max on the external field e 4 for three values of available gain ε L (labeled in the plot) near the resonance.This maximum is achieved in the outer gain shell (subscript "3"), near the Ag surface, in the plane of incidence defined by the k-vector and external field vector e 4 , at an angle 45° with respect to the former.The inset shows the field distribution for a gain level of ε L = 1.5ε thr and an external field e 4 /E sat = 0.0124, indicated by the filled circle at the upper curve.The shape of the filed distribution is similar for other parameters.(b) The dependence of the maximum field E max on the frequency detuning above threshold (ε L = 1.5ε thr ) for three values of external field e 4 labeled in the plot.
the Figs.9(a) and 9(c) in Appendix D. For small external fields,

Fig. 7 .
Fig. 7. Dependence of cross-sections of the large shell quadrupolar spaser on the external driving field e 4 in the presence of saturation.Numerical results are given near the generation frequency, ω = 1.0001ω thr .Other parameters are as in Fig. 6.(a) ε L = 0.5ε thr , (b) ε L = ε thr , (c) and (d) ε L = 1.5ε thr , with increasing and decreasing external field respectively, to show different branches of solutions.The dashed rectangle indicates the region where the structure can work as an amplifier, as discussed in the text.Note the difference in ordinate scale between all panels.The behavior can be better understood by comparing with Fig. 2(d) for the scattering and Figs.13(a), (b) for absorption and extinction.

Figure 6 (
Figure 6(b) shows the spectral dependence of the maximum field in the gain material.It is similar to the Fig. 5(b), and likewise similar to the cuts of Fig. 2(c) by the planes 4 s a t / e E const =, as labeled in the plot.As the frequency is gradually changed, one sees abrupt changes in the gain field at the points where the upper "spaser-like" solution appears or disappears.
. The incident field is polarized along the long axis.This spheroid is immersed into a spherical gain material (radius 2

Fig. 8 .
Fig. 8. Spheroidal dipolar spaser driven by the external field e 3 , numerical results.a) The dependence of the maximum field in the gain material (near the tip of the Ag spheroid) on the available gain ε L at the resonance ω = ω thr , for the external field e 3 = 0.001 E sat .The upper curve corresponds to the spherical gain material with the radius a 2 = 4c = 80 nm, while the lower one to a gain shell of constant thickness h 2 = 10 nm along all axes of the Ag spheroid.Long semi-axis c = 20 nm in both cases.The results are normalized to the corresponding (somewhat different) numerical ε thr and ω thr values.The insets illustrate the geometries of both structures.Further details about the geometries, material parameters and threshold values are given in the text.The two panels above the main plot show the normalized field amplitude (⏐(E)⏐/E sat ) in the plane of incidence defined by the k-vector and external field vector e 3 , for the gain values ε L = 3ε thr , indicated by the small filled circle and ellipse on the corresponding curves.The shape of the field distribution is similar for other parameter values.The dashed circle indicates the region studied in further detail in Figs.(b) and (c).(b) The region where (for the a 2 = 80 nm case) the high-field branch of solution ceases to exist at ω = ω thr due to a drift in spaser generation frequency with the available gain, ω gen (ε L ).(c) The numerical cut of the solution surface for ε L = 4.7135ε thr , near the termination point in Fig. (b).The arrow shows the shift of the spaser generation frequency for this level of gain.The crosses in Figs.(b) and (c) indicate the same solution point.

3 10 − to 2 10−
(b) and 8(a).Let us gradually increase the available gain level L remain on the thin "tilted pillar" of the upper spaser-like solution branch, while if gen enough values of available gain L ε , the upper branch solution should cease to exist.The numerically calculated field in the gain material for a certain maximum, then gradually decrease, and finally abruptly jump to its value in the lower branch solution.This effect can be seen in Fig.8(b), which is a blow-up of the region marked by the dashed circle in Fig.8(a).At first the slope of the max of solutions disappears.This is the expected behavior, if the thin tilted pillar of solutions in Figs.3(b) (or more pronounced in 11(a)) deviates from the exact thr ω ω = direction.To quantify this deviation, Fig. 8(c) shows the numerical results for the upper branch of solution frequency ω is varied.The central, spaser generation frequency is red-shifted to a value gen L 10 − ⋅ may seem small, but it is computationally significant, as the threshold frequency thr ω was determined numerically with the relative accuracy of about 6 1.3 10 − ⋅ .The relative red shift of frequency discussed in [50] is of the order of 3 .Smaller values in present work can be related to differences in the systems and numbers involved, as well as the close match between the plasmon and atomic resonances in our case.Note, that the real part in (1), in the center of the atomic line, available gain L ε , and it is the real part G ε ′ , what largely determines the generation frequency [17].Hence, the shift is of the order of 4 10 − only.
denominator in the last expression in ( previously.For comparison with other works, we recall standard definitions of population (longitudinal 1 T ) expression results in L ~0.116 ε the typical values of the field in the gain material are in the range of gain sat ẽ E (see Figs. 5, 6, 8, etc. dielectric function (1) of the core depends on s: on three main control parameters: excitation frequency ω , available gain L ε and external field 3

e
ideal resonance between the external field, atomic line and plasmon (and therefore the spaser), 36) further simplifies (see the derivation in Appendix C.3 below) to: , over-pumping parameter L δ , and non-resonant parameter C: with(34); the latter is more general, as it applies also to the detuning between the atomic and nanostructure resonances.Additional functional dependences are discussed at the end of Appendix C.3.

Figure 9 3 0
Figure 9 shows the dependence of the core field amplitude 1 e and the scattering crosssection sca σ on the available gain L ε and the external field 3 e .Figures 9(a) and 9(b) are drawn exactly at resonance thr ω ω =.One can see the freely generating spaser line (i.e.,

Figure 10
is a projection of the Figs.2(c)and 2(d).Figure 10(b) shows the L 3 ( , ) e ε plane, and is a projection of the Figs.9(c),(d); it represents a typical cusp catastrophe [57].On resonance thr ω ω = the bistability region extends up to the axis 3 0 e = .

3 e
from the main text.As the external field decreases, the tilted resonant 1 e -pillar in Fig.11(a) becomes narrower in frequency ( 3(b) and 3(d)), and the upper branch further approaches the freely generating spaser regime (arrow in Fig. 11(b)).With vanishing external field 3 0 e → , this pillar degenerates into a freely generating spaser line.Simultaneously, the "scattering pillar" becomes less tilted (Fig. 11(b)).This illustrates the continuous transition to the scattering singularity in linear case (Figs.3(c)), with vanishing external field 3 0 e → .The square root dependence of the core field on the available gain L ε

Fig. 11 .
Fig. 11.Dependence of core field amplitude |e 1 | and scattering cross-section σ sca on frequency ω and available gain ε L for the driven core-shell spaser, as in Figs.3(b) and 3(d), but for the trice smaller external field e 3 = 0.002E sat .(a) core field, arrow indicates the limiting case of a freely generating spaser, (b) scattering cross section.

Figure 12 and
Figure 12 and Fig. 13 show the behavior of the absorption (a), (c), and extinction (b), (d) cross-sections.Figure 12 further illustrates the changes in spectral behavior of cross-sections with the increasing external field, studied in Fig. 4. It shows the behavior at (a), (b), and above (c), (d) threshold, complementing Figs.2(b) and 2(d), where scattering is shown.In particular, at (or below) threshold gain levels ext 0 σ = can be observed only at the blue wing of the resonance

Figure 14 compares
the saturated and unsaturated absorption (a), (c), (e), and extinction (b), (d), (f), cross-sections in L ( , ) ω ε coordinates, similar to the Figs.3(c) and 3(d) for the scattering in the main text.Figures 14(c) and 14(f) use extended and reversed scale for the available gain L ε , to better illustrate the topology of the absorption and extinction surfaces.It is important to remember the attribution of different branches with respect to the core field.The upper branch, with the higher field and scattering (see Figs. 2(c) and 2(d) and 3(b) and 3(d))

Fig. 14 .
Fig. 14.Influence of saturation on the absorption and extinction cross-sections for a quasistatic dipolar core-shell spaser driven by the external field e 3 = 0.006E sat .Dependence on frequency ω and available gain ε L .Parameters are as in Fig. 3. (a), (c) and (e) show absorption σ abs , while (b), (d) and (f) show extinction cross section σ ext .(a) and (b) refer to the linear problem without saturation (s = 0), while (c)-(f) include gain saturation.(e) and (f) illustrate the behavior at higher available gain levels and have reverse orientation of the ε L axis, which better illustrates the global shape of the surfaces.
In this case, the polarizabilities of the spasers "s" and the neutral background agents "b" enter the Lorentz-Lorenz formula for the dielectric function ε together: equality we utilized the fact that the dielectric function b ε of the virgin background without the spasers obeys the combined dielectric function ε of the spasers + MM composite is found from here, loss compensation corresponds to the condition 0 the properties of the background MM being compensated, b ε , and on the concentration of spasers s lies exactly at the Lorentzian frequency, with was obtained: