Power Balance and Temperature in Optically Pumped Spasers and Nanolasers

Spasers and nanolasers produce a significant amount of heat, which impedes their realizability. We numerically investigate the farfield emission and thermal load in optically pumped spasers with a coupled electromagnetic/thermal model, including additional temperature discontinuities due to interfacial Kapitza resistance. This approach allows to explore multiple combinations of constitutive materials suitable for robust manufacturable spasers. Three main channels of heat generation are quantified: metal absorption at pumping and spasing wavelengths and nonradiative relaxations in the gain material. Out-radiated power becomes comparable with absorption for spasers of realistic dimensions. Two optimized spaser configurations emitting light near 520 nm are compared in detail: a prolate metal-core/gain-shell and an oblate gain-core/metal-shell. The metal-shell design, which with the increasing size transforms into a metal-clad nanolaser, achieves an internal light-extraction efficiency of 22.4%, and stably operates up to several hundred picoseconds, an order of magnitude longer than the metal-core spaser.


Purcell-enhanced spontaneous emission
The total spontaneous emission rate of the chromophores from levels 2 to 1 (γ 21 ) includes several channels. The electronic excitation in the dye system can either decay into a resonant (in the present case dipolar) surface plasmon, or it can decay into higher order (non-resonant, multipolar) surface plasmons. It can also decay radiatively (into far-field photons) or nonradiatively (into phonons). Spontaneous emission contributes to noise in a spaser, S1,S2 which we are not able to simulate within the current frequency-domain framework. Thus, we assume that our systems are "good" single-mode spasers, where the spontaneous emission primarily yields dipolar surface plasmons.
The spontaneous emission rate is Purcell-enhanced, S3 since the nanoparticle acts as a resonator with a very small mode volume. Quantum-mechanics predicts that the stimulated emission into a given single mode is always a factor n s stronger than the spontaneous emission into it (see Section 8.3 in Yariv, S4 Supplement 2 in Khurgin et al. S2 and Supplement Section 2 in Kewes et al. S5 ), where n s is the number of resonant surface plasmons: n s can be expressed in terms of the total electromagnetic (spasing) energy in the system (see §80 in Landau and Lifshitz S6 ), To get an estimate for γ 21 used in the rate equations results in the main text, it is convenient to average the number of transitions per unit volume over the gain material,

S2
where N 2 is the population density and n 2 the total population of level 2. Inserting the equation for the spasing rate W s , similar to Eq. (3) in the main text, yields (assuming spasing in the center of the line with ω s = ω 21 ) The ratio in parentheses quantifies the overlap between the density profiles of the electromagnetic energy w and the population inversion N 2 . The result has the dimension of inverse volume: The dimensionless factor f = W G /W is the fraction of electromagnetic energy stored in the gain material (neglecting small magnetic contributions) vs. total electromagnetic energy in the system, while V G is the modal volume in the gain material, modified by the distribution of population inversion. We assume that half of the energy is located in the gain material, 2f = 1, approximate V G by the gain volume, (which, for our geometries, is related to the volume of the nanoparticle), and use the averaged value γ 21 =γ 21 everywhere. This leads to our final approximation for the spontaneous emission rate: The last equality in Eq. (S7) was obtained expressing σ 21 via the bulk (unmodified) radiative decay rate in the host medium, γ 21,bulk (refractive index n h = √ ε h ), (see Eq. (20) in Arnold et al. S7 and references therein). It contains the quality factor of the atomic line Q L = ω/γ L , S3 rather than that of the resonator Q R = ω/γ R , which enters the conventional Purcell factor  in Saleh and Teich S8 ). The latter is derived for a narrow atomic line γ L γ R , while for organic dyes (or semiconductors) the situation is reversed, or the linewidths are comparable. In such situations, the emission linewidth becomes the dominant factor, see Eq. (37) in Khurgin,S9 or Eq. (40) in the Supplementary Materials to Khurgin and Sun,S2 which explicitly discusses this issue.
The first expression in Eq. (S7), (Eq. (5) in the Main Text) elucidates the corpuscular origin of the Purcell enhancement as an increase in collision frequency between the emitter/absorber and the photon, which "moves" with the velocity c/ √ ε h within the small modal volume. To incorporate all this into the model, going beyond the mesoscopic description, would be an overshot in accuracy, diluting the main message of the work, which lies in the analysis of retardation, light extraction and thermal effects, keeping the key parameters within the physically admissible range.
Similarly, quantum coherence effects are smeared out by both strong dephasing of the chromophores, and the fast decay of plasmons. This can be illustrated as follows. The quantum dynamics of a spaser can be described by three equations: see (4)-(6) in ref., S1 or (64)- (67) in ref., S10 or (34) in ref. S11 (the latter discusses also more elaborate models): Here, ρ and n are complex amplitudes of non-diagonal and diagonal elements of the density matrix, which describe polarization and population inversion; a is the complex amplitude of S4 (quasi-classical) plasmon number operator; ω, ω 21 and ω M are the frequencies of operation, atomic transition and plasmon, respectively. Γ 21 = γ L /2 is the decay of non-diagonal part of density matrix, (i.e., dephasing rate of chromophores, decay of macroscopic polarization), γ 21 is the depopulation rate of the upper level, including non-radiative and (Purcell-enhanced) radiative contributions, and γ M is the plasmon decay rate. g is the pumping rate parameter, Ω is the (normalized) Rabi frequency of the spasing transition, and n tot the total number of (identical) chromophores.
These equations contain quantum coherence effects. For example, if the spasing field a is enforced externally (and pumping g is correspondingly removed), all frequencies are equal and damping and dephasing are absent, the last terms in the first two equations in (S8) are oscillatory, with the Rabi frequency 2|aΩ|.
All three rates in equations (S8): the dephasing Γ 21 , Purcell-enhanced relaxation γ 21 and plasmon decay γ M , damp any oscillatory processes. These rates are in the range 3 × 10 11 to 2 × 10 14 s −1 . Thus, any quantum coherence oscillations will be suppressed at least within few picoseconds, which is well below the timescales of interest for practical spaser operation.
For such long times, calculations based on quantum density matrix (optical Bloch) equations for the spaser (section 2.1 in ref. S1 ), which do include (strongly damped) coherent effects, result in the (quasi-static) threshold expression (Eq. (82) in similar ref. S10 ), which is identical to the one derived from the purely electrodynamic considerations (using appropriate microscopic expression for ε gain ), as noted in ref. S12 (end of section 3 there). Post-threshold behavior in long pulses is also the same, as discussed in ref. S7

