Few-emitter lasing in single ultra-small nanocavities

Abstract Lasers are ubiquitous for information storage, processing, communications, sensing, biological research and medical applications. To decrease their energy and materials usage, a key quest is to miniaturise lasers down to nanocavities. Obtaining the smallest mode volumes demands plasmonic nanocavities, but for these, gain comes from only a single or few emitters. Until now, lasing in such devices was unobtainable due to low gain and high cavity losses. Here, we demonstrate a form of ‘few emitter lasing’ in a plasmonic nanocavity approaching the single-molecule emitter regime. The few-emitter lasing transition significantly broadens, and depends on the number of molecules and their individual locations. We show this non-standard few-emitter lasing can be understood by developing a theoretical approach extending previous weak-coupling theories. Our work paves the way for developing nanolaser applications as well as fundamental studies at the limit of few emitters.

Lasing occurs when stimulated emission into a cavity mode exceeds loss, leading to amplification.Typically, this causes a sharp change of slope in the emission vs input power.Such a sharp transition is analogous to a thermodynamic phase transition [10].As with phase transitions, there exists a 'system-size' parameter β, which determines how sharp the transition is [5].However, questions remain about the sharpness of this transition for lasers with few emitters, small cavities, and stronger light-matter coupling.Such questions are particularly important for nanocavities which confine light within sub-wavelength volumes [11], such as metasurfaces or plasmonic nanocavities that exploit collective electron oscillations in metallic nanostructures to achieve extreme confinement (V < 100 nm 3 ) [12,13].These enable coupling single emitters to light [14][15][16].Lasing in such plasmonic nanocavities presents new opportunities for miniaturisation and integration, but also raises new questions about the nature and conditions required for lasing in a regime which combines a small number of strongly-coupled emitters with lower quality resonators, Q ∼ 10.Our aim here is to understand this regime.
The sharpness of the lasing transition reflects how lasing enhances the conversion efficiency of input power into output light.Above threshold the efficiency is high, as stimulated emission directs almost all radiation into the cavity mode.A sharp transition requires low efficiency below threshold.This is captured by the parameter β, the ratio of input-output slopes below and above threshold.Small β indicates a sharp transition.When the light-matter coupling g is weak, β = g 2 /(g 2 + Γ ↓ Γ T ), where Γ T is the total emitter linewidth, and Γ ↓ the decay rate into non-cavity modes [5].Below threshold, β is the fraction of incoherent emission into the cavity and in this weak-coupling case, it fully determines the shape of the input-output curve.It is notable that β does not depend on the number of emitters, whereas the parameter determining the sharpness of thermodynamic phase FIG. 1. Plasmonic nanocavity with emitters.a, Nanoparticle-on-mirror (NPoM) cavity, formed by Au nanoparticle trapping few emitters in 0.9 nm gap above Au film.b, Nanogap region with N emitters of dephasing Γz excited by pump strength Γ ↑ , emitting into cavity with loss rate κ.c, Re-entrant lasing threshold, shown as colormap of photon number (yellow=high, purple=low) vs normalized coupling and pump strength, using κ = 0.1Γ ↓ , N = 10, and Γz = 10Γ ↓ .Dashed white line shows threshold defined by Eq. (3).d, Cross sections, as indicated by horizontal dashed lines in (c).At large coupling the traditional input-output curve is seen (red line) while at smaller coupling lasing is suppressed at higher pump power (blue line).
transitions is typically the system size.As we discuss below, allowing for strong light-matter coupling changes this behaviour considerably.Some extreme limits of lasing have previously been explored.The single-emitter limit [8] has been studied with atoms [17], superconducting circuits [18], and quantum dots [14,19].In this limit, the transition broadens as β becomes large.This is because lasing with N emitters requires N C > 1, where the cooperativity C = 4g 2 /(κΓ T ) depends also on the cavity loss rate κ.When N = 1, one needs C > 1, and so a sharp lasing threshold (β 1) is only possible if κ Γ ↓ , which is not typically the case.In this paper we explore lasing of a few organic molecules coupled to a plasmonic nanocavity.Despite low Q, emitters coupled to plasmonic modes can achieve lasing [2].However the smallest lasers must lase with only a few emitters, a goal so far thwarted, but attainable by using ultrasmall volume plasmonic nanocavities.These can be realised using bottom-up self-assembly, which we acheive via the nanoparticle-on-mirror (NPoM) geometry (Fig. 1a) of a Au nanoparticle on a thin dielectric spacer above an Au film [13].We trap light inside 1 nm-high ∼20 nm-wide gaps, which enhance incident optical fields by >250 while retaining ∼ 50% radiative efficiency.Enhancing both the pumping and the light-matter coupling now enables lasing, confirmed by enhanced coherence in the lasing state.To describe such experiments we must extend previous theoretical treatments [5] to consider the combination of stronger coupling, few emitters, bad cavities, pumping-induced noise, and inhomogeneity.
We model the light-matter interaction in this NPoM nanocavity using a standard master equation describing many two-level molecules coupled to light with strength g i .The molecules undergo incoherent pumping, loss, and pure dephasing at rates Γ ↑ , Γ ↓ , Γ z (see Methods, Fig. 1b).We then make a second-order cumulant expansion [4], giving coupled equations for photon number n, molecule-photon coherence, inter-molecular coherence, and molecular inversion.This approach is ideal for understanding how system size N controls the sharpness of the transition, since it captures finite size effects treating N as a parameter.It also includes the effects of spontaneous emission, and recovers the semiclassical theory of lasing [5,20] in appropriate limits.In contrast to weaklycoupled models of incoherent emission and absorption [5] our approach captures the back-action of the coherent light on the dynamics of the emitters.
Considering first the homogeneous case where g i = g, we solve for steady-state n with constant pump Γ ↑ giving with n 0 defined in the methods, and where the total linewidth is Γ T = Γ ↑ + Γ ↓ + 4Γ z .While Eq. ( 1) appears identical to that found in [5], a subtle and crucial difference exists: the terms Γ c and β both implicitly depend on Γ ↑ , because the cooperativity C depends on Γ T .This Γ ↑ dependence can often be ignored, when Γ T 4Γ z .However in the regime where the collective cooperativity is low (as here), it is important to keep this dependence.Physically, strong pumping introduces phase noise which ultimately kills lasing by suppressing the cooperativity.For typical lasers this effect only occurs at very high powers, however for plasmonic nanocavities with high losses, a switch-on and switch-off transition are seen, and these can merge where N C 1 [21].This is seen in Fig. 1(c,d), which show how the photon number depends on the normalised pump power and coupling strength.
A second crucial difference from [5] is the form of β in Eq. 2. In particular, this now depends on N , with β ∼ 1/N at large N .This encodes the anticipated effect that larger systems show smaller finite-size effects.
Experimentally, NPoMs with ∼0.9 nm gap are defined by a monolayer of cucurbit [7]uril molecules that can each encapsulate a single molecule of methylene blue [15] (for sample preparation see Methods).We vary the expected number of emitters N in the nanocavity by changing the ratio of emitting methylene blue to non-emitting cucurbit [7]uril molecules.Because the location of the emitting molecules is different in each NPoM, their overlap with the nanocavity mode and hence coupling g i varies.Using NPoM nanoparticle diameters of 60 nm tunes the dominant NPoM cavity mode to the gain peak (Fig. 2a), optimising the interaction of the emitted light with the nanocavity.
To create enough gain, ultrafast pulses (100 fs, 640 nm) are used to irradiate individual NPoM cavities whose emission is recorded by a low-noise spectrometer (for experimental setup see Methods).At average powers of P av = 4 µW (power density = 3.2 kWcm −2 ), the emission spectra for different numbers of molecules is similar and the total integrated intensity varies linearly with N (Fig. 2b,c).
Power dependent measurements on multiple NPoM constructs are performed with different N (Fig. 2d-g).For NPoM cavities with lower N , smooth transitions in the total emitted intensity are seen, becoming systematically sharper for higher N .In particular for N ∼ 2, no observable transition is seen and the total emission is linear before saturating at the highest powers.This is consistent with thresholdless lasing, approaching the single-emitter regime.For NPoMs with large N , wellpronounced thresholds are seen with super-linear emission, which then saturate to linear scaling as in conventional lasers.We fit Eq. ( 1) to the experimental data (for fitting procedure see Methods); our theory matches well across the range of N , indicating lasing in plasmonic nanocavities with a few emitters [22].Furthermore, the extracted β-factor varies inversely with N (Fig. 2h), as expected from our theoretical model (Eq.2).
The analytic form of Eq. ( 1) ignores the effects of random placement of molecules, which gives variable lightmatter coupling.This variation becomes important when the number of molecules is small, as it leads to significant variation between different NPoMs with the same number of molecules.With random positions one finds that the shape of the input-output curve (and thus the parameter β) depends on the distribution of couplings, as shown in Fig. 2(i,j).
Because the threshold is less defined at small N , it becomes important to find further experimental evidence for lasing.To this end, we examine the spectral evolution and coherence across the transition.As the average power P av increases from 4 µW to 250 µW (Fig. 3a), the emission (which is 10 nm red-shifted from the solution PL due to the cavity environment) significantly broadens, with increasing weight at shorter (bluer) wavelengths.This behaviour appears consistently in most NPoMs at all N , with emission even extending to the blue side of the pump laser wavelength (see SI).One explanation for this change in emission spectrum with power is a breakdown of Kasha's rule [23].Kasha's rule states photoluminescence is only from the lowest energy vibrational state, as excited states relax before luminescence.However, radiative decay is drastically sped up in NPoMs, reaching below 50 fs as the Purcell factor is ∼ 10 5 due to the large Q/V .When emission becomes stimulated, this can become faster than the vibrational lifetime (∼ 1 ps) [16], allowing direct emission from higher vibronic states, breaking Kasha's rule [24].
To track the coherence of the observed emission we explore spectral interference using a Michelson interferometer (Fig. 3b).The delay time τ between the two arms of the interferometer is greater than the temporal coherence time (τ > τ c ∼ h/Γ T ∼ 13 fs).Spectral fringes are observed both below and above threshold (Fig. 3) with a wavelength period 7 nm inversely proportional to τ .The extracted visibility from the spectral envelope of the fringes, decreases with wavelength below P th to 0.1 at 800 nm (shown for two different time delays τ = {33, 100} fs, in Fig. 3c).By contrast above threshold, visibility increases to 0.8 becoming constant over the entire spectrum.
As discussed in Ref. [25], by combining spectral filtering with coherence measurements, we show in Fig 3c,d how the spatial coherence of the light changes through the threshold.To better understand this increasing visibility, the far-field emission is measured in Fourier space (for optical setup see Methods).The dominant bright plasmonic mode S 10 emits at high angles both below and above threshold (seen as bright, purple-coloured, ring in Fig. 3e), around 40 − 50 • (Fig. 3f).Additional weaker emission at low angles-arising from the S 11 mode [26]is proportionately larger below threshold.This mode mostly outcouples light from molecules near the edges of the nanoparticle facet (Fig. 3f), but only with poor quantum efficiency < 5% [26].We define the emission ratio r = I ring /I centre between the intensity from the annular ring I ring and from the centre I centre .This measures the relative emission of these modes.By comparing the emission ratios above and below threshold r a,b , we see that r a /r b ∼1 for some NPoMs but >1 for many others.This indicates preferential occupation of one mode, as is typical in lasers, due to stimulated emission.
In summary, we show extreme confinement of optical fields exciting a few molecules in a nanometer gap that reveals lasing is possible with a few emitters, and that this can match only extended models.This detailed influence of cavity Q allows for future work exploring the photon statistics from such ultimately miniaturised lasers, and clarifying the fundamentals of the lasing nonlinear phase transition.

