Spatio-temporal Modeling of Lasing Action in Core-Shell Metallic Nanoparticles.

Nanoscale laser sources based on single metallic nanoparticles (spasers) have attracted significant interest for their fundamental implications and technological potential. Here we theoretically investigate the spatio-temporal dynamics of lasing action in core-shell metallic nanoparticles that include optically pumped four-level gain media. By using detailed semiclassical simulations based on a time-domain generalization of the finite-element method, we study the evolution of the lasing dynamics when going from a spherical case to an elongated nanorod configuration. Our calculations show that there exists an optimal nanoparticle elongation that exhibits significantly improved lasing threshold and slope efficiency over those obtained for its spherical counterpart. These results are accounted for in terms of a coupled-mode theory analysis of the variation with elongation of the light confinement properties of localized surface plasmons. This work could be of importance for further development of nanoscale light sources based on localized surface plasmon resonances.

The efficient generation of coherent light at the nanoscale has become one of the ultimate frontiers of laser research. In addition to their fundamental interest, nanometric laser sources have the potential for enabling novel applications in a variety of areas, from spectroscopy, sensing, and bioimaging to emerging quantum technologies. The advent of plasmonics[1] has opened up whole new prospects for realizing such light sources. Plasmonic resonances (electromagnetic excitations arising from the interaction of light with conduction electrons at metallic interfaces) allow confining and enhancing light fields at deep-subwavelength scales.[2] The associated increase of light–matter interaction strength can be exploited to boost the effective gain provided by an active medium adjacent to a suitably designed plasmonic structure. This physical mechanism not only has been shown to allow compensation of the significant absorption losses of metals in the visible regime[3−9] but can also lead to self-sustained laser oscillations at the nanoscale.[10−32]

Of special interest are plasmonic lasers based on localized surface plasmons (LSPs)[33−36] supported by single metallic nanoparticles. These structures (often referred to as spasers(10)) consist essentially of a core–shell metallic nanoparticle, designed so that its lowest frequency LSP is tuned to the emission wavelength of an externally pumped active medium incorporated into the shell. The versatility of this class of structures as a testbed to explore fundamental laser phenomena at the nanoscale, together with the availability of efficient approaches to tailor the nanoparticle response,[37−46] has stimulated extensive research work.[10,13,15,24−26,31] In this context, the full spatio-temporal modeling of spasers could contribute to further development and optimization of nanoscale light sources based on LSPs.

In this work we report on a theoretical study of lasing action in metallic nanoparticles coated by an optically pumped four-level gain media, focusing on the evolution of the lasing characteristics when going from the conventional spherical configuration to an elongated nanorod geometry. We investigate this problem by means of detailed three-dimensional (3D) semiclassical simulations based on a time-domain generalization of the finite-element method. These simulations are able to track the nonlinear spatio-temporal dynamics arising from the interaction of the gain medium with the highly nonuniform field distribution characterizing the system. Our calculations show how by departing from the conventional spherical configuration through the elongation of the nanoparticle, it is possible to significantly reduce the lasing threshold, while at the same time increase the laser slope efficiency. Moreover, we account for these findings by using a semianalytical coupled-mode theory analysis.

Results and Discussion

Figure a renders a schematic view of the core–shell nanoparticles under study. The core consists of a metallic nanorod, which can be seen as a cylindrical central body of length Lrod and radius Rin, capped by two hemispherical ends of radius equal to that of the cylindrical middle part. This core is coated by a dielectric shell (of index nh) featuring the same external shape as the core; that is, the length Lrod is maintained but the hemispherical caps now have a radius Rout > Rin. This shell is assumed to contain four-level organic dye molecules, which when optically pumped can provide optical gain to the system. In addition, we consider that the whole structure is embedded in vacuum.

Figure 1