Heat source in the gain material
The gain material is heated by non-radiative relaxation and decay of electronic states into various vibronic excitations of the dye chromophores. The corresponding heat source follows S5 from the corresponding terms in an idealized 4-level system: Q G =hω 20 γ 20,nr N 2 + (hω 30 γ 30,nr +hω 32 γ 32 )N 3 +hω 10 γ 10 N 1 .
Here, we additionally included non-radiative relaxations from levels 3 → 0 and 2 → 0 (which are not shown in Fig. 1). The term withhω 20 is more appropriate for non-radiative transitions thanhω 21 . Indeed, in a typical 4-level dye, the vibrational-rotational structure is very complex and includes many eigenmodes. As a result, the electronic excitation energy will primarily thermalize through available vibronic states directly to the lowest energy level, which is 0.

S6
The stationary solution for a continuously heated, spherical nanoparticle with Kapitza resistance In the ambient, the stationary heat equation is with the solution where T 2 is the stationary surface temperature rise (with respect to the background temperature T 0 ) on the ambient side and a is the nanoparticle radius. The flux J over the boundary S is continuous and fulfills Eq. (11), where k a is the thermal conductivity of the ambient and T 1 is the stationary surface temperature rise on the particle side of the interface. Then, the relative change in temperature at the boundary is which is inversely proportional to the nanoparticle radius a.