A. Theoretical modelling
In this section we describe the theoretical approach we use to calculate the input-output curve for lasing with few strongly-coupled emitters.We first introduce the model, then discuss the cumulant approach, and how the resulting equations can be solved with and without disorder in the matter-light coupling.

Model
We model the system as a plasmonic mode coupled to a collection of N two-level systems: here, a ( †) is the annihilation (creation) operator for the plasmonic mode of energy ω while σ α i are the Pauli matrices representing the two-level system at site i.We use a Lindblad master equation to account for the following incoherent processes with their corresponding rates: photon loss (κ), incoherent pump (Γ ↑ ), non-radiative losses (Γ ↓ ) and pure dephasing (Γ z ).We thus have with D[X] = XρX † + 1 2 X † Xρ + ρX † X .The model we consider does not explicitly include the vibronic progression present in organic molecules [4].The incoherent pumping rate Γ ↑ we consider should however be understood as an effective rate, arising from the combination of coherent pumping exciting a higher vibrational state followed by vibrational relaxation (see Fig. 1).This simplified model is chosen since, as discussed below, it allows for closed-form expressions for the lasing threhsold.Exploring how more complex treatment of the molecular spectrum affects these results is an important question for future work.

Cumulant equations
Various approaches exist to model a system described by the above equations [27].In this work, as discussed in the main text, our focus is on understanding how system size N enters into determining the sharpness of the transition.For this, it is useful to choose an approach which captures finite size effects, treats N as a parameter, and captures the semiclassical effects of spontaneous emission.The ideal approach for this is to use second order cumulants [4,28], which provide a systematic expansion in 1/N .
From the density matrix equation of motion given above, we can write down equations of motion for the second order moments of the system, using cummlant expansion to break higher order moments into a combinations of first and second order moments.The non-zero moments are: n = a † a , P i = Im aσ + i , S i = σ z i and where we require the two operators in the last moment to act on different molecules.Assuming the resonant case ω = , the equations of motion end up taking the form: where Γ T = 4Γ z +Γ ↑ +Γ ↓ is the total molecular linewidth and S 0 = Γ ↑ +Γ ↓ is the inversion set by the pump.Note that in writing these equations, we have allowed molecule-dependent coupling strengths g i .