(a) Schematic view of the considered core–shell metallic nanoparticle. It consists of a gold nanorod coated by a dielectric shell (of index nh = 1.46). The shell is doped with organic dye molecules, which upon optical pumping (indicated by a red arrow) can provide enough optical gain to yield laser emission (blue arrows). The geometrical parameters defining the system, as well the reference system used in this work, are also included. (b) Main panel: Extinction cross section spectra calculated for nanoparticles with different values of the nanorod length Lrod, going from Lrod = 0 (spherical case; see red line) to Lrod = 100 nm (elongated nanorod configuration; black line). For Lrod = 0, we chose Rin = 7 nm and Rout = 22 nm. In the rest of configurations, we kept Rin = 7 nm but reduced Rout as Lrod grows, so that the volume of the dielectric shell is maintained for all cases. The system is illuminated by a plane wave linearly polarized along the long axis of the nanoparticle (see sketch in the bottom corner). Inset: Same as main panel but now the incident illumination is polarized perpendicularly to the long axis of the nanoparticle (see sketch included in the inset).

We start by analyzing the optical response of the passive system (i.e., the linear response of the system assuming a negligible concentration of dye molecules in the shell). Main panel of Figure b displays the calculated extinction cross section (ECS) spectra for nanoparticles with different values of elongation, ranging from Lrod = 0 (i.e., the spherical case; see red line) to Lrod = 100 nm (elongated nanorod configuration; black line). In all cases, we assumed a gold core and an nh = 1.46 shell. For the spherical case, we chose Rin = 7 nm and Rout = 22 nm, so that the case Lrod = 0 corresponds to the experimentally realized spaser considered in ref (13). For nanoparticles featuring larger values of Lrod, we fixed Rin = 7 nm and reduced continuously the value of Rout as Lrod is increased, so that the volume of the dielectric shell is the same for all considered configurations. As we show below, this is an important ingredient when analyzing the lasing characteristics of the active structures. In all cases, the system is illuminated by a plane wave, linearly polarized along the long axis of the nanoparticle and incident normally to that axis. For comparison, the inset of Figure b displays the ECS spectra for the same configurations shown in the main panel, but now the incident plane wave is polarized perpendicularly to the long axis of the nanoparticle. These linear calculations have been performed with the conventional frequency-domain implementation of the finite-element method (FEM),[47] in which a Drude–Lorentz form—fitted to available experimental data[48,49]—is used to model the metallic region of the system. In addition, we use absorbing boundary conditions placed in the far field to simulate the response of a single isolated nanoparticle.

As observed in main panel of Figure b, for the long-axis incident polarization, the ECS spectra are dominated by a broad resonant peak, whose location and maximum value increase linearly with Lrod. This peak corresponds to the dipolar LSP resonance (usually referred to as longitudinal LSP), which has been extensively studied for its subwavelength-enhanced field profile.[33−36] The near-field distributions calculated at resonance for three representative values of Lrod are shown in Figure a (left, center, and right panels display the cases Lrod = 0, Lrod = 20 nm, and Lrod = 60 nm, respectively). As seen, as Lrod grows, the local field enhancement around the caps of the nanorods (lobes of the dipolar profile in the spherical case) increases significantly. The local field enhancement obtained by departing from the spherical geometry is accompanied by an increase of the temporal confinement of the electromagnetic (EM) fields in the system. From separate simulations, based on computing the decay rate of the total EM energy of the system,[50,51] we obtained that the Q-factor of the system grows from Q ≈ 8 for the spherical case to Q ≈ 20 for the larger values of Lrod. As expected for these nanoparticle sizes, these Q-factors are dominated by absorption losses; we obtained radiative Q-factors of about 2 orders of magnitude larger than absorption ones. A more detailed account of the evolution of Q with Lrod will be given below in the context of the coupled-mode theory analysis.

Figure 2

(a) Electric-field profiles (|E|, normalized to the incident field) for the structure displayed in Figure a, as computed for different values of the nanoparticle elongation, Lrod. Left, center, and right panels correspond to Lrod = 0, Lrod = 20 nm, and Lrod = 60 nm, respectively. The rest of the parameters defining the structure are the same as those considered in Figure b. All cases have been calculated at the corresponding resonant wavelengths and assuming that the incident illumination is linearly polarized along the long axis of the nanoparticle. (b) Same as (a), but now the incident illumination is polarizared along the short axis of the nanoparticle. A sketch of the incident polarization and the reference system are included in the leftmost part of both panels. The incident wave propagates along the y-direction (perpendicular to the page) for both polarizations.