shell structures
Before solving Maxwell's equations for a spaser numerically, it is often a good idea to look at their quasistatic approximation (i. e., neglecting retardation). It can be used when the S7 structure under study is much smaller than the typical wavelength. Practically, this means that the wave vector k is neglected in Maxwell's equations. The resulting equations can often be solved analytically and do not contain magnetic fields.
Bohren and Huffman S14 derived a quasistatic solution for the polarizability α of a confocal spheroidal core-shell structure in an infinite ambient (see Section 5.4, Eq. (5.35)), where ε 1 , ε 2 and ε 3 are the dielectric functions of core, shell and ambient, V is the total volume of the structure and h < 1 is the volume fraction of the core spheroid to the entire structure. The function L i = L(e i ) depends on (i) the eccentricities e 1 and e 2 of the core and shell spheroids, (ii) the polarization of the incident electric field and (iii) on whether the spheroids are oblate or prolate. Note that the definition of the eccentricity depends on the shape of the spheroid, e = 1 − c 2 a 2 with c < a (oblate), where c is the semi-axis parallel to the axis of revolution. If the incident electric field is parallel to the axis of revolution, then , If the polarization of the incident field is perpendicular to the axis of revolution, then where L z is the appropriate function from Eqs. (S19). By neglecting saturation effects (which is usually justified near the threshold) and tuning the spaser to the 2 → 1 transition of the chromophores (meaning ω s = ω 21 ), the gain dielectric function (1) simplifies to at the threshold. The spasing threshold corresponds to the polarization going to infinity in the unsaturated case, meaning that the fields in the structure become very high for very small incident fields. In other words, sustained oscillations become mathematically possible as a solution of a homogeneous equation, i. e., without incident field at all. This yields the following condition for the denominator of Eq. (S17), which can be solved to estimate the spasing threshold. Since the complex-valued condition (S22) should be fulfilled at the threshold wavelength ω thr , we can numerically look for solution-pairs (ω thr , ε L,thr ). By doing a sweep over several aspect ratios, we can estimate which aspect ratio is favorable in terms of spasing frequency and threshold. For gain-core/metal-shell structures, there exist three solution branches per polarization of the incident field that all fulfill condition (S22). They represent different surface plasmon oscillations: S15 the highest-wavelength solution branch corresponds to the surface charges at the inner and outer surface of the metal shell oscillating in phase ("bonding" plasmon), while for the lowest-wavelength branch the charges oscillate in anti-phase ("antibonding" plasmon). It can be shown that the middle-wavelength branch always has the highest gain threshold and is thus of no interest for us. When plotted over the aspect ratio, the three solution branches only exist up until some aspect ratio κ 0 . At κ = κ 0 , two branches merge and vanish simultaneously: for κ > κ 0 , only one branch remains (see e. g., Fig. S1b). Thus, finding a single solution point is usually not sufficient for a complete analysis of a spaser. Table S1: Quasistatic results for eight different core/shell spaser configurations with shell thickness (along all semi-axes) equal to major semi-axis h = a. The core aspect ratio is varied from 1 (spherical) to 8 (spheroidal), larger aspect ratios would lead to unrealistically thin structures. The shape of the spasers can be either prolate (cigar-shaped) or oblate (pancakeshaped) spheroidal. The pumping field can either be parallel (E z) or perpendicular (E x) to the axis of revolution z. In the reference FEM calculations (see Main Text), the values for the shell thickness and major semi-axis are h = a = 23 nm for the prolate metal-core/gain-shell and h = a = 30 nm for the oblate gain-core/metal-shell spaser. Every gain-core/metal-shell configuration has three solution branches for the quasistatic spasing condition (S22), each corresponding to different surface plasmon oscillations. S15 The middlewavelength branch was omitted in the table, since it always has much higher gain threshold.
For a system with a small gain core and a thick metal shell, however, quasistatic predictions can be very accurate. Table S1 shows the quasistatic results of eight core/shell spaser configurations with equal major semi-axis and shell thickness a = h. The aspect ratio is varied from 1-8 to avoid unrealistically thin structures. As long as geometrical similarity a = h is preserved, these results hold for quasistatic structures of arbitrary size. Here, we imply a = h = 30 nm for an oblate spaser, and a = h = 23 nm for a prolate one. Strictly speaking, Eq. (S22) only holds for confocal (quasistatic) spheroids -however, it still gives a reasonable approximation for Figure S1: Quasistatic calculations for three spaser configurations: prolate metal-core/gainshell E z (red), oblate metal-core/gain-shell E x (green) and oblate gain-core/metalshell E z (blue). (a) Gain thresholds and (b) generation wavelengths over the range of realistic aspect ratios. The dotted grey lines indicate the aspect ratios with a generation wavelength of 520 nm, where the gain thresholds have local minima. For the gain-core/metalshell structure, the gain thresholds of the first and second branches are out of range in the plot (see Table S1).
our geometries.
Three spaser configurations have low gain thresholds (ε L,thr < 0.15) and generation frequencies in the visible range (see Figs. S1a-b). We choose to simulate the oblate gain-core/metalshell structure with E z and the prolate metal-core/gain-shell structure with E z. The spasers are tuned to a generation wavelength of 520 nm: there, the gain threshold has a local minimum and the (quasistatic) aspect ratios are realistic.
Due to appreciable retardation effects, the quasistatic approximation provides a rough es- Table S2: Exact numeric threshold values and core aspect ratios for the spaser configurations under study.
Spaser configuration κ ε L,thr λ thr Oblate gain-core/metal-shell 6.25 0.1133 519.96 nm Prolate metal-core/gain-shell 1.94 0.1159 519.78 nm timation of the gain threshold, generation frequency and necessary shape of a core-shell spaser. Using FEM simulations, we can find the exact numeric values, see Table S2. The final geometries are shown in Fig. 2 in the Main Text.