Homogeneous limit
We will start by considering the case where g i = g, i.e.where all the molecules couple identically to the light mode.If we use Eqs.(6)(7)(8)(9) to adiabatically eliminate everything but the photon occupation n we end up with the quadratic equation for the steady state: If we make the approximation that Γ T + κ 1 − 1 N ≈ Γ T + κ, we can write the solution on the form Introducing the cooperativity C = 4g 2 /κΓ T we can write the quantities appearing in Eq. ( 1) as: As noted in the main text, the dependence of Γ T on Γ ↑ means the input-output curve takes a more complex form than discussed by [5].

Effects of inhomogeneity
In reality the electric field strength at the location of the molecules (and thus their coupling to light) will not be identical.This becomes especially important when the molecules are few and the mode volume is very small, as is the case in this paper.We next discuss the effects of such disorder.To do this we reintroduce an i-dependence of the coupling strength g i in the cumulant equations.In this case, to find the steady state, it is convenient to first solve for P i .This must satisfy the quadratic equation 0 = A i P 2 i + B i P i + C i with coefficients that depend on the other P j via the moments Π (n) = i g n i P i .The coefficients take the form: For a given set {g i } we can solve these equations iteratively.We first guess an initial {P i }, then evaluate the moments Π (n) , n = 0, 1, 2 to find the coefficients A i , B i and C i .We then solve the quadratic equation to give a new set {P i }.This process can then be iterated to convergence.We find it necessary to use successive over relaxation to improve the stability of this.In Fig. 2(i,j) we show the result of this, where we choose a set of g i drawn from a normal distribution of a given mean and variance.In plotting these we in fact adjust the mean of this distribution so as to keep the root mean squared g constant, sincemthe behaviour of strongly-coupled systems typically depends on i g 2 i .Figure 2(i) shows a set of input output curves; one sees that with disorder, each realisation leads to a slightly different curve.Figure 2(j) is a histogram showing the frequency of slope ratios β extracted from these curves, which is therefore the resulting probability distribution of this parameter.