In contrast, for the short-axis incident polarization we obtain much smaller peak values of ECS (see inset of Figure b). These resonances correspond to the so-called dipolar transversal LSP resonances.[1] Note that dipolar longitudinal and transversal LSP resonances can be seen as arising from degeneracy breaking (via symmetry breaking introduced by elongation) of their counterpart modes in the spherical case. As observed in the inset of Figure b, the spectral position of the transversal LSP resonant peak remains almost insensitive to Lrod. In addition, our calculations show that the Q-factors of these resonances are also essentially independent of Lrod (specifically, we obtained that Q decreases monotonically from Q ≈ 8 to Q ≈ 5 when Lrod is increased from Lrod = 0 to Lrod = 100 nm). The described characteristics of longitudinal LSPs, together with the possibility of tuning their spectral response just by varying Lrod, make longitudinal resonances ideal candidates for creating nanoparticle-based laser systems. Therefore, in the rest of the work, we focus on this class of plasmonic resonances.

We turn now to investigate how lasing action emerges in the considered systems when a significant concentration of optically pumped dye molecules is incorporated in the shell. Before proceeding with a detailed discussion of the corresponding simulation results, we briefly summarize the numerical framework used to tackle this problem (a more detailed account can be found in refs (52−54) and in the Methods section). A semiclassical simulation of the class of plasmonic lasers considered here poses two main challenges from a numerical standpoint. First, the highly nonlinear dynamics of the problem involves two very different time scales, namely, the one associated with the optical-frequency plasmon oscillations (on the order of fs) and the one corresponding to the evolution of electronic populations of the dye molecules forming the gain medium (on the order of ns). To overcome this difficulty, the employed theoretical framework exploits the fact that two of the characteristic frequencies of the system are known, i.e., the absorption and emission electronic transitions of dye molecules (denoted by ωa and ωe, respectively). This allows writing the following ansatz for the functional form of the vector potential of the system: A(r, t) = A(r, t) exp(−iωat) + A(r, t) exp(−iωet), where A(r, t) (i = a, e) are slowly varying complex amplitudes (once A(r, t) is known, the E-field profile can be obtained just by using E = −∂A/∂t; note that the temporal gauge condition is assumed in this approach). This ansatz does not introduce any constraint on the spatial profile of the nanoparticle modes or their harmonic oscillation frequencies (frequency-pulling effects are reflected in additional harmonic time dependences of A(r, t), i.e, A(r, t) = Ã(r, t) exp(iΔωt), with Δω being the corresponding frequency shift with respect to ω). We also assume a similar expansion for P(r, t) inside the gain medium; in this case the slowly varying amplitudes P(r, t) represent the stimulated absorption (Pa(r, t)) and emission (Pe(r, t)) of photons (see Methods).

The second important challenge in these simulations is the highly nonuniform field distribution that characterizes the studied structures at the subwavelength scale. In this regard, the adaptive meshing of the FEM method is ideally suited to solve the problem. To benefit from this type of meshing (usually only employed in frequency-domain formulations), we used a FEM algorithm[47] to solve, in the time domain, the so-called weak form of the field equation governing the spatio-temporal dynamics of A.[55] We emphasize that the presence in the system of optically pumped dye molecules (which we model as four-level emitters) makes this field equation nonlinear. This nonlinearity arises from the coupling of P to the corresponding electronic population inversions of the dye molecules, these electronic populations being, in turn, driven by the product E(∂P/∂t)[56] (see Methods). Although in this work we focus on a particular type of four-level organic dye molecules, the described approach is general and can be used to simulate other classes of optically pumped gain media (see for instance ref (53), where a host–guest molecular gain system assisted by Förster energy transfer is considered). Therefore, we expect the conclusions of this work to be general and applicable to spasers based on different types of gain media. Finally, we note that the applied formalism is a purely semiclassical approach that cannot account for quantum fluctuations associated with spontaneous emission or thermal effects. Therefore, the numerical framework used in this work cannot provide accurate predictions of the evolution of the laser linewidth with the pumping strength. A fully quantum generalization of our approach to incorporate quantum fluctuations is out of the scope of the present work and will be addressed elsewhere.