Temperature-dependent Drude model
Since the thermal and electromagnetic problems are coupled via a temperature-dependent dielectric function of the metal component in the spaser, a proper description of the thermal behavior of the latter is necessary. This has been extensively measured for various metals.
For example, an overview of the temperature dependence of the dielectric function of gold can be found in Perner S16 (Section 2.1). For silver, measurements by Sundari et al. S17 exist, but since their imaginary part is questionably high and their real part changes with temperature, we choose not to use their data. A different approach would be to correct a reliably measured dielectric function using a temperature-dependent Drude model, .
Here, ε JC Ag are the well-known measurements by Johnson and Christy S18 (measured at room temperature T 0 ) and ε Dr Ag is a Drude interpolation of the data with temperature-dependent plasma and collision frequency. Do not confuse the Drude plasma frequency ω P with the S12 pumping frequency of the spaser ω p .
The temperature-dependence of the plasma frequency is determined by the density and effective mass of free electrons. S19,S20 Assuming an isotropic material, the density ρ changes due to thermal expansion with where α is the linear thermal expansion coefficient. Since α is on the order of 10 −5 K −1 for metals, the total change of ω P over the available temperature range is negligible. According to Reddy et al.,S20 the change in ω P due to a temperature-dependent electronic effective mass is less than 10 % from 300 K to 900 K. Since both of these effects are small, we omit the temperature-dependence of the plasma frequency for all further purposes.
The thermal behavior of the Drude collision frequency depends on the electron-phonon interaction. According to Ujihara, S21 the temperature-dependence of electron-phonon scattering can be described via where γ(T 0 ) is the collision frequency at room temperature calculated via Drude interpolation and θ is the Debye temperature. Using the temperature-dependent measurements for the Debye temperature by Simerská S22 (θ(300 K) ≈ 210 K and θ(1000 K) ≈ 190 K), the temperature-dependent dielectric function for silver can be calculated. The parameters of the Drude fit arehω P = 9.169 eV,hγ = 0.021 eV and ε ∞ = 3.58. The data was fitted with a least-square method in the range 190-1900 nm, which emphasizes longer wavelengths. The results agree with the measurements of the temperature-dependent dielectric function of silver thin films by Reddy et al. S20

S13
Because of ω γ, Eq. (S23) can be Taylor-expanded as: This shows that the changes in the real part of the dielectric function are negligible. The integrand of f in Eq. (S27) can also be Taylor-expanded at elevated temperatures, which yields for γ: Thus, the imaginary part of ε Ag increases approximately linearly with temperature.
Quasistatic estimation for the temperature dependence of the spasing field at the threshold generation frequency In a quasistatic core/shell spaser, the electric field E in the gain core is almost constant and is given by Eq. (34) in Arnold et al.: S7 Considering it on resonance (ω = ω thr ) and using the two-material approximation for the spasing threshold (Eq. (11) in Arnold et al. S12 ), yields

S14
The spasing field in the metal, E s , is linearly related via boundary conditions to the (approximately constant) electric field in the core of the structure. Thus, the expected temperature dependence of the metal spasing field E s at the threshold generation frequency ω thr is: This is Eq. (14) from the main text. Due to approximations used, especially the assumptions of the gain core and the two-materials threshold (S31), it can serve as a guideline only. With three materials, one can find the threshold from the zero of the denominator in (S17) (expression (S22)). This replaces h by h +f , where f is a cumbersome combination of parameters, which, however, does not change the qualitative trends. Numerical results confirm that a similar decrease of metal field with temperature persists also for the prolate metal-core/gain-shell geometry. However, such analytical estimations are less justified there due to variations in saturation across the gain material.
Cooling time Figure S2 shows the same data as in Fig. 3, but in a linear scale for both axes and on a far longer time scale including a purely thermal cooling simulation. The metal-core/gain-shell structure has much shorter cooling times. This is related to the higher temperature gradients leading to a larger heat flux to the ambient; besides, in non-stationary regime, the cooling time for an arbitrary structure is typically similar to its heating time.

Polystyrene as gain host material
In Kristanz, S13 we also simulated similar spasers with polystyrene (PS) as gain host material. PS was recently used as a matrix for organic dyes in microresonators, S23-S25 and Figure S2: The average temperature T in the respective spaser components (metal and gain). Solid curves are for the oblate gain-core/metal-shell, dashed curves for the prolate metal-core/gain-shell spaser. The maximum operation time is 650 ps and 60 ps, respectively. The dotted grey lines indicate the end of spaser operation and the start of a purely thermal cooling simulation.
was also modeled in other spaser simulations. S7,S12 The PS-based spasers turn out to be much more heat sensitive than the silica-based ones because of lower thermal stability of PS (melting point is 510-540 K instead of 1986 K for silica S26,S27 ) and its inferior thermophysical properties (thermal conductivity is 0.155 vs. 2.4 W/(m K), thermal diffusivity is 1.2 × 10 −3 vs. 8.5 × 10 −3 cm 2 /s for silica S27-S32 ). As a result, the PS spasers can only operate for 110 ps (oblate gain-core/metal-shell) and 45 ps (prolate metal-core/gain-shell), before the PS melting point is reached. This compares unfavorably with the 650 ps and 60 ps for the silica case.