B. Experimental setup
The details of the setup can be found in Ref. [16] and we give a brief description here.For the power dependent measurements, we used a variable neutral density filter to control the intensity of 100 fs pulses at wavelength 640 nm generated from a tunable optical parametric oscillator (SpectraPhysics, Inspire), pumped (Spect-raPhysics, Maitai) at 820 nm with a repetition rate of 80 MHz.A brightfield/darkfield microscope objective (Olympus, 100×, numerical aperture = 0.9) focused the attenuated pulses on a single NPoM nanocavity to excite the emitters.The emitted light was collected in the reflection direction, filtered, and measured using a grating spectrometer (Andor, Shamrock 303i) and an electron multiplying charged-coupled detector (EMCCD, Andor iDus).
The farfield Fourier space imaging was performed by demagnifying the collected emission by the objective and then focussing it in the focal plane of a lens before the entrance of a monochromator slit.Darkfield spectroscopy was performed on a home-built microscope that illuminates the sample with white light from a halogen lamp (Philips, 100W) at angle of 60 • through the darkfield microscope.The scattered light was collected through the microscope objective and the scattered spectrum was measured using a fibre-coupled grating spectrometer (OceanOptics).

C. Sample preparation
The sample substrate is a Au mirror on a Si wafer that was fabricated using a template-stripping process.Thermal evaporation was used to coat a Si wafer with a 100 nm thick Au film at an average evaporation rate of 0.5 Angstrom/s.Small pieces of another Si wafer (area 4×4 mm 2 ) were glued to the evaporated Au using epoxy (Epo-Tek 377), left overnight at 150 • C and then slowly cooled down to room temperature.On applying a slight force, the silicon pieces easily peeled off and were covered by an atomically-flat Au film.A 1 mM solution of methylene blue (Sigma Aldrich) and a 1 mM solution of cucurbit [7]uril (CB) (Sigma Aldrich) were mixed together, thus allowing the encapsulation of the methylene blue guest molecules inside the CB cavities.A freshly stripped substrate with Au film was submerged in the solution overnight, thoroughly rinsed with de-ionized water, and then blown dry with nitrogen, leaving behind a self-assembled monolayer (SAM) of CB encapsulating a single molecule of methylene blue [15].Au nanoparticles (BBI Solutions, diameter 60 nm) were drop cast onto the SAM on the Au substrate for 10 s.The excess Au nanoparticle solution was rinsed with water and blown dry using nitrogen.This leaves sparsely spaced Au nanoparticles deposited on the CB layer on the Au film, forming the NPoM nanocavity.To control the number of emitters in the nanogap, we diluted the 1 mM methy-lene blue solution by factors of 2 to 20 and kept the CB concentration fixed at 1mM.Using the surface packing density of CB of 0.24 molecules.nm 2 [26] and a facet size of 15 nm, we estimated the average number of molecules in the nanogap to vary between 2 and 50 for different dilution ratios.