Figure a and b summarize the lasing dynamics and lasing characteristics of the active counterparts of two of the structures considered in Figures and 2 (Figure a and b correspond to Lrod = 0 and Lrod = 60 nm, respectively). Specifically, the main panels of Figure a and b show the time evolution of the ωe component (the lasing component) of the E-field amplitude of the system, spatially averaged over the computational domain of our simulations (we denote this quantity as ⟨|E(r, t)|⟩). For comparison, the results for several values of the pump intensity |Epump|2 are also included in these panels (red, green, and blue correspond respectively to |Epump|2 = 9|Eth|2, |Epump|2 = 16|Eth|2, and |Epump|2 = 25|Eth|2, with |Eth|2 being the threshold pump intensity needed to reach the lasing regime for the corresponding value of Lrod). The normalization constant E0 in Figure a and b stands for the saturation electric field of the gain medium,[57], where τ21 and σa are the stimulated emission lifetime and the absorption cross section of the lasing transition, respectively (see Methods); ℏ is the reduced Planck constant, ϵ0 is the vacuum permittivity, and c is the speed of light in a vacuum. In all cases, we assume that the structure is illuminated by a plane wave of wavelength λp = 490 nm, linearly polarized along the long axis of the nanoparticle. In addition, all laser simulations shown in this work have been obtained by tuning ωe to the corresponding LSP dipolar resonance obtained for each Lrod value. The rest of the parameters modeling the emitters are the following: τ21 = 1.6 ns, τ32 = τ10 = 100 fs, Γa = Γe = 1/(10 fs), and σa = σe = 2.55 × 10–16 cm2 (see Methods for the definition of each parameter). These values correspond to the emitters used in ref (13). Finally, the total density of molecules in the gain medium is N0 = 3.7 × 1026 m–3. We recall that the geometry of the considered configurations is such that the volume of the active shell surrounding the metallic core is maintained for all considered values of Lrod.

Figure 3

(a) Time evolution of the spatially averaged laser field amplitude ⟨|E(r, t)|⟩ (normalized to the saturation electric field, E0), as computed for a spherical spaser based on the active counterpart of the core–shell nanoparticle considered in Figure a. The results for three different values of the amplitude of the pump (|Epump|) are displayed (red, green, and blue lines correspond respectively to |Epump| = 3|Eth|, |Epump| = 4|Eth|, and |Epump| = 5|Eth|, where |Eth| is the corresponding pump field amplitude at the lasing threshold). In all cases the pump wavelength is λp = 490 nm, whereas the lasing wavelength is λe = 530 nm. (b) Same as panel (a) but now for the active counterpart of the structure displayed in Figure c (an elongated core–shell nanoparticle with Lrod = 60 nm). In this case, λe = 890 nm (λp is the same as in panel (a)). Both panels include schematic views of the corresponding core–shell nanoparticles in the top-left part. Insets of both panels display the steady-state values for the averaged field intensity ⟨|E(r, t)|2⟩ as a function of the pump intensity |Epump|2 (the colored circles represent the steady-state values of the corresponding cases with the same color shown in the main panels). For Lrod = 0, |Eth| = 1.5E0 is obtained, while |Eth| = 0.3E0 for Lrod = 60 nm.

As observed in Figure a and b (main panels), for both Lrod = 0 and Lrod = 60 nm we obtain the canonical features of lasing dynamics, consisting of a series of sudden spikes in the laser field signal that settles down to a steady-state value for long times.[56] As seen, the case Lrod = 60 nm displays larger steady-state values than those corresponding to the spherical case, which suggests that the elongated geometry exhibits improved lasing properties over the spherical one. This is confirmed by comparing the results of insets of Figure a and b, which show the steady-state values of ⟨|Ee(r, t)|2⟩ as a function of |Epump|2 for Lrod = 0 and Lrod = 60 nm, respectively. The observed linear dependence above threshold confirms that the two considered configurations are indeed lasing. In addition, from the displayed results it can be deduced that the nanoparticle with Lrod = 60 nm features a slope efficiency factor of approximately 300 times larger than that obtained for Lrod = 0 (we obtained values for the slope efficiencies of 2.1 × 10–24 and 6.5 × 10–22 m3 for Lrod = 0 and Lrod = 60 nm, respectively). At the same time the lasing threshold is reduced by a factor of about 24 in the elongated case with respect to the spherical one (values of 0.5 × 1012 and 11.9 × 1012 V2/m2 were obtained for the elongated and the spherical cases, respectively). A detailed analysis of the variation of these enhancement factors with Lrod is provided below.

