## Abstract

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

## 1. Introduction

Surface plasmon amplification by stimulated emission (spaser) has been extensively discussed in the literature over the last decade theoretically [1–5 ] and experimentally [6, 7 ]. Similar concepts were discussed for nano-lasers [8–12 ] and loss compensation in metamaterials [13–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 cross-sections, 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 multi-level rate equations coupled with the quasi-classical description of local polarizations [26–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 cross-sections 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.

## 2. Quasi-static core-shell dipolar spaser

#### 2.1 Analytical results

### 2.1.1 Core 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 ${\epsilon}_{\text{G}}={{\epsilon}^{\prime}}_{\text{G}}+i{{\epsilon}^{\u2033}}_{\text{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.

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 ${e}_{1}$ inside the core can be derived from electrostatic equations (Appendix B), and is related to the external field ${e}_{3}$ in the following way:

*D*is the determinant of the coefficients matrix for the linear equations originating from the boundary conditions at two interfaces. With a saturating ${\epsilon}_{1}={\epsilon}_{\text{G}}({e}_{1})$, as defined by Eq. (1), relation (2) is a non-linear equation for the core field ${e}_{1}$. For small values of available gain ${\epsilon}_{\text{L}}$ it has a single solution, while above a certain threshold value ${\epsilon}_{\text{L}}>{\epsilon}_{\text{thr}}$, a bistability appears near the spaser generation frequency ${\omega}_{\text{thr}}$. At threshold, $D=0$ and threshold values for the dipolar quasi-static core-shell spaser are given by the solutions $({\omega}_{\text{thr}},\text{\hspace{0.17em}}{\epsilon}_{\text{thr}})$ of the complex-valued equation:

The following parameters were used in the calculations: Shell material is silver with ${\epsilon}_{2}={\epsilon}_{\infty}-\frac{{\omega}_{\text{p}}^{2}}{\omega (\omega +i\gamma )}$, where ${\epsilon}_{\infty}=4.9,\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\omega}_{\text{p}}=9.5\text{\hspace{0.17em}}\text{eV,}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\gamma \text{=0}\text{.05}\text{\hspace{0.17em}}\text{eV}$. 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 $-{{\epsilon}^{\u2033}}_{2}/{{\epsilon}^{\prime}}_{2}$ ratio [17], around $\omega \approx 2.479\text{\hspace{0.17em}}\text{eV}\text{\hspace{0.17em}}\text{\hspace{0.17em}}(\lambda \approx 500\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{nm})$. Core and shell sizes and the parameters of the gain Lorentzian (1) were chosen to make the resonance close to this value, namely: core radius ${a}_{1}=10\text{\hspace{0.17em}}\text{nm}$, shell thickness ${h}_{2}=2.7\text{\hspace{0.17em}}\text{nm}$, core material polystyrene (PS) with background dielectric function ${\epsilon}_{\text{h}}=2.6$, Lorentzian center frequency and FWHM: ${\omega}_{\text{L}}=2.5\text{\hspace{0.17em}}\text{eV}\text{\hspace{0.17em}}\text{\hspace{0.17em}}({\lambda}_{\text{L}}=495.94\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{nm})\text{,}\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\gamma}_{\text{L}}\text{=0}\text{.25}\text{\hspace{0.17em}}\text{eV}$. Such a broad emission line is typical for organic dyes. The outer material is air with ${\epsilon}_{3}=1$. This results in the quasi-static threshold values calculated from the condition (3):

*ω*and external driving (or seeding) field ${e}_{3}$ at threshold gain strength ${\epsilon}_{\text{L}}={\epsilon}_{\text{thr}}$. It is convenient to normalize the control parameters $\omega /{\omega}_{\mathrm{thr}}$, ${e}_{3}/{E}_{\text{sat}}$ and ${\epsilon}_{\text{L}}/{\epsilon}_{\text{thr}}$, because the functional dependences remain similar between different geometries, while material properties and spaser geometry largely influence the numerical pre-factors and threshold values.

The core field is single valued, and the amplitude monotonously increases with the external field . Post-threshold behavior is shown in Fig. 2(c) for ${\epsilon}_{\text{L}}/{\epsilon}_{\text{thr}}=1.5$, where three solutions exist for low values of external field and detuning. Exactly at resonance $\omega ={\omega}_{\text{thr}}$, 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):