D. Fitting of theoretical curves to experimental data
In Fig. 2(d-g) we present results of fitting the theoretical expression Eq. ( 1) to the experimental input-output curves.Here we describe the approach used to this fitting.
To perform such a fit, one must relate the theoretical input Γ ↑ and output n to the experimental power in and power out.This requires a description of how light couples into and out of the coupling.The in-coupling is taken to be proportional to the coupling strength squared, Γ ↑ ∝ g 2 ×Power density.This is because the external pump couples to the emitters through the same dipole moment and mode profile as controls the lightmatter coupling.As we assume this incoupling is an incoherent absorption process, Fermi's golden rule implies a rate proportional to g 2 .Including this effect was crucial in producing an accurate fit; without this factor our fitting led to spurious correlations between fitting parameters.Out-coupling is simpler, as this involves how the light in the cavity couples to the far-field.This was assumed to be independent of other parameters and considered as an intensity scaling factor.
With the above assumptions, we then performed a standard least squares fit to find the remaining model parameters.Some of these parameters are taken as global parameters-common to all experiments; these were Γ z , Γ ↓ , κ and were first optimized (see SI).The remaining fitting parameters for each NPoM are g, an intensity scaling factor, and N which only varies by 20%, due to variations in the number of emitters.The results of the fittings show that there is no correlation between the fitting parameters and that the out-coupling efficiency is linearly proportional to the output intensity, as expected (for more details, see SI).NPoM does not show any corresponding peaks between the SERS and emission (Fig S4 ), implying that the ground state vibrational levels are not involved in the process.It would therefore be interesting in future to investigate the impact of the various vibronic states on the light-matter interactions in this lasing regime [31].