Of special interest is the comparison of the above numerical results for the spherical case with those reported experimentally by Noginov et al.[13] We found that for this particular spaser a density of the molecules of at least 5 times larger than the one used in ref (13) is required in the shell of the structure in order to reach lasing action at the pump intensities levels considered in that work (the number of dye molecules per nanoparticle estimated in ref (13) is 2.7 × 103, which, assuming that the molecules are uniformly distributed over the silica shell, yields a concentration of N0 = 6.3 × 1025 m–3). This is clearly visualized in the inset of Figure a, which corresponds to pump intensities similar to those considered in ref (13), but a concentration of dye molecules of N0 = 3.7 × 1026 m–3. Equivalently, we did not observe any signature of laser action at those pump intensity levels for the value of the concentration of dye molecules reported by Noginov and coauthors. This discrepancy can be ascribed to the fact that an additional optical feedback mechanism, beyond the one associated with the LSP supported by a single nanoparticle (such as the reillumination among neighboring nanoparticles present in the experimental realization), could be playing an important role in the experimental laser observations reported in ref (13). Alternatively, the active molecules in the experimental configuration could be not uniformly distributed over the whole silica shell volume, but concentrated in a layer of smaller volume close to the metallic core. Finally, we also note that the difference between our simulations and the experimental data of ref (13) cannot be accounted for in terms of the so-called Purcell effect.[58] This effect would actually lead to an effective reduction of the available gain close to the metallic surfaces of the system,[53] which, in turn, would increase the difference between theoretical predictions and the experimental results. Instead, below we propose a novel route, based on tailoring the shape of the nanoparticles, which allows reaching lasing action at (or even below) the above-mentioned pump intensity levels using dye molecule concentrations as low as those used in ref (13).

To further investigate how lasing action is realized in the considered class of systems, we study how the above-described lasing dynamics is linked to the population inversion of the corresponding lasing transition. For definiteness, we focus on the largest pump intensity considered in Figure b (|Epump|2 = 25|Eth|2). Figure a displays an enlarged view, between t = 0 and t = 200 ps, of the time evolution of the calculated ⟨|E(r, t)|⟩ (gray line, right y-axis), along with the dynamics of the corresponding spatially averaged population inversion ⟨ΔN(r, t)⟩ (normalized to N0; see red line, left y-axis). As expected, before the first lasing spike occurs, the averaged population inversion grows almost linearly with time. This corresponds to the regime in which the population of the upper level of the relevant lasing transition is increasing (the system is accumulating population inversion), and the whole system effectively behaves as an optical amplifier. This first operating regime lasts until t ≈ 50 ps. At that time the population inversion becomes large enough so its associated optical gain can overcome all the losses present in the system (both radiative and ohmic). As observed in Figure a, it is then when the first burst of laser light generation takes place. This burst, in turn, leads to a significant depletion of the population inversion (a significant amount of the upper-level population of the laser transition decays via stimulated emission), leading to a dramatic drop of the laser signal. After that, it starts a subsequent recovery of the population inversion, until, again (when enough population inversion is accumulated), a second spike of the laser signal occurs, accompanied by the corresponding drop in the population inversion. This series of bursts and subsequent drops of the population inversion takes place sequentially (for larger times smoother spikes and drops of the lasing signal and population inversion are obtained) until the steady state of the laser is finally reached.

Figure 4

(a) Time evolution of the spatially averaged population inversion density, ΔN (normalized to the total density of dye molecules, N0), for the spaser considered in Figure b, as computed for |Epump| = 5|Eth| (red line, left axis). For comparison, the dynamics of the corresponding spatially averaged lasing field ⟨|E(r, t)|⟩ (normalized to the saturation electric field, E0) is also included (right axis, gray line). (b) Spatial distribution of the population inversion inside the active shell computed at four representative times (labeled as tA, tB, tC, and tD in both panels).