The freely generating spaser without an external field, ${e}_{3}=0$, corresponds to the tip of the conical surface in Fig. 2(c), marked by the arrow. Its core field ${e}_{1}$ and generation frequency ${\omega}_{\text{gen}}$ are derived in Appendix C.1, Eq. (34). Under the assumptions discussed there, in the quasi-static case, the generation frequency remains constant, ${\omega}_{\text{gen}}({\epsilon}_{\text{L}})={\omega}_{\text{thr}}$, while the intensity *s* is linear with the available gain ${\epsilon}_{\text{L}}$ above threshold, in agreement with findings in [20, 22, 23
]:

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 analysiscarried 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.

### 2.1.2 Core field and optical cross-sections as a function of frequency and available gain

Let us now study the dependence of core field ${e}_{1}$ on the pair $(\omega ,\text{\hspace{0.17em}}{\epsilon}_{\text{L}})$ for different external fields ${e}_{3}$. Such a situation often occurs in experiments and is studied in Fig. 3
. Without saturation $(s=0)$, there exists a gain level ${\epsilon}_{\text{thr}}$ and a generation (threshold) frequency ${\omega}_{\text{thr}}$, where all quantities become singular (Figs. 3(a), 3(c)). 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–39
]. A similar effect was observed for particles where a gain in the form of ${\epsilon}^{\u2033}<0$ was introduced into the dielectric components of an anisotropic hyperbolic permittivity [40]. In contrast, saturation $(s={\left|{e}_{1}\right|}^{2}/{E}_{\text{sat}}^{2})$ leads to the following qualitative changes. 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 $\omega ={\omega}_{\text{thr}}$ monotonously increase with the available gain${\epsilon}_{\text{L}}$ and smoothly pass the threshold value ${\epsilon}_{\text{L}}={\epsilon}_{\text{thr}}$ without any singularity. Both tilted pillars narrow with increasing ${\epsilon}_{\text{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 ${e}_{1}$ and the scattering cross-section demonstrate abrupt changes as a function of frequency. The transition to the linear problem with ${e}_{3}\to 0$ is discussed in Appendix D (Fig. 11).

The behavior of other cross-sections is more complex. Figure 4 compares all 3 optical cross-sections exactly at threshold gain ${\epsilon}_{\text{L}}={\epsilon}_{\text{thr}}$ for the unsaturated and saturated cases. Unsaturated, linear cross-sections with $s=0$ 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 $s={\left|{e}_{1}\right|}^{2}/{E}_{\text{sat}}^{2}$ depend on the values of external field ${e}_{3}$, labeled in the panel (b). 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 cross-section 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 ${e}_{3}$ increases, and the extinction remains positive everywhere to the red side of the resonance. At the blue side, ${\sigma}_{\text{ext}}\le 0$ in the region indicated by the dashed ellipse; this region depends on the external field ${e}_{3}$. For small fields ${e}_{3}$, weakly saturated gain can overcompensate the losses, while for larger ${e}_{3}$, the gain saturates, and ${\sigma}_{\text{ext}}>0$ everywhere (dash-dotted curve in Fig. 4(a)). The boundary of the region ${\sigma}_{\text{ext}}\le 0$ is often used as a proxy for loss compensation by spasers in metamaterials, which exactly at threshold ${\epsilon}_{\text{L}}={\epsilon}_{\text{thr}}$ requires $\omega >{\omega}_{\text{thr}}$. Above threshold, for ${\epsilon}_{\text{L}}>{\epsilon}_{\text{thr}}$, this restriction is lifted, and regions with ${\sigma}_{\text{ext}}<0$ exist on both sides of the resonance (Fig. 12(d) in Appendix D), see also Refs. [16, 20, 25, 41 ]. Figure 4 uses ${\epsilon}_{\text{L}}={\epsilon}_{\text{thr}}$, 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 .

#### 2.2 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–44 ]. Equating the corresponding denominator to zero, one obtains the threshold values with retardation:

With retardation, the field ${e}_{1}$ 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 ${\epsilon}_{\text{L}}/{\epsilon}_{\begin{array}{l}\text{thr}\\ \text{Mie}\end{array}}$ and $\omega /{\omega}_{\begin{array}{l}\text{thr}\\ \text{Mie}\end{array}}$ will be similar to the non-retarded case with respect to ${\epsilon}_{\text{L}}/{\epsilon}_{\begin{array}{l}\text{thr}\\ \text{QS}\end{array}},\text{\hspace{0.17em}}\text{\hspace{0.17em}}\omega /{\omega}_{\begin{array}{l}\text{thr}\\ \text{QS}\end{array}}$. 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.

Figure 5(a) shows the dependence of the core field amplitude $\left|{e}_{1}\right|$ on the external field ${e}_{3}$, on resonance, $\omega ={\omega}_{\text{thr}}$, at 3 gain levels ${\epsilon}_{\text{L}}/{\epsilon}_{\text{thr}}=0.5,\text{\hspace{0.17em}}\text{\hspace{0.17em}}1,\text{\hspace{0.17em}}\text{\hspace{0.17em}}1.5$ labeled in the plot. The curve with ${\epsilon}_{\text{L}}/{\epsilon}_{\text{thr}}=1$ corresponds to the cut of the Fig. 2(a) along the line $\omega /{\omega}_{\text{thr}}=1$, while the bistable curve with ${\epsilon}_{\text{L}}/{\epsilon}_{\text{thr}}=1.5$ corresponds to the similar cut of the Fig. 2(c). As expected from the approximation (5), the relation between the fields ${e}_{1}$ and ${e}_{3}$ is close to a cubic function, which has 1 solution up to the threshold and 3 solutions above it. Solid curves in Fig. 5(b) show the influence of frequency detuning for several values of the incident field. They correspond to cuts of the Fig. 2(c) at fixed ${e}_{3}/{E}_{\text{sat}}$ values.

Different solution branches above threshold can be explored numerically, if one enters the bistable region from opposite sides, using the previous solutions as starting points for the iterations. For example, for ${\epsilon}_{\text{L}}=1.5{\epsilon}_{\text{thr}}$, if the external field is gradually increased from ${e}_{3}=0$, the lower branch in Fig. 5(a) can be calculated up to ${e}_{3}=0.0042{E}_{\text{sat}}$. If the external field is gradually decreased from the value ${e}_{3}=0.006{E}_{\text{sat}}$ (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 ${e}_{3}$, 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.

## 3. 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 ${a}_{1}=50\text{\hspace{0.17em}}\text{nm}$, Ag shell thickness ${h}_{2}=4.3\text{\hspace{0.17em}}\text{nm}$, outer gain shell ${h}_{3}=50\text{\hspace{0.17em}}\text{nm}$, gain material polystyrene (PS) with the dielectric function ${\epsilon}_{\text{h}}=2.6$, Lorentzian center ${\omega}_{\text{L}}=1.6369\text{\hspace{0.17em}}\text{eV}$ $({\lambda}_{\text{L}}=757.42\text{\hspace{0.17em}}\text{nm)}$ and FWHM ${\gamma}_{\text{L}}\text{=0}\text{.25}\text{\hspace{0.17em}}\text{eV}$. The material outside the spaser is air with ${\epsilon}_{4}=1$. For the sake of generality, we assume a spatially homogeneous available gain, i.e., ${\epsilon}_{\text{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–44 ] lies exactly at the Lorentzian frequency, with

*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.

Similarly to Fig. 5, the system demonstrates bistability above the threshold for the ${\epsilon}_{\text{L}}=1.5{\epsilon}_{\text{thr}}$ curve in the Fig. 6(a), and all curves in Fig. 6(b). The lower branch disappears at external field ${e}_{4}\approx 0.0124{E}_{\text{sat}}$, 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 ${e}_{4}=0$, but disappears around ${e}_{4}=0.00138{E}_{\text{sat}}$. This is due to slight detuning from the (numerical) threshold frequency ${\omega}_{\text{thr,Num}}$ from (10), which translates into the relative detuning of $\omega /{\omega}_{\text{thr,Num}}=1.0001$, see also the Figs. 9(a) and 9(c) in Appendix D. For small external fields, ${e}_{4}<0.00138{E}_{\text{sat}}$, the lower branch is the only existing solution. Seeding of a spaser by a week, 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, ${e}_{4}>0.00138{E}_{\text{sat}}$, we expect the upper, spaser-like, solution to be realized in practice.Figure 6(b) shows the spectral dependence of the maximum field in the gainmaterial. It is similar to the Fig. 5(b), and likewise similar to the cuts of Fig. 2(c) by the planes ${e}_{4}/{E}_{\text{sat}}=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.

Figure 7 shows the dependence of all cross-sections on the external field at three values of available gain: ${\epsilon}_{\text{L}}=0.5,\text{\hspace{0.17em}}\text{\hspace{0.17em}}1,\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{and}\text{\hspace{0.17em}}\text{\hspace{0.17em}}1.5{\epsilon}_{\text{thr}}$. 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 ${\omega}_{\text{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 ${e}_{4}<0.00138{E}_{\text{sat}}$ (as in Fig. 6(a)), the only existing solution has both ${\sigma}_{\text{abs}}<0$ and ${\sigma}_{\text{ext}}<0$ (see are marked dashed in Figs. 7(c) and 7(d)). For these external fields ${e}_{4}$, the considered CS structure is expected to operate as a stable amplifier.

## 4. Prolate spheroid dipolar spaser

#### 4.1 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 semi-axis is $c=20\text{\hspace{0.17em}}\text{nm}$, and of the shorter ones is $a=b=5.56\text{\hspace{0.17em}}\text{nm}$, which corresponds to an aspect ratio $c/a=3.6$. The incident field is polarized along the long axis. This spheroid is immersed into a spherical gain material (radius ${a}_{2}=4c=80\text{\hspace{0.17em}}\text{nm}$) with ${\epsilon}_{2}$ defined by the Lorentzian (1) with the central frequency ${\omega}_{\text{L}}=1.6123\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{eV}\text{\hspace{0.17em}}({\lambda}_{\text{L}}=768.99\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{nm)}$ and FWHM ${\gamma}_{\text{L}}\text{=0}\text{.25}\text{\hspace{0.17em}}\text{eV}$. 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 ] is:

The upper curve in Fig. 8(a) shows the numerical dependence of the maximum field in the gain material near the tip of the Ag spheroid on the available gain ${\epsilon}_{\text{L}}$. These results are at the resonance $\omega ={\omega}_{\text{thr}}$, for the external field ${e}_{3}/{E}_{\text{sat}}=0.001$. Because this value of the external field is low, the results closely emulate the field in the freely generating spaser. The upper curve is qualitatively similar to the cut of the Fig. 3(b) by the plane $\omega ={\omega}_{\text{thr}}$, and the corresponding “tilted pillar” is narrow for ${e}_{3}/{E}_{\text{sat}}=0.001$. The maximum gain field ${E}_{\mathrm{max}}/{E}_{\text{sat}}$ is significantly larger than the almost constant core field for quasi-static spaser inFigure 3(b). Strong field enhancement and localization near the tip of spheroid can be seen in the left Fig. of the upper panel for ${\epsilon}_{\text{L}}=3{\epsilon}_{\text{thr}}$ (compare with right panel in Fig. 1).

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 ${a}_{2}=80\text{\hspace{0.17em}}\text{nm}$ 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 ${h}_{2}=10\text{\hspace{0.17em}}\text{nm}$ along all semi-axes. 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 $c=20\text{\hspace{0.17em}}\text{nm}$, and shorter ones $a=b=4.9\text{\hspace{0.17em}}\text{nm}$, which corresponds to the aspect ratio $c/a=4.08$. This results in:

An interesting observation is the difference in the behavior of the maximum field in two cases. While for the thin gain shell (red curve) ${E}_{\text{max}}/{E}_{\text{sat}}\propto \sqrt{{\epsilon}_{\text{L}}-{\epsilon}_{\text{thr}}}$ 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.

#### 4.2 Changes in the generation frequency with the available gain

The quasi-static spaser generation frequency does not depend on the available gain level ${\epsilon}_{\text{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 ${\epsilon}_{\text{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 ${e}_{\text{ext}}=0$ 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 ${e}_{4}/{E}_{\text{sat}}<<0.0015$ in the Fig. 6(b)). An even better numerical recipe is provided by the Figs. 3(b) and 8(a) . Let us gradually increase the available gain level ${\epsilon}_{\text{L}}$, while staying on resonance $\omega ={\omega}_{\text{thr}}$ at small values of ${e}_{\text{ext}}/{E}_{\text{sat}}<<1$. If the generation frequency does not change, ${\omega}_{\text{gen}}({\epsilon}_{\text{L}})={\omega}_{\text{thr}}$, we will remain on the thin “tilted pillar” of the upper spaser-like solution branch, while if ${\omega}_{\text{gen}}({\epsilon}_{\text{L}})\ne {\omega}_{\text{thr}}$, for large enough values of available gain ${\epsilon}_{\text{L}}$, the upper branch solution should cease to exist. The numerically calculated field in the gain material for a fixed $\omega ={\omega}_{\text{thr}}\ne {\omega}_{\text{gen}}({\epsilon}_{\text{L}})$ should first reach 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 ${E}_{\mathrm{max}}(\omega ={\omega}_{\text{thr}},\text{\hspace{0.17em}}{\epsilon}_{\text{L}})$ starts to decrease, reaches a maximum and finally, around the value ${\epsilon}_{\text{L}}=4.7135{\epsilon}_{\text{thr}}$ the upper branch 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 $\omega ={\omega}_{\text{thr}}$ direction. To quantify this deviation, Fig. 8(c) shows the numerical results for the upper branch of solution ${E}_{\mathrm{max}}(\omega ,\text{\hspace{0.17em}}{\epsilon}_{\text{L}}=4.7135{\epsilon}_{\text{thr}})$, when the frequency *ω* is varied. The central, spaser generation frequency is red-shifted to a value ${\omega}_{\text{gen}}({\epsilon}_{\text{L}}=4.7135{\epsilon}_{\text{thr}})=0.999884{\omega}_{\text{thr}}$, indicated by an arrow. This relative shift of $1.12\cdot {10}^{-4}$ may seem small, but it is computationally significant, as the threshold frequency ${\omega}_{\text{thr}}$ was determined numerically with the relative accuracy of about $1.3\cdot {10}^{-6}$. The relative red shift of frequency discussed in [50] is of the order of ${10}^{-3}$ to ${10}^{-2}$. 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, ${{\epsilon}^{\prime}}_{\text{G}}({\omega}_{\text{L}})={\epsilon}_{\text{h}}$ remains independent on the available gain ${\epsilon}_{\text{L}}$, and it is the real part ${{\epsilon}^{\prime}}_{\text{G}}$, what largely determines the generation frequency [17]. Hence, the shift is of the order of ${10}^{-4}$ only.

## 5. 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 ${\epsilon}_{\text{G}}({\left|E(r)\right|}^{2})$ 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 cross-sections 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. In 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 ${10}^{-4}$ 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 ${N}_{i}$ 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 ${N}_{3,1}\approx 0$, ${N}_{0}+{N}_{2}\approx {N}_{\text{tot}}$ (total concentration in cm^{−3}), and the rate equations reduce to a single equation for the population inversion of lasing transition 2→1, ${N}_{2}-{N}_{1}\approx {N}_{2}$

*W*. The approximation (15) holds for excitation rates smaller than depopulation rates of levels 3 and 1: ${W}_{\text{p}},{W}_{\text{s}}<<{\gamma}_{32,}{\gamma}_{10}~{10}^{12}{\text{s}}^{\text{-1}}$, which translates into intensities:

*pumping*saturation intensity, similar to the one estimated in (19) for the emission line: ${I}_{\text{p,sat}}=\frac{{\gamma}_{\text{N}}\hslash {\omega}_{03}}{{\sigma}_{03}}~6.6\cdot {10}^{5}{\text{W/cm}}^{2}$.

Let us now express the saturation field via the (absolute value of) matrix element for the dipole moment ${\mu}_{12}$. We recall that the emission cross-section for the Lorentzian line with the FWHM ${\gamma}_{\text{L}}$ is [19] Eq. (7.5.62), or [51], Eq. (2.4.18):

*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

The amplitude of the Lorentzian gain profile ${\epsilon}_{L}$ in the expression (1) is proportional to the atomic polarizability $\alpha $ and to the equilibrium population inversion ${N}_{2}(0)$, in CGS units ${\epsilon}_{\text{L}}=4\pi \alpha {N}_{2}$. 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 $4\pi {\epsilon}_{0}$ in SI units):

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 ${\left(\frac{{\epsilon}_{\text{h}}+2}{3}\right)}^{2}$ to the gain Lorentzian and to the saturation term *s* in the denominator of expression (1). This is equivalent to the following rescaling of ${\epsilon}_{\text{L}}$ and saturation intensity or field indicated by the additional subscript “LL”:

*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 Fig. 8. Using Purcell-enhanced saturation field estimated after the Eq. (23) and the modal volume $V\approx {10}^{4}\text{\hspace{0.17em}}{\text{nm}}^{3}$, Eq. (25) gives plasmon numbers ${n}_{\text{pl}}~1$ near the threshold. This is consistent with the predictions of quantum-mechanical 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 is${E}_{3}$, 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].

## Appendix C Spaser threshold, free generation, and near-threshold response

#### C.1. Spaser threshold and generation without an external field

Let us first analyze the freely generating spaser without an external field, ${e}_{3}=0$. 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 $s={\left|{e}_{1}\right|}^{2}/{E}_{\text{sat}}^{2}$, because the saturated dielectric function (1) of the core depends on *s*: ${\epsilon}_{1}={\epsilon}_{\text{G}}(\omega ,{\epsilon}_{\text{L}},s)$. The phase of the core field is of no importance for the generation analysis. The consistency condition $D=0$ can be resolved with respect to ${\epsilon}_{1}$:

*s*and the generation frequency ${\omega}_{\text{gen}}$, considered as a function of all parameters and available gain ${\epsilon}_{\text{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 $s\to 0$. This defines the threshold frequency ${\omega}_{\text{thr}}$ and available gain level ${\epsilon}_{\text{L}}={\epsilon}_{\text{thr}}$, such thatAbove the threshold, the Eq. (30) should hold for $s\ne 0$, so that at the generation frequency ${\omega}_{\text{gen}}$:In several important cases, we can easily find a physical solution of this equation. If we choose ${\omega}_{\text{gen}}={\omega}_{\text{thr}}$, the complex number on the r.h.s. in (31)-(32) remains unchanged. If in addition there exists a real normalized intensity $s={\left|{e}_{1}\right|}^{2}/{E}_{\text{sat}}^{2}$ such, that the complex number on the l.h.s. does not change as well, the Eq. (32) is fulfilled automatically. This requires:

*s*. For Lorentzian gain (1) this results in:

*s*together with ${\omega}_{\text{gen}}={\omega}_{\text{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 ${\epsilon}_{\text{L}}$. The slope of this linear dependence is defined by (34). Note, that the threshold parameters ${\omega}_{\text{thr}}$ and ${\epsilon}_{\text{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 ${\epsilon}_{\text{G}}(\omega ,{\epsilon}_{\text{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 ${\omega}_{\text{gen}}({\epsilon}_{\text{L}}>{\epsilon}_{\text{thr}})\ne {\omega}_{\text{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 ${e}_{3}\ne 0$, all quantities remain finite even when the spaser is driven at its native generation frequency $\omega ={\omega}_{\text{thr}}$. The nonlinear equation (2) defines the dependence of the core field ${e}_{1}={e}_{1}({\epsilon}_{\text{L}},\omega ,{e}_{3})$ and other quantities on three main control parameters: excitation frequency $\omega $, available gain ${\epsilon}_{\text{L}}$ and external field ${e}_{3}$. There exists a phase difference between ${e}_{3}$ and ${e}_{1}$, or other quantities. As ${e}_{1}$ enters all expressions non-linearly, it is usually more convenient to treat ${e}_{1}$ as real and calculate the phases “backwards”.

To obtain general dependences, we rewrite (2) with respect to its denominator *D*:

*D*on the l.h.s. Here, all dielectric functions ${\epsilon}_{1,2,3}$ should be expanded in frequency, while ${\epsilon}_{1}$, defined by (1), additionally depends on the available gain ${\epsilon}_{\text{L}}$ and the normalized intensity $s$. The expressions for all the derivatives taken at threshold $\omega ={\omega}_{\text{thr}},\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\epsilon}_{\text{L}}={\epsilon}_{\text{thr}},\text{\hspace{0.17em}}\text{\hspace{0.17em}}s=0$, are listed in Appendix C.3. The equality (35) is then approximated by:

*C*:

*C*is a complex number, and a phase difference exists between ${e}_{1}$ and ${e}_{3}$. Exactly at the threshold ${\delta}_{\text{L}}=0$, and the core field from (37) becomes:

*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 freely generating spaser is described by the relation (37) with zero external field ${e}_{3}=0$, which for an ideal resonance results in ${e}_{1}^{*}=\sqrt{{\delta}_{\text{L}}}$, in consistence 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.

#### 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:

^{1st}term disappears, while the

^{2nd}term is a complex, but non-singular expression depending on the dispersion of all materials and structure geometry. The outer derivative in the

^{3rd}term is also non-singular:

^{1st}term there disappears, while $\partial D/\partial {\epsilon}_{1}$, defined by (42), can be simplified if the term ${\epsilon}_{2}+2{\epsilon}_{3}$ is taken from the threshold condition (30). This results in:

## Appendix D Additional Figs. for the driven quasi-static spaser

#### D.1. Auxiliary Figs. for the core field and scattering cross-section

Figure 9 shows the dependence of the core field amplitude $\left|{e}_{1}\right|$ and the scattering cross-section ${\sigma}_{\text{sca}}$ on the available gain ${\epsilon}_{\text{L}}$ and the external field ${e}_{3}$. Figures 9(a) and 9(b) are drawn exactly at resonance $\omega ={\omega}_{\text{thr}}$. One can see the freely generating spaser line (i.e., ${e}_{3}=0$), indicated by the arrow in the Fig. 9(a), which corresponds to the square root behavior (34). The scattering cross-section diverges exactly at threshold for small external fields in accordance with the Eq. (40). Small red detuning from the resonance frequency $\omega =0.999{\omega}_{\text{thr}}$ is shown in Figs 9(c) and 9(d). It removes the free spaser curve (the bistability region does not extend to zero of external field ${e}_{3}=0$), and makes all cross-sections finite everywhere. The behavior of the scattering is similar for small blue detuning $\omega =1.001{\omega}_{\text{thr}}$ (not shown).

Above threshold, there exists a region of bistability in the parameter space, where 3 solutions for the core field ${e}_{1}$ in Eq. (2) exist. This region is shown in the Fig 10 . Figure 10(a) is drawn in the $(\omega ,\text{\hspace{0.17em}}{e}_{3})$ plane, and is a projection of the Figs. 2(c)and 2(d). Figure 10(b) shows the $({\epsilon}_{\text{L}},\text{\hspace{0.17em}}{e}_{3})$ plane, and is a projection of the Figs. 9(c),(d); it represents a typical cusp catastrophe [57]. On resonance $\omega ={\omega}_{\text{thr}}$ the bistability region extends up to the axis ${e}_{3}=0$.

Figure 11
shows the influence of the external field ${e}_{3}$ on the core field amplitude $\left|{e}_{1}\right|$ and on the scattering cross-section ${\sigma}_{\text{sca}}$ for the externally driven quasi-static core-shell spaser near the threshold; it should be compared with the Fig. 3 from the main text. As the external field ${e}_{3}$ decreases, the tilted resonant ${e}_{1}$-pillar in Fig. 11(a) becomes narrower in frequency (${e}_{3}=0.002{E}_{\text{sat}}$ in Fig. 11 vs. ${e}_{3}=0.006{E}_{\text{sat}}$ in Figs. 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 ${e}_{3}\to 0$, 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 ${e}_{3}\to 0$. The square root dependence of the core field on the available gain ${\epsilon}_{\text{L}}$ in a freely generating spaser ${e}_{3}\to 0$ 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).

#### D.2. Behavior of absorption and extinction and loss compensation

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 ${\sigma}_{\text{ext}}=0$ can be observed only at the blue wing of the resonance $\omega >{\omega}_{\text{thr}}$ (Fig. 12(b)), while above threshold (Fig. 12(d)), ${\sigma}_{\text{ext}}=0$ can be observed at both sides of the resonance. Figure 13 shows the dependence on available gain at resonance $\omega ={\omega}_{\text{thr}}$ (a), (b), and for a small blue detuning $\omega =1.001{\omega}_{\text{thr}}$, (c), (d), complementing Figs. 9(b), (d) for the scattering. One can see that there always exists a region where ${\sigma}_{\text{ext}}=0$ when${\epsilon}_{\text{L}}>{\epsilon}_{\text{thr}}$, as indicated by the arrow in Fig. 13(d). The situation at negative (red) detuning is similar above the threshold.

Figure 14 compares the saturated and unsaturated absorption (a), (c), (e), and extinction (b), (d), (f), cross-sections in $(\omega ,\text{\hspace{0.17em}}{\epsilon}_{\text{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 ${\epsilon}_{\text{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) ), 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).

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, whichcan 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 ${\sigma}_{\text{ext}}=0$ can be achieved on both sides of the resonance. The condition ${\sigma}_{\text{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 ${\sigma}_{\text{ext}}=0$, which is equivalent to a condition ${{d}^{\u2033}}_{3}=0$ 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 ${\epsilon}_{\text{L}}$, external field ${e}_{3}$ and frequency $\omega $, 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 self-consistent 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}_{\text{s}}[c{m}^{-3}]$, and the (CGS) “atomic” polarizability

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:

*ε*of the spasers + MM composite is found from here, loss compensation corresponds to the conditionThis depends on the properties of the background MM being compensated, ${\epsilon}_{\text{b}}$, and on the concentration of spasers ${N}_{\text{s}}$. 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.

## Acknowledgments

We thank G. Kristanz (JKU Linz) for technical assistance. Funding from the European Research Council (ERC Starting Grant 257158 ‘Active NP’), and partial support from the ONR Grant (MURI N00014-13-0649) are gratefully acknowledged.

## 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). [CrossRef] [PubMed]

**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). [CrossRef] [PubMed]

**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). [CrossRef]

**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). [CrossRef] [PubMed]

**5. **J. B. Khurgin and G. Sun, “Comparative analysis of spasers, vertical-cavity surface-emitting lasers and surface-plasmon-emitting diodes,” Nat. Photonics **8**(6), 468–473 (2014). [CrossRef]

**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). [CrossRef] [PubMed]

**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). [CrossRef] [PubMed]

**8. **N. M. Lawandy, “Localized surface plasmon singularities in amplifying media,” Appl. Phys. Lett. **85**(21), 5040–5042 (2004). [CrossRef]

**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). [CrossRef]

**10. **A. Y. Smuk and N. M. Lawandy, “Spheroidal particle plasmons in amplifying media,” Appl. Phys. B **84**(1-2), 125–129 (2006). [CrossRef]

**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). [CrossRef] [PubMed]

**12. **I. E. Protsenko, “Quantum theory of dipole nanolasers,” J. Russ. Laser Res. **33**(6), 559–577 (2012). [CrossRef]

**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). [CrossRef]

**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). [CrossRef] [PubMed]

**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). [CrossRef] [PubMed]

**16. **M. I. Stockman, “Spaser action, loss compensation, and stability in plasmonic systems with gain,” Phys. Rev. Lett. **106**(15), 156802 (2011). [CrossRef] [PubMed]

**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). [CrossRef] [PubMed]

**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). [CrossRef] [PubMed]

**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). [CrossRef]

**22. **M. I. Stockman, “The spaser as a nanoscale quantum generator and ultrafast amplifier,” J. Opt. **12**(2), 024004 (2010). [CrossRef]

**23. **V. M. Parfenyev and S. S. Vergeles, “Quantum theory of a spaser-based nanolaser,” Opt. Express **22**(11), 13671–13679 (2014). [CrossRef] [PubMed]

**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). [CrossRef]

**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). [CrossRef] [PubMed]

**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). [CrossRef] [PubMed]

**27. **A. Pusch, S. Wuestner, J. M. Hamm, K. L. Tsakmakidis, and O. Hess, “Coherent amplification and noise in gain-enhanced nanoplasmonic metamaterials: a Maxwell-Bloch Langevin approach,” ACS Nano **6**(3), 2420–2431 (2012). [CrossRef] [PubMed]

**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). [CrossRef] [PubMed]

**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). [CrossRef] [PubMed]

**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). [CrossRef] [PubMed]