C. Fitting results
As described in the Methods section, we fit the theoretical result to the experimental data to determine relevant parameters.We divide parameters into three classes: A. Parameters which we fix globally over all NPOMs.In this class is Γ z , which is assumed an intrinsic property of the molecules.
B. Parameters which are constant over NPOMs but allowed to vary according to the concentration of emitters (and thus target N ).In this class are Γ ↓ , κ.
C. Parameters which vary between each NPOM.In this class is the value g and an output intensity scaling which we denote S out .This output scaling relates the measured output counts C out to the theoretical photon number n vias C out = S out n.In addition, we allow the value N to vary, within 20% of the target N , to account for variation in the number of emitters actually coupled.
and finding the area in parameter space defined by this inequality allows us to estimate errors.Here, k eff and ∆ are the reduced number of values and p-value from a χ 2 distribution table.This area is shown by grey contour in Fig. S5.
In performing this fitting, we take the theoretical parameter Γ ↑ = g 2 × (P owerdensity), where Power density is the measured value in units of Wcm −2 .The theoretical behaviour depends only on dimensionless ratios of the parameters Γ ↑ , Γ ↓ , Γ z , κ, g, in order to produce the dimensionless photon number n.As such, the somewhat unusual units for Γ ↑ do not cause a problem.Physical units could be restored by introducing an extra scaling parameter between power density and Γ ↑ , however since the theoretical fits depend only on ratios of parameters, the fit would be degenerate in this parameter.
A subset of the measured input-output curves appear to show different behaviour, with a saturation of the output power occurring below the lasing threshold.As such, these datasets cannot be fitted to the theoretical model.The origin of this different form of behaviour may be physical or instrumental, which needs further investigation.The fit to these datasets are considered invalid, since the extracted fitting parameters give erroneously high values.Thus, we consider only fits where χ 2 R ≤ χ2 R .These fits account for < 30% of the entire dataset shown in Fig. S2.There is no correlation between the fitting parameters in the class C [Figs.S6(a-c)], for χ 2 R ≤ χ2 R .This lack of correlation is expected, since the coupling between emitters and cavity mode, and the out-coupling of the cavity mode to the detector should be independent parameters.There is however a correlation between the outcoupling intensity scaling and the emission below threshold [Figs.S6(d-f)], which implies that an efficient outcoupling efficiency of NPoM cavity results into a high emission intensity.

D. Shape of input-output curve
As mentioned in the main text, the shape of the input-output curve is not fixed when one goes beyond the weakcoupling laser theory of Ref. [5].In the weak coupling limit, the parameter β fully determines the shape of the curve, while other parameters merely rescale the photon number or input power.In contrast, for strong couplingparticularly combined with inhomogeneous coupling strengths g i -the shape can vary.
One way to characterise the variation of input-output curve shape is to consider the relationship between slope ratio β and the 'width' of the transition.As discussed in the main text, we define β by the ratio of slopes below and above  threshold.We define transition width by considering the curvature of the input-output curve; the curvature should vanish in both the normal and lasing phases (where output is proportional to input).As such, we can use the peak curvature to numerically locate the critical power P c , and the half-width half-maximum of curvature as a transition width ∆P .We use the right half-width to avoid problems when the transition width becomes comparable to the critical power.With these defined, we can calculate a dimensionless parameter P c /(∆P ) which is a characteristic of the lasing transition.Figure S7 shows the relation between the dimensionless parameters β and P c /(∆P ).[5].The points correspond to calculations with disordered couplings gi, choosing gi drawn from a distribution with standard deviations as indicated in the legend.Notably, with this disorder, there is no longer a one-to-one relation between transition width and slope ratio: i.e. the shape of the input-output curve varies.Parameters used (matching Fig. 2(i,j)): Γz = Γ ↓ , κ = 1.74Γ ↓ .In the homogeneous case, a range of values of g are considered to map out the range of β.In the inhomogeneous cases, we fix the root mean square of gi to be 1.5Γ ↓ , and sample distributions with widths as indicated in the legend.
If the input-output curve were of fixed shape, there would exist a one-to-one relationship between β and P c /(∆P ).We see first that the relationship predicted by our model without disorder differs from that of Ref. [5].We note further (not shown) that this relationship is now parameter dependent, and would change depending on, e.g.Γ z /Γ ↓ , κ/Γ ↓ .More importantly, the one-to-one relation clearly breaks down when we consider inhomogeneous couplings, with varying g i .This demonstrates that the input-output curve does not in general have a fixed shape.