The microscopic origin of the above-described macroscopic behavior can be understood by examining the spatial profile of the corresponding population inversion distribution, ΔN(r, t). Figure b shows cross sections (along the zx-plane) of ΔN(r, t) for the configuration considered in Figure a, as computed at four representative values of t (labeled as tA, tB, tC, and tD in Figures a and b). As expected, before the onset of lasing action (t = tA), the population inversion follows the profile of the longitudinal LSP supported by the nanoparticle (see rightmost panel of Figure a). This distribution starts changing quickly once we enter into the lasing regime. Specifically, after the lasing onset, the regions of high electric-field intensity are the ones that experience faster depopulation of the lasing transition (these regions feature enhanced stimulated emission rates), while regions with low intensities retain most of the population inversion they accumulated before lasing action starts. This is precisely what can be observed in the sequence of ΔN(r, t) for t = tB, t = tC, and t = D displayed in Figure b. As shown, as time grows, the population inversion distribution starts increasingly resembling the complementary profile of the above-mentioned LSP field profile (we find maxima of the population inversion at the minima of the field and vice versa). This particular population inversion distribution settles down in the steady state, giving rise to a subwavelength-scale highly nonuniform gain distribution whose spatial average leads to the effective lasing response of the system for long times. This behavior can be seen as a novel instance at the subwavelength scale of the well-known lasing processes observed in traditional active systems.[56]

Next, we study how the steady-state characteristics of the system depend on the elongation of the nanoparticle. Figure a renders the simulation results for the lasing threshold (blue squares; left y-axis) and the slope efficiency (red squares; right y-axis) as a function of Lrod (the rest of the geometrical and gain medium parameters are the same as those used in Figure ). Both magnitudes have been normalized to the corresponding values for the spherical configuration. As seen in Figure a, the normalized lasing threshold decreases rapidly from 1 to 0.06 when the elongation of the nanoparticle is increased from Lrod = 0 to Lrod = 20 nm. However, for values of Lrod greater than 20 nm, the lasing threshold becomes much less sensitive to elongation. Interestingly, a closer examination of the threshold for large Lrod reveals a nonmonotonic dependence, leading to an optimal elongation value of Lrod ≈ 50 nm (see inset of Figure a). On the other hand, the normalized slope efficiency also displays a rapid variation with Lrod (the slope efficiency grows with Lrod as shown in Figure a), followed again by a weak nonmonotonic dependence on the elongation. In this case, the plateau region and the location of the optimal value for Lrod occurs at larger values of elongation than those obtained for the lasing threshold (an optimal value of Lrod ≈ 70 nm is obtained for the slope efficiency). Despite this difference, the remarkable broad range of Lrod values in which both the slope efficiency and the lasing threshold are near their optimal values allows for the existence of an interval (defined by 60 nm < Lrod < 80 nm) in which the two magnitudes are simultaneously improved by factors close to 300 and 24, respectively. We believe this result could have a significant influence for further engineering and optimization of spasers based on metallic nanoparticles.

Figure 5

(a) Simulation results for the steady-state values of the lasing threshold (blue squares, left axis) and the slope efficiency (red squares, right axis) as a function of the nanoparticle elongation Lrod. The same volume of the active shell is assumed for all Lrod values. The rest of geometrical and gain medium parameters are the same as in Figure . Both the lasing threshold and the slope efficiency are normalized to the corresponding values obtained for the spherical configuration (Lrod = 0). The inset shows an enlarged view of the normalized lasing threshold between Lrod = 20 nm and Lrod = 100 nm. (b) Coupled-mode theory results for the normalized slope efficiency calculated for the same configurations as in Figure a (red squares). The panel also includes the results for two alternative conditions to the same-volume condition applied to the shell: green squares correspond to assuming the same shell thickness for all Lrod values, whereas black squares correspond to imposing the same molecule number for all elongations. (c) Same as (b) but now for the normalized lasing threshold. Inset of panel (b) shows the calculated quality factor of the lasing localized surface plasmon mode as a function of Lrod, whereas the inset of panel (c) displays the corresponding value for the spatial confinement factor of that mode, ξ1 (normalized to its maximum value, ξ1), as a function of elongation.