**36. **H. P. Zhang, J. Zhou, W. B. Zou, and M. He, “Surface plasmon amplification characteristics of an active three-layer 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). [CrossRef]

**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). [CrossRef]

**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). [CrossRef]

**40. **X. Fan, Z. Shen, and B. Luk’yanchuk, “Huge light scattering from active anisotropic spherical particles,” Opt. Express **18**(24), 24868–24880 (2010). [CrossRef] [PubMed]

**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). [CrossRef]

**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).

**45. **P. B. Johnson and R. W. Christy, “Optical constants of noble metals,” Phys. Rev. B **6**(12), 4370–4379 (1972). [CrossRef]

**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). [CrossRef] [PubMed]

**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). [CrossRef]

**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). [CrossRef]

**50. **V. M. Parfenyev and S. S. Vergeles, “Intensity-dependent frequency shift in surface plasmon amplification by stimulated emission of radiation,” Phys. Rev. A **86**(4), 043824 (2012). [CrossRef]

**51. **O. Svelto and D. C. Hanna, *Principles of Lasers* (Springer, 2009).

**52. **R. W. Boyd, *Nonlinear Optics* (Academic Press, 2008).

**53. **I. A. Fedorov, V. M. Parfenyev, S. S. Vergeles, G. T. Tartakovsky, and A. K. Sarychev, “Allowable Number of Plasmons in Nanoparticle,” JETP Lett. **100**(8), 530–534 (2014). [CrossRef]

**54. **J. B. Khurgin, “Ultimate limit of field confinement by surface plasmon polaritons,” Faraday Discuss. **178**, 109–122 (2015). [CrossRef] [PubMed]

**55. **J. B. Khurgin and G. Sun, “How small can “Nano” be in a “Nanolaser”?” Nanophotonics-Berlin **1**, 3–8 (2012).

**56. **L. D. Landau, E. M. Lifshit︠s︡, and L. P. Pitaevskiĭ, *Electrodynamics of Continuous Media* (Pergamon, 1984).

**57. **T. Poston and I. Stewart, *Catastrophe Theory and its Applications* (Dover Publications, 1996).