FIG. 2 .
FIG. 2. Emission dependence on excitation power and number of emitters.a, Absorption (green) and emission (orange) spectra of methylene blue in solution, and typical NPoM darkfield scattering spectrum (grey).b, Dashed lines: Example emission spectra for different expected number of emitters N in NPoM.Solid curves are averaged over > 50 NPoMs.c, Spectral integrated intensity vs N (blue circles) and linear fit (solid curve).Error bar is the standard error.d-g, Spectral integrated intensity vs power density for different N , colours correspond to different NPoMs.Dashed coloured lines are theoretical fits (see Methods) and the dashed grey lines are guide-to-the eye linear trends, which clearly deviate from the N > 2 data.h, Extracted slope ratio β vs N with predicted theoretical trend (solid curve), bars give range.i,j, Theoretical inputoutput curves, and extracted distribution of β, showing effect of disorder in couplings gi.Black line shows the homogeneous limit, gi = g = 1.5Γ ↓ for N = 10 emitters, with Γz = Γ ↓ , κ = 1.74Γ ↓ .The red and pink curves (and corresponding distributions of β) arise by sampling gi from a Gaussian distribution with root mean square g = 1.5Γ ↓ , and standard deviation g/3, g respectively.The thick black line corresponds to the homogeneous case, where no disorder is present.

FIG. 3 .
FIG. 3. Emission spectra and coherence of nanolaser.a, Emission spectra for increasing excitation powers P , with dark-field (DF) and solution PL. b, Schematic of Michelson interferometer to measure spectral coherence.Emission from NPoM is split in two by 50:50 beamsplitter (BS), reflected from mirrors (M1,M2) with delay in one arm and recombined, filtered (F) and sent to spectrometer (SP).c, Spectral fringes with average powers below (upper panel) and above (lower panel) lasing threshold (P th ).d, Visibility vs wavelength at different delays τ .e, Far-field Fourier space images of emission below and above threshold, divided into annular ring and centre.The integrated intensities in these regions, Iring and Icentre allow us to define an intensity ratio r = Iring/Icentre both below threshold (r b ) and above threshold (ra).Numerical aperture 0.9 limits collection angle < 64 • .f, NPoM radiative modes S10 and S11 that emit at high and low angles, respectively.g, Histogram of relative intensities in the high-angle ring ra/r b .
FIG. S2. (a-c) Input-output curves of several NPoM samples (> 50, open red circles) each for different number of molecules N .The filled red circles are the averages over all the NPoM samples

FIG. S4 .
FIG. S4.Comparison of surface-enhance Raman scattering using a continuous wave (CW) excitation (wavelength = 633 nm) and emission from NPoM at different pulsed laser powers.
FIG. S5.Reduced χ 2 maps obtained by optimizing the values of Γ ↓ , κ for a fixed value of Γz = 10 for N=5 (a) and N=10 (b).The green cross shows the minimum value of the χ 2R and the grey curves are the contour lines of the 1σ and 2σ confidence intervals.As discussed in the text, the units of these parameters are arbitrary, and only (dimensionless) ratios matter for the fitting.

FIG
FIG. S6.(a-c)The extracted coupling strength g vs intensity scaling, showing uncorrelation between the parameters as indicated by the guide-to-the-eye blue line.(d-f) Intensity scaling vs emission below threshold and the orange curve is a linear fit.
FIG. S7.Variation of shape of input-output curve.(a) Definition of transition width, extracted from the ratio of critical power Pc as defined by peak curvature, to right half-width half-maximum of curvature ∆P .Here curvature is the second derivative of input-output curve.(b) Relation between transition width, Pc/(∆P ), and slope-ratio β.The thick black line indicates the homogeneous limit of our model (no disorder in couplings.The gray line instead shows the weak-coupling formula of Ref.[5].The points correspond to calculations with disordered couplings gi, choosing gi drawn from a distribution with standard deviations as indicated in the legend.Notably, with this disorder, there is no longer a one-to-one relation between transition width and slope ratio: i.e. the shape of the input-output curve varies.Parameters used (matching Fig.2(i,j)): Γz = Γ ↓ , κ = 1.74Γ ↓ .In the homogeneous case, a range of values of g are considered to map out the range of β.In the inhomogeneous cases, we fix the root mean square of gi to be 1.5Γ ↓ , and sample distributions with widths as indicated in the legend.