Much in the same way as occurs in conventional macroscopic laser sources,[59] we expect that the above-described numerical results can be qualitatively understood in terms of the temporal and spatial light confinement properties of the studied structure. In particular, we expect the quality factor of the longitudinal LSP supported by the nanoparticle (describing the temporal confinement properties of the system) and the fraction of field energy residing in the gain medium (the so-called confinement factor in conventional laser rate equation approaches[59]) are the key parameters to account for the evolution of the system lasing characteristics with elongation. To obtain specific analytical expressions of the dependence of both the laser threshold and the slope efficiency on these two magnitudes, we apply a coupled-mode theory (CMT) analysis to the problem. This framework, based on a first-order perturbation-theory approach to Maxwell’s equations,[60−62] can provide analytical expressions of the lasing characteristics of plasmonic lasers. Specifically, following refs (61) and (62), we can write the slope efficiency (S) and the lasing threshold (Rth) of the considered spasers as S = 2ℏωeN0τ(ξ1/ξ2) and Rth = 2α/(N0τξ1). Here τ = 2Q/ωe is the lifetime of the lasing mode, whereas α is a parameter that contains the details of the lasing transition, α = Γe/(τ21Ke) (where Ke is the coupling strength of the polarization to the local electric field; see Methods). The parameters ξ1 and ξ2 quantify the degree of spatial confinement of the lasing mode; ξ1 = (1/2U)∫ dr |E(r)|2 and ξ2 = (2/U)(∫ dr |E(r)|4/∫ dr |Ee(r)|2), where Ee(r) is the spatial profile of the mode, U represents the corresponding total EM energy, and VA is the region of the system where the gain medium resides. Note that ξ1 represents the fraction of the total field energy contained in the active region of the structure (it can be seen as the CMT counterpart of the so-called confinement factor appearing in conventional laser rate equation approaches[59]).

Red squares in Figure b and c display respectively the CMT results for the normalized slope efficiency and normalized lasing threshold calculated for the same configuration as in Figure a. As seen, although CMT cannot accurately reproduce the actual values obtained in the full simulations, it does capture the overall dependence of the lasing characteristics on Lrod (including the location of the optimal elongation value in each case). Moreover, we found numerically that the ratio ξ2/ξ1 decreases monotonically as a function of Lrod. In view of the expression for S provided above, this result implies that the observed nonmonotonic dependence of the slope efficiency on the nanoparticle elongation relies on the evolution with Lrod of the quality factor Q. This is further corroborated in the inset of Figure b, which shows the calculated Q values as a function of Lrod.

On the other hand, from the fact that Rth is inversely proportional to the product ξ1Q (see analytical expression above), we deduce that both spatial and temporal light confinement properties play an equally important role in reducing the lasing threshold as the nanoparticle elongation increases. The calculated dependence of ξ1 on Lrod is rendered in the inset of Figure c. As observed, the rapid reduction of the lasing threshold with Lrod is actually the result of the joined growth of ξ1 and Q for small Lrod values, whereas the presence of an optimal Lrod value now emerges as a consequence of the combination of the nonmonotonic behaviors of Q and ξ1 as a function of elongation. This analysis also provides insight into the dependence on Lrod of the lasing threshold for large elongation values. In particular, as deduced from the above analytical expression, the dependence of Rth on the nanorod elongation enters through the corresponding dependence on Lrod of the product Qξ1/ωe (the rest of the parameters in the analytical expression of Rth do not depend on Lrod). This shows that for large Lrod values the trend displayed in Figure c (red squares) is the net result of the combination of the temporal and light confinement properties of the structure, along with the associated shift in frequencies of the considered LSP resonances. Finally, we note that the above CMT generalizes previous analytical approaches based on quasi-static analytical descriptions.[15]

We have also applied CMT to two cases in which the constraint of maintaining the volume of the active shell for all Lrod values is replaced by two alternative conditions that could be of relevance from a experimental standpoint. Namely, green squares in Figures b and c show the numerical results for the slope efficiency and the lasing threshold computed by fixing the thickness (instead of volume) of the active shell for all considered Lrod. Black squares in the same panels correspond to similar CMT calculations but now assuming the same number of dye molecules for all elongation values (the concentration of molecules is the same in all cases). Despite the relative variation in the magnitude of the slope efficiency and the lasing threshold (which can be ascribed to the difference in gain medium volume among configurations), the overall trend of the lasing characteristics is maintained, which supports the general character of our findings.

Conclusions

In conclusion, by using detailed simulations based on a time-domain generalization of the FEM method, we have analyzed the spatio-temporal dynamics of lasing action in spasers based on core–shell metallic nanoparticles. We have particularly focused on studying how the lasing characteristics of this class of structures are influenced by the nanoparticle shape. We have found that both the lasing threshold and the slope efficiency of conventional spherical spasers can be significantly improved simply by elongating the nanoparticle. In this context, we have also found that the enhancement of the laser characteristics is maintained across a broad range of elongation values. Moreover, we have used an analytical coupled-mode theory to explain these findings in terms of the spatial and temporal light confinement properties of the LSP modes supported by the nanoparticles. We expect this work to stimulate further theoretical and numerical investigations on laser light generation assisted by localized surface plasmons and particularly on how nanoparticle shape optimization can be used to tailor lasing emission at the subwavelength scale. In particular, future work includes the extension of the reported time-domain semiclassical approach to a fully quantum-mechanical model able to account for quantum confinement and electron mean free path effects. Due to the versatility of metallic nanoparticles in a number of contexts beyond lasing, we believe this work could be of relevance across a broad spectrum of different areas, including molecular sensing, photovoltaics, nanoscale microscopy, and emerging quantum technologies.

Methods

To compute the optical response of the considered active plasmonic systems, we use a semiclassical numerical framework based on solving the field equation for the vector potential A(r, t),where P = P(r, t) is the time-dependent part of the polarization (the time-independent contributions to the polarization are accounted for through the parameter ϵr). Within the metallic regions of the system, P(r, t) is computed through the conventional Drude–Lorentz form. Inside the active medium, P(r, t) has two different contributions, P(r, t) = Pa(r, t) + Pe(r, t). These contributions arise from the stimulated absorption (Pa(r, t)) and emission (Pe(r, t)) of photons in the gain medium. In the particular case of a gain medium including four-level quantum emitters (with absorption and emission electronic transitions centered at ωa and ωe, respectively), the dynamics of P(r, t) (i = a, e) is governed by the following Lorentzian equation:with Γ and K being respectively the linewidth of the corresponding transition and the coupling strength of the polarization P to the local electric field, E = −∂A/∂t. The functions ΔN = ΔN(r, t) represent the population inversion densities of the absorption and emission transitions, ΔN(r, t) = N3(r, t) – N0(r, t) and ΔNe(r, t) = N2(r, t) – N1r, t). N(r, t) (with j = 1, ..., 4) are the population densities of each energy level of the quantum emitters. The temporal evolution of these populations is determined by the following rate equations:where τ are the nonradiative decays between the jth and kth energy levels and ⟨·⟩ denotes the time average over a period Te = 2π/ωe.

Thus, by solving the coupled set of nonlinear equations given by eqs –6, we obtain the whole spatio-temporal dynamics of the studied systems, including all the characteristics of their eventual laser emission. In order to reduce the computational requirements of this problem, we introduce the two additional steps summarized in the main text. First, we exploit the fact that two characteristic frequencies of the system are known (namely, ωa and ωe). This allows us to write A(r, t) = Aa(r, t) exp(−iωat) + Ae(r, t) exp(−iωet), where A(r, t) are slowly varying complex amplitudes. We also assume that a similar expansion holds for P(r, t). These expansions enable tracing the fast optical oscillations out of the problem, and, consequently, reduce the total simulation time. To allow for the building-up of energy at ωe, an initial seed excitation (of amplitude much smaller than that of the pump) is assumed in the temporal profile of Ae. We have checked numerically that our results do not depend on the amplitude or duration of that seed excitation.

Second, we rewrite eq in the so-called weak form. In order to do that, we multiply the field equation governing A(r, t) by a test function Ã, then we integrate the result over the simulation domain Ω, and finally we perform the resulting integral by parts. This yields[55]where ∂Ω represents the boundary of the simulation domain and n is a unitary along the direction normal to that boundary. Note that the physical boundary conditions of the studied problem are imposed through the last term of the right-hand side of the above equation. By construction, this weak-form formulation makes the considered numerical problem ideally suited to be solved with a FEM method.[47] This enables bringing the benefits of the adaptative nonuniform meshing inherent to FEM methods to time-domain simulations of active plasmonic systems.