Analysis of a multi-mode plasmonic nano-laser with an inhomogeneous distribution of molecular emitters

TL;DR

This paper extends laser theory to analyze a plasmonic nano-laser with inhomogeneous molecular distribution, revealing how molecular arrangement and driving strength influence lasing and plasmon excitation limits.

Contribution

It introduces a model accounting for molecular inhomogeneity and mode-correlation in a plasmonic nano-laser, providing insights into lasing conditions and excitation limits.

Findings

Large number of strongly driven molecules is needed for lasing.

Molecular inhomogeneity affects plasmon mode excitation.

Compact molecular arrangements can limit plasmon excitation.

Abstract

We extend Lamb's reduced density matrix laser theory to analyze the inhomogeneous molecular couplings and the mode-correlation in a plasmonic nano-laser consisting of a gold sphere and many dye molecules interacting with a driving optical field and with the quantized plasmon modes. The molecular inhomogeneity is accounted for by simulating their random distribution around the sphere. Our analysis shows that in order to obtain lasing we must employ a large number of strongly driven molecules to compensate strong damping of the plasmon modes. The compact molecular arrangement, however, can lead to molecular energy-shifts and thus reduce the excitation of the plasmon modes and ultimately suggests a maximum limit for the plasmon excitation for any specific system.

Tables

Table 1. Table 1: Physical parameters (for explanation see text)

$\hslash {\omega }_{\mathrm{pl}}$	$2.6$ eV	$\hslash {\omega }_{eg}$	$2.6$ eV
$\hslash {\gamma }_{\mathrm{pl}}$	$100$ meV	$\hslash {\omega }_{fg}$	$2.7$ eV
$d_{\mathrm{pl}}$	$2925$ D	$d_{gf}$	$16$ D
$E_0$	$1.2\times {10}^8$ V/m	$d_{ge}$	$14.4$ D
$\hslash {\omega }_0$	$2.7$ eV	$\hslash k_{f\to e}$	$100$ meV
		others	$0$ meV

Equations

$$\frac{\partial }{\partial t}\hat{\rho }=-\frac{i}{\hslash }\left[H_{\mathrm{pl}}+H_{\mathrm{e}}+V_{\mathrm{pl}-\mathrm{e}}+V_{\mathrm{e}}\left(t\right),\hat{\rho }\right]-\mathcal{D}\left[\hat{\rho }\right],$$

$$\mathcal{D}\left[\hat{\rho }\right]=\left(1/2\right)\sum _u{k}_u\left({\left[{\hat{L}}_u^{+}{\hat{L}}_u,\hat{\rho }\right]}_{+}-2{\hat{L}}_u\hat{\rho }{\hat{L}}_u^{+}\right).$$

$$\frac{\partial }{\partial t}{P}_{\mu }=-\sum _{j=1}^3\left[\left({\gamma }_j{\mu }_j+\kappa {}_{\mu }^{\left(j\right)}\right)P_{\mu }-\sum _{k=1}^3{\eta }_{\mu }^{\left(jk\right)}{P}_{{\mu }_k-1}\right]$$

$$+\sum _{j=1}^3\left[\left({\gamma }_j\left({\mu }_j+1\right)+{\kappa }_{{\mu }_j+1}^{\left(j\right)}\right)P_{{\mu }_j+1}-\sum _{k=1}^3{\eta }_{{\mu }_j+1}^{\left(jk\right)}{P}_{{\mu }_j+1{\mu }_k-1}\right],$$

$$P_{\mu }=\frac{{\sum }_{k=1}^3\left({\sum }_{j=1}^3{\eta }_{\mu }^{\left(jk\right)}\right)P_{{\mu }_k-1}}{{\sum }_{j=1}^3\left({\gamma }_j{\mu }_j+\kappa {}_{\mu }^{\left(j\right)}\right)}.$$

$$\frac{\partial }{\partial t}{\rho }_{\mu \nu }=-\sum _{j=1}^3\left(i{\omega }_j\left({\mu }_j-{\nu }_j\right)+{\gamma }_j\left[\left({\mu }_j+{\nu }_j\right)/2\right]\right){\rho }_{\mu \nu }$$

$$+\sum _{j=1}^3{\gamma }_j\sqrt{\left({\mu }_j+1\right)\left({\nu }_j+1\right)}{\rho }_{{\mu }_j+1{\nu }_j+1}$$

$$\displaystyle-i\sum_{n=1}^{N_{\mathrm{e}}}\sum_{j=1}^{3}v_{ge}^{\left(jn\right)}\Big{[}\sqrt{\mu_{j}}\rho_{e\mu_{j}-1,g\nu}^{\left(n\right)}-\sqrt{\nu_{j}}\rho_{g\mu,e\nu_{j}-1}^{\left(n\right)}$$

$$\displaystyle-\left(\sqrt{\nu_{j}+1}\rho_{e\mu,g\nu_{j}+1}^{\left(n\right)}-\sqrt{\mu_{j}+1}\rho_{g\mu_{j}+1,e\nu}^{\left(n\right)}\right)\Big{]}.$$

$$\displaystyle-\sum_{j}\Big{(}\left[i\omega_{\mu\nu}-\gamma_{j}\left(\mu_{j}+\nu_{j}\right)/2\right]\rho_{a\mu,b\nu}^{\left(n\right)}$$

$$\displaystyle+\gamma_{j}\sqrt{\left(\mu_{j}+1\right)\left(\nu_{j}+1\right)}\rho_{a\mu_{j}+1,b\nu_{j}+1}^{\left(n\right)}\Big{)}\to-i\tilde{\omega}_{\mu\nu}\rho_{a\mu,b\nu}^{\left(n\right)}.$$

$$\tilde{\omega}_{\mu\nu}=\sum_{j}\Big{(}\omega_{j}\left(\mu_{j}-\nu_{j}\right)-i\gamma_{j}\left[\left(\mu_{j}+\nu_{j}\right)/2-\sqrt{\mu_{j}\nu_{j}}\right]\Big{)}.$$

$$\frac{\partial }{\partial t}{\rho }_{g\mu ,g\nu }^{\left(n\right)}=-i{\tilde{\omega }}_{\mu \nu }{\rho }_{g\mu ,g\nu }^{\left(n\right)}-\left(k_{g\to f}^{\left(n\right)}+k_{g\to e}^{\left(n\right)}\right){\rho }_{g\mu ,g\nu }^{\left(n\right)}$$

$$+k_{f\to g}^{\left(n\right)}{\rho }_{f\mu ,f\nu }^{\left(n\right)}+k_{e\to g}^{\left(n\right)}{\rho }_{e\mu ,e\nu }^{\left(n\right)}-i{v}_{gf}^{\left(n\right)}\left({\tilde{\rho }}_{f\mu ,g\nu }^{\left(n\right)}-{\tilde{\rho }}_{g\mu ,f\nu }^{\left(n\right)}\right)$$

$$+i\sum _{j=1}^3{v}_{ge}^{\left(jn\right)}\left(\sqrt{{\nu }_j}{\rho }_{g\mu ,e{\nu }_j-1}^{\left(n\right)}-\sqrt{{\mu }_j}{\rho }_{e{\mu }_j-1,g\nu }^{\left(n\right)}\right),$$

$$\frac{\partial }{\partial t}{\rho }_{f\mu ,f\nu }^{\left(n\right)}=-i{\tilde{\omega }}_{\mu \nu }{\rho }_{f\mu ,f\nu }^{\left(n\right)}-\left(k_{f\to g}^{\left(n\right)}+k_{f\to e}^{\left(n\right)}\right){\rho }_{f\mu ,f\nu }^{\left(n\right)}$$

$$+k_{g\to f}^{\left(n\right)}{\rho }_{g\mu ,g\nu }^{\left(n\right)}+k_{e\to f}^{\left(n\right)}{\rho }_{e\mu ,e\nu }^{\left(n\right)}+i{v}_{gf}^{\left(n\right)}\left({\tilde{\rho }}_{f\mu ,g\nu }^{\left(n\right)}-{\tilde{\rho }}_{g\mu ,f\nu }^{\left(n\right)}\right),$$

$$\frac{\partial }{\partial t}{\rho }_{e{\mu }_j-1,e{\nu }_j-1}^{\left(n\right)}=-i{\tilde{\omega }}_{{\mu }_j-1{\nu }_j-1}{\rho }_{e{\mu }_j-1,e{\nu }_j-1}^{\left(n\right)}+k_{g\to e}^{\left(n\right)}{\rho }_{g{\mu }_j-1,g{\nu }_j-1}^{\left(n\right)}$$

$$-\left(k_{e\to g}^{\left(n\right)}+k_{e\to f}^{\left(n\right)}\right){\rho }_{e{\mu }_j-1,e{\nu }_j-1}^{\left(n\right)}+k_{f\to e}^{\left(n\right)}{\rho }_{f{\mu }_j-1,f{\nu }_j-1}^{\left(n\right)}$$

$$+i{v}_{ge}^{\left(jn\right)}\left(\sqrt{{\nu }_j}{\rho }_{e{\mu }_j-1,g\nu }^{\left(n\right)}-\sqrt{{\mu }_j}{\rho }_{g\mu ,e{\nu }_j-1}^{\left(n\right)}\right).$$

$$\frac{\partial }{\partial t}{\tilde{\rho }}_{g\mu ,f\nu }^{\left(n\right)}=-i\left({\tilde{\omega }}_{gf}^{\left(n\right)}+{\tilde{\omega }}_{\mu \nu }\right){\tilde{\rho }}_{g\mu ,f\nu }^{\left(n\right)}$$

$$-i\sum _{j=1}^3{v}_{ge}^{\left(jn\right)}\sqrt{{\mu }_j}{\tilde{\rho }}_{e{\mu }_j-1,f\nu }^{\left(n\right)}-i{v}_{gf}^{\left(n\right)}\left({\rho }_{f\mu ,f\nu }^{\left(n\right)}-{\rho }_{g\mu ,g\nu }^{\left(n\right)}\right),$$

$$\frac{\partial }{\partial t}{\tilde{\rho }}_{f\mu ,g\nu }^{\left(n\right)}=i\left({\tilde{\omega }}_{gf}^{\left(n\right)\ast }-{\tilde{\omega }}_{\mu \nu }\right){\tilde{\rho }}_{f\mu ,g\nu }^{\left(n\right)}$$

$$+i\sum _{j=1}^3{v}_{ge}^{\left(jn\right)}\sqrt{v_j}{\tilde{\rho }}_{f\mu ,e{\nu }_j-1}^{\left(n\right)}+i{v}_{gf}^{\left(n\right)}\left({\rho }_{f\mu ,f\nu }^{\left(n\right)}-{\rho }_{g\mu ,g\nu }^{\left(n\right)}\right),$$

$$\frac{\partial }{\partial t}{\rho }_{g\mu ,e{\nu }_j-1}^{\left(n\right)}=i\left({\tilde{\omega }}_{eg}^{\left(n\right)\ast }-{\tilde{\omega }}_{\mu {\nu }_j-1}\right){\rho }_{g\mu ,e{\nu }_j-1}^{\left(n\right)}-i{v}_{gf}^{\left(n\right)}{\rho }_{f\mu ,e{\nu }_j-1}^{\left(n\right)}$$

$$+i{v}_{ge}^{\left(jn\right)}\left(\sqrt{{\nu }_j}{\rho }_{g\mu ,g\nu }-\sqrt{{\mu }_j}{\rho }_{e{\mu }_j-1,e{\nu }_j-1}^{\left(n\right)}\right),$$

$$\frac{\partial }{\partial t}{\rho }_{e{\mu }_j-1,g\nu }^{\left(n\right)}=-i\left({\tilde{\omega }}_{eg}^{\left(n\right)}+{\tilde{\omega }}_{{\mu }_j-1\nu }\right){\rho }_{e{\mu }_j-1,g\nu }^{\left(n\right)}+i{v}_{gf}^{\left(n\right)}{\rho }_{e{\mu }_j-1,f\nu }^{\left(n\right)}$$

$$+i{v}_{ge}^{\left(jn\right)}\left(\sqrt{{\nu }_j}{\rho }_{e{\mu }_j-1,e{\nu }_j-1}^{\left(n\right)}-\sqrt{{\mu }_j}{\rho }_{g\mu ,g\nu }\right),$$

$$\frac{\partial }{\partial t}{\tilde{\rho }}_{e{\mu }_j-1,f\nu }^{\left(n\right)}=-i\left({\tilde{\omega }}_{ef}^{\left(n\right)}+{\tilde{\omega }}_{{\mu }_j-1\nu }\right){\tilde{\rho }}_{e{\mu }_j-1,f\nu }^{\left(n\right)}$$

$$+i{v}_{gf}^{\left(n\right)}{\rho }_{e{\mu }_j-1,g\nu }^{\left(n\right)}-i{v}_{ge}^{\left(jn\right)}\sqrt{{\mu }_j}{\rho }_{g\mu ,f\nu }^{\left(n\right)},$$

$$\frac{\partial }{\partial t}{\rho }_{f\mu ,e{\nu }_j-1}^{\left(n\right)}=i\left({\tilde{\omega }}_{ef}^{\left(n\right)\ast }-{\tilde{\omega }}_{\mu {\nu }_j-1}\right){\rho }_{f\mu ,e{\nu }_j-1}^{\left(n\right)}$$

$$-i{v}_{gf}^{\left(n\right)}{\tilde{\rho }}_{g\mu ,e{\nu }_j-1}^{\left(n\right)}+i{v}_{ge}^{\left(jn\right)}\sqrt{{\nu }_j}{\tilde{\rho }}_{f\mu ,g\nu }^{\left(n\right)}.$$

$$\left({\tilde{\omega }}_{eg}^{\left(n\right)\ast }-{\tilde{\omega }}_{\mu {\nu }_j-1}\right){\rho }_{g\mu ,e{\nu }_j-1}^{\left(n\right)}=v_{gf}^{\left(n\right)}{\rho }_{f\mu ,e{\nu }_j-1}^{\left(n\right)}$$

$$-v_{ge}^{\left(jn\right)}\left(\sqrt{{\nu }_j}{\rho }_{g\mu ,g\nu }-\sqrt{{\mu }_j}{\rho }_{e{\mu }_j-1,e{\nu }_j-1}^{\left(n\right)}\right),$$

$$\left({\tilde{\omega }}_{eg}^{\left(n\right)}+{\tilde{\omega }}_{{\mu }_j-1\nu }\right){\rho }_{e{\mu }_j-1,g\nu }^{\left(n\right)}=v_{gf}^{\left(n\right)}{\rho }_{e{\mu }_j-1,f\nu }^{\left(n\right)}$$

$$-v_{ge}^{\left(jn\right)}\left(\sqrt{{\mu }_j}{\rho }_{g\mu ,g\nu }-\sqrt{{\nu }_j}{\rho }_{e{\mu }_j-1,e{\nu }_j-1}^{\left(n\right)}\right),$$

$$\left({\tilde{\omega }}_{ef}^{\left(n\right)}+{\tilde{\omega }}_{{\mu }_j-1\nu }\right){\tilde{\rho }}_{e{\mu }_j-1,f\nu }^{\left(n\right)}=v_{gf}^{\left(n\right)}{\rho }_{e{\mu }_j-1,g\nu }^{\left(n\right)}-v_{ge}^{\left(jn\right)}\sqrt{{\mu }_j}{\rho }_{g\mu ,f\nu }^{\left(n\right)},$$

$$\left({\tilde{\omega }}_{ef}^{\left(n\right)\ast }-{\tilde{\omega }}_{\mu {\nu }_j-1}\right){\rho }_{f\mu ,e{\nu }_j-1}^{\left(n\right)}=v_{gf}^{\left(n\right)}{\tilde{\rho }}_{g\mu ,e{\nu }_j-1}^{\left(n\right)}-v_{ge}^{\left(jn\right)}\sqrt{{\nu }_j}{\tilde{\rho }}_{f\mu ,g\nu }^{\left(n\right)}.$$

$${\tilde{\rho }}_{e{\mu }_j-1,f\nu }^{\left(n\right)}=-{\Xi }_{\mu \nu }^{\left(jn\right)}{v}_{ge}^{\left(jn\right)}\sqrt{{\mu }_j}{\rho }_{g\mu ,f\nu }^{\left(n\right)}$$

Full text

I Introduction

The interaction between metals and light has been investigated for more than a century with Maxwell’s electromagnetic theory. One essential insight obtained is that the electromagnetic (EM) field is enhanced and localized around metal nano-particles (MNP) and on the interfaces between metallic films and dielectrics MPelton due to the excitation of surface plasmons involving collective oscillations of conductance electrons in the metal. The enhancement boosts the interaction between quantum emitters and the EM field MSTame ; RMMa-0 ; MPelton ; PBerini and thus leads to enhanced absorption NICade ; YZelinskyy , emission PAnger ; YZhang-4 and Raman scattering MFleischmann ; PJohansson . This can be utilized to improve the sensibility of spectroscopic instruments SYDing and the efficiency of LEDs XFGu ; NGao and solar cells HAAtwater ; LJWu .

The localization introduces EM modes with mode volumes that are not limited by the wavelength of free-space light YYin ; PBerini . These modes can be excited if externally pumped quantum emitters are placed near MNPs or metallic films. Under suitable conditions, the energy loss of those modes can be even compensated and the system can achieve lasing. This phenomenon known as SPASER, was proposed by Bergman and Stockman Bergman and verified firstly by Noginov, et. al. MANoginov with an experiment involving a gold nano-sphere and many dye molecules. Since then many experimental demonstrations have been reported with structures like semiconductor wires JHo ; RFOulton ; CYWu ; YJLu-1 ; YJLu ; YHou ; QZhang ; TPPHSi ; YHChou ; BTChou /squares RMMa ; RMMa-1 on metallic films, semiconductor pillars MTHill ; MPNezhad ; SHkwon ; JHLee ; MKha ; KDing /dots AMatsudira ; CYLu /wires CYLu-1 ; SWChang inside metallic cavities as well as dye molecules in periodically arranged MNP arrays JYSuh ; WZhou ; AYang ; AYang-1 ; AHSchokker .

In order to theoretically describe these systems, we have to determine the lasing modes and consider how the gain medium transfers energy to these modes. The modes can be analyzed by solving Maxwell’s equations analytically YHChou ; SWChang-1 or numerically JHo ; RFOulton ; CYWu ; YJLu-1 ; YJLu ; YHou ; QZhang ; TPPHSi ; YHChou ; BTChou ; MTHill ; MPNezhad ; SHkwon ; JHLee ; MKha ; KDing ; JYSuh ; WZhou ; AYang ; AYang-1 ; AHSchokker ; NLi . The energy transfer requires us to model the gain medium as a random spatial distribution of multi-level emitters. The multi-level model allows us to couple some levels with external driving fields or electron reservoirs to describe pumping mechanism and couple other levels with the lasing modes. It means that we have to deal with a complex theoretical problem involving many emitters, many levels and many modes. To reduce the complexity, semi-classical theories have been developed and utilized, for example, rate equations AMatsudira ; CYLu ; CYLu-1 and Maxwell-Bloch equations WZhou ; AYang ; AYang-1 . Because of the mean-field approximations involved, these theories, however, yield no statistical information about the the lasing modes and thus advanced full quantum laser theories are needed. Unfortunately, so far those quantum theories treat the emitters as identical two-level systems MRichter ; VMPar and are thus incapable of dealing with randomly distributed emitters.

Most existing theories can be viewed as effective descriptions, where the emitters are treated only in an average sense. They can reproduce main characteristics measured in experiments because the inhomogeneity of the emitters becomes irrelevant if huge amounts of them are involved. However, this may not be the case in the plasmonic nano-laser. Because of the strong inhomogeneous subwavelength distribution of the near-field in the nano-laser, the spatial distribution of emitters can significantly affect how they convert the pumping energy into the plasmon energy and even determine whether the systems can achieve lasing or not. By analyzing this influence, we can achieve more insights about the systems and more importantly understand how to improve the system performance by engineering the spatial distribution.

In this article, we provide a systematic analysis of the inhomogeneity of molecular emitters in a nano-laser of Fig. 1(a), which resembles the one studied in MANoginov . As the basis for our analysis, we will firstly describe our theoretical model in Sec. II. To account for multi-molecules and multi-modes, we have extended the density matrix laser theory of Lamb MSargent in our model. By following the procedure developed in YZhang ; YZhang-1 , we first establish a reduced density operator equation for the entire system and then derive a quantum master equation for the reduced density matrix (RDM) of the plasmon modes. Our extended theory allows us to analyze how the molecular distribution affects the plasmon statistics and the molecule-induced mode-correlations. This analysis is presented in Sec. III. In the end, we summarize our findings and present an outlook for future work.

II Theoretical Model

As indicated by Fig.1 (a), we consider a random arrangement of molecular emitters separated by more than $2.5$ nm from the surface of a gold nano-sphere of 10 nm radius. The separation guarantees that electron tunneling is suppressed KJSavage and the molecule-MNP coupling is dominated by Coulomb coupling. The molecules are assumed to be resonant with the dipole plasmons of the sphere. The higher multipole plasmons have minor influence on the system YZhang-3 and contribute only as an off-resonant reservoir to the excited-state decay rate of the molecules JGersten .

II.1 Reduced Density Operator Equation

For a nano-sphere, there are three degenerate dipole plasmon modes with transition dipole moments pointing along three orthogonal axes. Therefore, we can label them by $j=x,y,z$ or $j=1,2,3$. These modes can be described as quantum harmonic oscillators with Hamiltonian $H_{\mathrm{pl}}=\sum_{j}\hbar\omega_{j}C_{j}^{+}C_{j}$, where $C_{j}^{+}$ and $C_{j}$ are creation and annihilation operators and $\hbar\omega_{j}=\hbar\omega_{\mathrm{pl}}$ is their excitation energy Gweick . We describe the molecules as three-level systems with the internal energy level scheme and transitions shown in Fig.1(b). The molecular Hamiltonian reads $H_{{\rm e}}=\sum_{n=1}^{N_{{\rm e}}}\sum_{a_{n}}E_{na}\left|a_{n}\right\rangle\left\langle a_{n}\right|$ where ground states $\left|a_{n}=g_{n}\right\rangle$, first $\left|a_{n}=e_{n}\right\rangle$ and second $\left|a_{n}=f_{n}\right\rangle$ excited states have the energies $E_{ng}$,$E_{ne}$,$E_{nf}$, respectively ThreeLevels . We assume that the plasmon modes are resonant with the ground-to-first excited state transition, cf. the blue arrow in Fig.1(b), and introduces the coupling Hamiltonian $V_{{\rm pl-e}}=\hbar\sum_{n=1}^{N_{{\rm e}}}v_{ge}^{\left(jn\right)}\left(C_{j}^{+}\left|g_{n}\right\rangle\left\langle e_{n}\right|+{\rm h.c.}\right)$ in the rotating wave approximation. Here, the coefficient $\hbar v_{ge}^{\left(jn\right)}=\left[\mathbf{d}_{ge}^{\left(n\right)}\cdot\mathbf{d}_{j}-3\left(\mathbf{d}_{ge}^{\left(n\right)}\cdot\hat{\mathbf{x}}_{n}\right)\left(\mathbf{d}_{j}\cdot\hat{\mathbf{x}}_{n}\right)\right]/\left|\mathbf{X}_{n}\right|^{3}$ is determined by the transition dipole moment $\mathbf{d}_{ge}^{\left(n\right)}$ of the molecules, $\mathrm{\mathbf{d}}_{j}=d_{\mathrm{pl}}\mathrm{\mathbf{e}}_{j}$ of the plasmon modes as well as the distances $X_{n}$ and directional unit vectors $\hat{\mathbf{x}}_{n}$ connecting the $n^{\rm th}$ molecule and the sphere-center. We assume that a classical driving field is resonant with the ground-to-second excited state transition, cf. the red arrow in Fig.1(b), and introduces the coupling Hamiltonian $V_{{\rm e}}\left(t\right)=\hbar\sum_{n=1}^{N_{{\rm e}}}v_{gf}^{\left(n\right)}\left(e^{i\omega_{0}t}\left|g_{n}\right\rangle\left\langle f_{n}\right|+{\rm h.c.}\right)$ in the rotating wave approximation. Here, the coefficient $\hbar v_{gf}^{\left(n\right)}=\mathbf{d}_{gf}^{\left(n\right)}\cdot\mathbf{n}E_{0}$ is determined by another molecular transition dipole moment $\mathbf{d}_{gf}^{\left(n\right)}$ and the driving field is specified by a frequency $\omega_{0}$, a polarization vector $\mathbf{n}$ and an amplitude $E_{0}$ NFO . Here, we consider continuous optical excitation and thus $E_{0}$ is time-independent.

The density operator $\hat{\rho}$ for the quantized plasmon modes and the molecular emitters obeys the following quantum master equation

[TABLE]

where the system dissipation is accounted for by the Lindblad terms:

[TABLE]

The damping of the plasmon modes is included by terms with $k_{u}=\gamma_{j}=\gamma_{\mathrm{pl}}$, $\hat{L}_{u}=C_{j}$ for each mode $j$. The decay processes of the molecules are included by terms with $k_{u}=k_{a\to b}^{\left(n\right)}$, $\hat{L}_{u}=\left|b_{n}\right\rangle\left\langle a_{n}\right|$ for $E_{na}>E_{nb}$ for each molecule, cf. the black arrows in Fig.1(b). For the sake of simplicity, we ignore pure molecular dephasing.

II.2 Plasmon Reduced Density Matrix Equation

The main goal of our analysis is to determine the plasmon state populations (probabilities) and correlations as quantified by the reduced density matrix (RDM) with elements $\rho_{\mu\nu}\equiv\text{tr}_{\text{S}}\left\{\hat{\rho}\left(t\right)\left|\nu\right\rangle\left\langle\mu\right|\right\}$, where $\text{tr}_{\text{S}}$ denotes the trace over the system and $\left|\mu\right\rangle\equiv\left|\mu_{x}\right\rangle\left|\mu_{y}\right\rangle\left|\mu_{z}\right\rangle$ and $\left|\nu\right\rangle\equiv\left|\nu_{x}\right\rangle\left|\nu_{y}\right\rangle\left|\nu_{z}\right\rangle$ denote product states of the plasmon occupation number Fock states. From Eq. (1), we observe that $\rho_{\mu\nu}$ depends on the molecule-plasmon correlations $\rho_{e\mu_{j}-1,g\nu}^{\left(n\right)}$, $\rho_{g\mu,e\nu_{j}-1}^{\left(n\right)}$ and $\rho_{e\mu,g\nu_{j}+1}^{\left(n\right)}$, $\rho_{g\mu_{j}+1,e\nu}^{\left(n\right)}$ with a general definition $\rho_{a\mu,b\nu}^{\left(n\right)}\equiv\text{tr}_{\text{S}}\left\{\hat{\rho}\left(t\right)\left|b_{n}\right\rangle\left\langle a_{n}\right|\times\left|\nu\right\rangle\left\langle\mu\right|\right\}$, cf. Appendix B.

The equations for the correlations also follow from Eq. (1). These equations result in dependence between the plasmon RDM and the correlations, shown in Fig. 1(c), which is caused by the couplings and the dissipation rates in the master equation (1). This dependence also indicates our procedure to solve those inter-dependent equations: Because both molecular and plasmonic dissipation rates contribute to the decay of the correlations, they must decay faster and thus may adiabatically MOScully follow the plasmon RDM elements which are only affected by the plasmon damping. Because of the molecular dissipation, the correlations represented within the blue dashed box of Fig. 1(c) depend on the correlations outside the box. Fortunately, they all can be expressed as functions of the plasmon RDM because of the symmetry hidden in the coupling Hamiltonians. Finally, we back substitute these expressions and obtain closed dynamic equations for the plasmon RDM, where the molecules contribute by several coefficients, cf. Eq.(87) in Appendix B.

The diagonal elements $P_{\mu}\equiv\rho_{\mu\mu}$ are the populations (the probabilities) of the plasmon number states $\left|\mu\right\rangle$ while the off-diagonal elements $\rho_{\mu\nu}$ ($\mu\neq\nu$) represent the coherence of the plasmons. Here, we focus on the populations by solving the equations for those diagonal elements:

[TABLE]

where the rates $\kappa_{\mu}^{\left(j\right)}\equiv-\sum_{n=1}^{N_{\mathrm{e}}}\alpha{}_{\mu\mu}^{\left(jn\right)}$ and $\eta_{\mu}^{\left(jk\right)}\equiv-\sum_{n=1}^{N_{\mathrm{e}}}\beta_{\mu\mu}^{\left(jkn\right)}$ include contributions from individual molecule $\alpha{}_{\mu\mu}^{\left(jn\right)}$ and $\beta_{\mu\mu}^{\left(jkn\right)}$ respectively, cf. Eqs. (88) and (89) in the Appendix B. Here, $\mu_{j}\pm 1$ for $j=x$ indicates $(\mu_{x}\pm 1,\mu_{y},\mu_{z})$ and $\mu_{j}+1\mu_{k}-1$ for $j=x$ and $k=y$ denotes $(\mu_{x}+1,\mu_{y}-1,\mu_{z})$. Since the former rates decrease the population of higher plasmon states but increase that of lower states, they can be interpreted as molecule-induced plasmon damping rates. Since the latter rates have the opposite effect on the population, they can be interpreted as molecule-induced plasmon pumping rates. The latter rates depend on two plasmon mode indices and thus account for correlation between different plasmon modes induced by the molecules. These rates can be considered as extended Einstein’s AB coefficients accounting for the multi-plasmon modes, the molecular pumping mechanism and the molecular inhomogeneity.

At steady-state the second line of Eq. (3) is recovered if we replace $\mu_{j}$ by $\mu_{j}+1$ on the right side of the first line, which suggests a recursion relation of the populations. We obtain such a relation by setting the first line to zero:

[TABLE]

Together with the normalization condition $\sum_{\mu}P_{\mu}=1$, the above relation can be utilized to easily calculate the populations according to the procedure outlined in Fig. 8 in Appendix B. Although $P_{\mu_{x}\mu_{y}\mu_{z}}$ contains all the information about the three dipole plasmons, it is more intuitive to consider physical quantities related to one or two dipole plasmon. We calculate them by tracing out one mode to get $P_{\mu_{x}\mu_{y}}=\sum_{\mu_{z}}P_{\mu_{x}\mu_{y}\mu_{z}}$, $P_{\mu_{y}\mu_{z}}=\sum_{\mu_{x}}P_{\mu_{x}\mu_{y}\mu_{z}}$ and $P_{\mu_{x}\mu_{z}}=\sum_{\mu_{y}}P_{\mu_{x}\mu_{y}\mu_{z}}$ (the joint population of two modes) or by tracing out two modes to get $P_{\mu_{x}}=\sum_{\mu_{y}\mu_{z}}P_{\mu_{x}\mu_{y}\mu_{z}}$, $P_{\mu_{y}}=\sum_{\mu_{x}\mu_{z}}P_{\mu_{x}\mu_{y}\mu_{z}}$ and $P_{\mu_{z}}=\sum_{\mu_{x}\mu_{y}}P_{\mu_{x}\mu_{y}\mu_{z}}$ (the reduced population of one mode). We can also quantify the strength of plasmon excitation with the so-called plasmon mean numbers: $N_{j}\equiv\sum_{\mu_{j}}\mu_{j}P_{\mu_{j}}$ and the plasmon statistics with the so-called (steady-state) second order correlation functions: $g_{j}^{\left(2\right)}(0)\equiv\sum_{\mu_{j}}\mu_{j}\left(\mu_{j}-1\right)P_{\mu_{j}}$. To analyze how the individual molecule contributes to the plasmon excitation, we can calculate the population of individual molecular states: $P_{g}^{\left(n\right)}\equiv\sum_{\mu}\rho_{g\mu,g\mu}^{\left(n\right)}$, $P_{f}^{\left(n\right)}\equiv\sum_{\mu}\rho_{f\mu,f\mu}^{\left(n\right)}$ and $P_{e}^{\left(n\right)}\equiv\sum_{\mu}\rho_{e\mu,e\mu}^{\left(n\right)}$, which are actually determined by $P_{\mu_{x}\mu_{y}\mu_{z}}$, cf. Eqs. (123), (124) and (125) in Appendix C.

III Results

The above theoretical model provides clues about how the molecular inhomogeneity may affect the system performance. The molecular inhomogeneity mainly originates from their positions and orientations of their transition dipole moments, which leads to that all the molecules interact with the three modes simultaneously but with random strengths. Since this situation is too complex, we shift its discussion to the end and first consider a special configuration where all the molecules are located along the equator of the gold sphere, cf. Fig.2 (a).

III.1 Configuration with Single Dipole Plasmon Mode

For the configuration in Fig. 2 (a), the molecule-plasmon coupling is reduced to $\hbar v_{ge}^{\left(jn\right)}=\delta_{j,z}\left(\pm d_{ge}^{(n)}d_{\rm pl}/X_{n}^{3}\right)$ with positive (negative) sign for the molecules oriented along the positive (negative) z-axis. Obviously, the dipole plasmon x- and y-mode are not involved and thus can be ignored in the following analysis. The coupling depends inversely on the cubic of the molecule-sphere center distance $X_{n}=a_{\rm MNP}+d_{n}$, cf. the black squares in Fig. 2 (b). Here, $a_{\rm MNP}$ is the radius of the sphere and $d_{n}$ the distance to the sphere-surface. In contrast, the driving field coupling $\hbar v_{gf}^{\left(n\right)}$ on the molecules depends only on the molecular orientations but not the positions, cf. the red triangles in Fig. 2 (b). If all the molecules have the same distance to the sphere-surface, they are equivalent and the resulting ideal system has been already analyzed in YZhang-1 . There, we focused on the pumping mechanism and found the optimal parameters of the system leading to the strongest plasmon excitation, cf. Table 1 in Appendix A. These parameters will be used as reference parameters for the following simulations.

To compensate the strong plasmon damping, the number of molecules coupled strongly with the plasmons is an essential parameter. It was demonstrated in the experiment XGMeng that the system properties like emission wavelength, intensity and pumping threshold strongly depend on the concentration (number) of the molecules. Here, we analyze this dependence from three aspects: density of molecules, spatial extension of molecular layer and molecular level shift.

As indicated by Fig.2 (b), the molecules close to the sphere couple strongly with the plasmons. Therefore, those molecules contribute more to the plasmon excitation than other molecules. By increasing the molecular density, we increase the number of molecules and thus the plasmon excitation. This is clearly reflected in Fig. 2 (c) by the increased populations $P_{\mu_{z}}$ of higher plasmon excited states and the increased plasmon mean number $N_{z}$ (the black dots and curve in the inset) with increasing number of molecules $N_{\mathrm{e}}$ from $10$ to $240$. For larger $N_{\mathrm{e}}$, $P_{\mu_{z}}$ resemble Poisson distributions indicating the formation of a coherent state and $N_{z}$ approaches $25$, which is much larger than unity and thus indicates that the system is lasing. This conclusion is further confirmed by the $g_{z}^{\left(2\right)}(0)$-function, cf. the red dots and red curve in the inset of Fig. 2 (c), which approaches unity for large $N_{\mathrm{e}}$, i.e. the Poisson limit. The fluctuation of the dots around the curves in the inset is caused by different molecular distribution in each simulation and may thus represent fluctuations encountered in experiments.

To understand why the increasing molecular density can increase the plasmon excitation, here, we analyze the state population for every molecule $P_{a}^{\left(n\right)}$, cf. Fig. 2(d). First, we notice that the molecule-plasmon coupling $V_{{\rm pl-e}}$ leads to reversible processes (spontaneous emission, stimulated emission and absorption of the plasmons) since it enters into our master equation as a coherent coupling, cf. Eq. (1). These processes tend to balance the population of the molecular excited states $P_{e}^{\left(n\right)}$ and ground states $P_{g}^{\left(n\right)}$. This leads to the reduced $P_{e}^{\left(n\right)}$, cf. the red curves and arrow, and the increased $P_{g}^{\left(n\right)}$, cf. the black curves and arrow, with increasing $N_{\mathrm{e}}$. In addition, because the reduced coupling with increasing distance $d_{n}$ (cf. in Fig. 2(b)) reduces the rate of the processes, the $P_{e}^{\left(n\right)}$ increase and $P_{g}^{\left(n\right)}$ decrease with increasing $d_{n}$. The population of the higher excited state $P_{f}^{\left(n\right)}$ is mainly determined by the strong decay rate from this state to the middle excited state and thus is always smaller than the other populations.

In the following, we consider the effect of varying the spatial extent (width $W$) of the molecular layer, cf. Fig. 2(e), which can also be studied in experiments like MANoginov ; XGMeng by precisely controlling the synthesis time of the molecular layer. In this case, the molecules are added far away from the sphere surface and will thus contribute less to the plasmon excitation because of the reduced molecule-plasmon coupling, cf. Fig. 2 (b). As a result, the plasmon mean number $N_{z}$ and the $g_{z}^{\left(2\right)}\left(0\right)$-function saturate for large $W$ as displayed by Fig. 2(e). In addition, we find that the data points are close to the fitted curve for small $W$ but fluctuate a lot for large $W$. This can be easily understood with the change of the molecular state population $P_{a}^{\left(n\right)}$, cf. Fig. 2(f). When $W$ increases from 1 nm to $30$ nm, $P_{g}^{\left(n\right)},P_{e}^{\left(n\right)},P_{f}^{\left(n\right)}$ change dramatically, cf. the zoomed inset, since all molecules contribute to the plasmon excitation. Therefore, $N_{z}$ increases and the molecular inhomogeneity has less influence on the plasmon excitation. However, when $W$ increases further, $P_{g}^{\left(n\right)},P_{e}^{\left(n\right)},P_{f}^{\left(n\right)}$ change less and now the molecules distributed in the region near to the sphere will significantly affect $N_{z}$.

The compact molecular arrangement around the sphere implies that the molecules may directly interact with each other through Coulomb coupling. If the molecular concentration is very high, electron transfer between molecules becomes possible. Although this process may be relevant here, it is however beyond the scope of our theory. For not too high concentration, the electron transfer can be ignored but direct energy exchange between excited molecular dipoles can lead to energy-shift (inhomogeneous broadening). In principle, such effect can be accurately described by directly incorporating the inter-molecular energy exchange coupling into the system Hamiltonian in the master equation (1). However, here, we follow an easier, phenomenological way to account for such effect by introducing random energy shift $\delta E^{\left(n\right)}$ to individual molecule with a Gaussian distribution HHaken , cf. Fig. 3 (a), $p\left(\delta E\right)=\left(1/\sqrt{2\pi}\sigma\right)\exp\left\{-\left(\delta E\right)^{2}/\left(2\sigma^{2}\right)\right\}$ (with standard deviation $\sigma$). It means that the transition energies are modified as $\hbar\omega_{eg}^{\left(n\right)}=\hbar\omega_{eg}+\delta E^{\left(n\right)}$ and $\hbar\omega_{fg}^{\left(n\right)}=\hbar\omega_{fg}+\delta E^{\left(n\right)}$, compared to the values in Table I in Appendix A.

The consequence of energy-shift is to perturb the perfect resonant condition for the molecular pumping and the molecule-plasmon energy transfer assumed previously. This is reflected by the irregular change of the state populations $P_{a}^{\left(n\right)}$ for the molecules at similar distances to the sphere surface, cf. the dots in Fig. 3 (b). However, since the majority of molecules has no or small energy shift as shown in Fig. 3 (a), the populations $P_{a}^{\left(n\right)}$ in Fig. 3 (b) still roughly follow the same trend observed in Fig. 2 (d), cf. the solid lines. The broadening of the transition energies is also reflected in the shift of the plasmon state population $P_{\mu_{z}}$ to lower states, a reduced plasmon mean number $N_{z}$ and an increased $g_{z}^{(2)}(0)$-function with increasing deviation $\sigma$ of the energy-shift from 0 meV to 100 meV, cf. Fig. 3 (c). These features can be understood by analyzing the contribution of individual molecule through their state populations $P_{a}^{\left(n\right)}$. As shown in Fig. 3 (d), the populations of the molecular middle excited states $P_{e}^{\left(n\right)}$ increase while those of the ground states $P_{g}^{\left(n\right)}$ decrease with increasing $\sigma$. These results reflect that the molecules are on average less affected by the plasmons and thus contribute less to the plasmon excitation. The features described above indicate that the lasing performance is strongly affected by the inhomogeneous molecular energy-shift. Finally, we point out that the standard deviation $\sigma$ characterizes energy-shifts due to intra-molecular interactions and should hence depends on the molecular concentration.

III.2 Configuration with Two and Three Dipole Plasmon Modes

We now consider the more complex situation where the molecular transition dipole moments $\mathbf{d}_{ge}^{\left(n\right)}$ and $\mathbf{d}_{gf}^{\left(n\right)}$ orient randomly in the x-y plane, cf. Fig. 4(a). In this case, the molecules couple with the dipole plasmon $x-$ and $y-$mode simultaneously with random strength, cf. the upper panel of Fig. 4(b). In addition, the driving field coupling becomes also random as shown in the lower panel of Fig. 4(b) because it also depends on the orientations. This implies that the molecules at similar distance to the sphere-surface experience different couplings and this leads to the random population of molecular states $P_{a}^{\left(n\right)}$ as displayed in Fig. 4(c). However, because the maximum of the molecule-plasmon coupling decreases with increasing distance to the sphere-surface $d^{\left(n\right)}$, the distance-dependent averaged $P_{a}^{\left(n\right)}$ show a similar behavior as in Fig. 2(d). The co-existence of the x- and -y mode is directly illustrated by the joint plasmon state population $P_{\mu_{x}\mu_{y}}$ for a system with $N_{\mathrm{e}}=500$ molecules, cf. Fig. 4(d). Here, to better visualize $P_{\mu_{x}\mu_{y}}$, it is shown as a smooth surface. The population has a maximum around $\mu_{x}=25$ and $\mu_{y}=25$, which indicates that the both plasmon modes are excited to the same strength. In addition, we have also analyzed the plasmon state populations $P_{\mu_{x}},P_{\mu_{y}}$ , the plasmon mean numbers $N_{x},N_{y}$ and the $g_{x}^{\left(2\right)}\left(0\right)$- and $g_{y}^{\left(2\right)}\left(0\right)$-function for systems with increasing molecular density (number of molecules) in Fig. 6(a-c) and with increasing deviation $\sigma$ of molecular energy shift in Fig. 6(d-e) in Appendix A. Basically, they show similar features like those in Fig. 2(c) and Fig. 3 (c) respectively.

Finally, let us turn to the realistic configuration of Fig. 1(a). In this case, the randomly distributed molecules in three dimensions couple with the three plasmon modes in a similar pattern, cf. Fig. 5 (a) (see also the coupling of individual molecule in Fig. 7 (a) in Appendix A). This implies that all the plasmon modes will be excited by the molecules in similar way and this is reflected by the joint populations $P_{\mu_{x}\mu_{y}}$, $P_{\mu_{y}\mu_{z}}$ and $P_{\mu_{x}\mu_{z}}$ with a peak around $\left(10,10\right)$ for a system with $N_{\rm e}=800$ molecules, cf. Fig. 5 (b) (see also Fig. 7 (b,c) in Appendix A). In addition, we also find the increased population $P_{\mu_{x}}$, $P_{\mu_{y}}$, $P_{\mu_{z}}$ of higher plasmon states, the increased plasmon mean number $N_{x},N_{y},N_{z}$ as well as the reduced $g_{x}^{\left(2\right)}\left(0\right),g_{y}^{\left(2\right)}\left(0\right),g_{z}^{\left(2\right)}\left(0\right)$-functions with increasing number of molecules $N_{\mathrm{e}}$ ( cf. Fig. 7(d,e,f) in Appendix A).

In order to achieve same plasmon excitation per mode, we must double (triple) the number of molecules in the case with two (three) modes compared to the single mode case. incidentally, our results show that the polarization of the driving field alone does not cause significant asymmetry between the excitation of the three plasmon modes.

IV Conclusions

In summary, we have developed a quantum laser theory based on reduced density matrix equation and applied it to a plasmonic nano-laser consisting of a gold nano-sphere and many dye molecules. Our study reveals that the molecular inhomogeneity and the multi-plasmon modes make strong molecular pumping necessary to compensate strong plasmon damping and to achieve lasing. By increasing the molecular density, the plasmon excitation increases, but molecular energy-shifts due to inter-molecular interaction may ultimately reduce the plasmon excitation.

In this article, we modeled the molecular emitters as three-level systems. However, the procedure illustrated can be readily applied to the emitters with arbitrary level structure, which will be necessary to study the influence of other intrinsic processes of the emitters on the laser performance. For example, by introducing more intermediate molecular vibrational levels, in principle, we can study how the intra-molecular vibrational energy redistribution and the temperature of the environment affect the system performance. This extended theory may be utilized to analyze the experiments AYang , where the varying excitation energy of lattice plasmons due to changing the surrounding material affects the dye molecules used by picking up the molecular energy levels resonant to the plasmons. This study will not only provide more insights about the interplay of the plasmons and the gain material but may also suggest how to optimize the system performance.

Appendix A System Parameters and Other Results

In Table 1, we list reference parameters for our simulations. We consider a gold nano-sphere with a radius of $10$ nm. The corresponding dipole plasmons have an excitation energy of $\hbar\omega_{\mathrm{pl}}=2.6$ eV, a damping rate $\hbar\gamma_{\mathrm{c}}=100$ meV and an optical transition dipole moment $d_{\mathrm{pl}}=2925$ D. The classical driving field has the photon energy $\hbar\omega_{0}=2.7$ eV and the amplitude $E_{0}=1.2\times 10^{8}$ V/m consistent with the values used in the experiments WZhou ; AYang ; AYang-1 ; XGMeng ; YJLu-1 ; YJLu ; RMMa ; CYWu ; QZhang ; KDing . The molecules have the transition energy $\hbar\omega_{eg}=2.6$ eV and the transition dipole moments $d_{gf}=16$ D and $d_{ge}=14.4$ D. The molecular transition energy $\hbar\omega_{fg}=2.7$ eV is off-resonant from the higher multipole plasmons YZhang-3 . The decay rate is $\hbar k_{f\to e}^{\left(n\right)}=100$ meV and we assume the other rates $k_{f\to g}^{\left(n\right)}$ and $k_{e\to g}^{\left(n\right)}$ can be ignored compared to the former decay rate.

In Fig. 6 and Fig. 7, we supplement the results presented in the main text with numerical results for system configurations with two and three dipole plasmons.

Appendix B Derivation of Plasmon Reduced Density Matrix Equation

In this section, we derive the master equation for the plasmon reduced density matrix (RDM) $\rho_{\mu\nu}$. From the definition of $\rho_{\mu\nu}$ and Eq. (1), we get the following equation:

[TABLE]

To simiplify notation, we consider $\rho_{\mu\nu}$ as a reference matrix element and denote the dependent matrix elements. For example, $\rho_{\mu_{j}+1\nu_{j}+1}$ differ from the reference element by increasing only the quantum number $\mu_{j}$ and $\nu_{j}$ by one. We see that $\rho_{\mu\nu}$ depends on the molecule-plasmon correlations: $\rho_{a\mu,b\nu}^{\left(n\right)}\equiv\text{tr}_{\text{S}}\left\{\hat{\rho}\left(t\right)\left|b_{n}\right\rangle\left\langle a_{n}\right|\times\left|\nu\right\rangle\left\langle\mu\right|\right\}$, more precisely, $\rho_{e\mu_{j}-1,g\nu}^{\left(n\right)}$, where the labels differ from the ones of $\rho_{e\mu,g\nu}^{\left(n\right)}$ only by subtracting $\mu_{j}$ by unity.

In the procedure to achieve an equation only for $\rho_{\mu\nu}$, the most crucial step is to analyze the equations for $\rho_{a\mu,b\nu}^{\left(n\right)}$ and express them as functions of $\rho_{\mu\nu}$. In the following, we present the equations for the population-like (coherence-like) correlations $\rho_{a\mu,a\nu}^{\left(n\right)}$ ($\rho_{a\mu,b\nu}^{\left(n\right)}$ with $a\neq b$). It turns out that these equations depend on the terms like $\gamma_{j}\sqrt{\left(\mu_{j}+1\right)\left(v_{j}+1\right)}\rho_{a\mu_{j}+1,b\nu_{j}+1}^{\left(n\right)}$ due to the damping from higher plasmon states. To avoid this dependence, we carry out the following replacement in all the equations for $\rho_{a\mu,b\nu}^{\left(n\right)}$:

[TABLE]

Here, we have introduced the complex transition frequency:

[TABLE]

Now, we start with the equations for the population-like correlation $\rho_{a\mu,a\nu}^{\left(n\right)}$ (density matrix elements with same molecular states):

[TABLE]

[TABLE]

[TABLE]

Then, we present the equations for the coherence-like correlations $\rho_{a\mu,b\nu}^{\left(n\right)}$ (density matrix elements with different molecular states) appearing in Eqs. (8), (9) and (10):

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

In the above equations, we have introduced the complex transition frequencies: $\tilde{\omega}_{eg}^{\left(n\right)}=\omega_{eg}^{\left(n\right)}-i\gamma_{eg}^{\left(n\right)}$ with the dephasing rate $\gamma_{eg}^{\left(n\right)}=\left(k_{e\to g}^{\left(n\right)}+k_{e\to f}^{\left(n\right)}+k_{g\to e}^{\left(n\right)}+k_{g\to f}^{\left(n\right)}\right)/2$ and $\tilde{\omega}_{ef}^{\left(n\right)}=\omega_{ef}^{\left(n\right)}+\omega_{0}-i\gamma_{ef}^{\left(n\right)}$ with $\gamma_{ef}^{\left(n\right)}=\left(k_{e\to g}^{\left(n\right)}+k_{e\to f}^{\left(n\right)}+k_{f\to g}^{\left(n\right)}+k_{f\to e}^{\left(n\right)}\right)/2$ as well as $\tilde{\omega}_{gf}^{\left(n\right)}=\omega_{gf}^{\left(n\right)}+\omega_{0}-i\gamma_{gf}^{\left(n\right)}$ with $\gamma_{gf}^{\left(n\right)}\equiv\left(k_{f\to e}^{\left(n\right)}+k_{f\to g}^{\left(n\right)}+k_{g\to e}^{\left(n\right)}+k_{g\to f}^{\left(n\right)}\right)/2$. In addition, we would like to point out that the pure dephasing rate of the emitters can be readily included into these dephasing rates. Because of the coupling with the driving field, we have introduced the following slowly varying correlations $\tilde{\rho}_{e\mu_{j}-1,g\nu}^{\left(n\right)}\equiv e^{-i\omega_{0}t}\rho_{e\mu_{j}-1,g\nu}^{\left(n\right)}$, $\tilde{\rho}_{g\mu,e\nu_{j}-1}^{\left(n\right)}\equiv e^{i\omega_{0}t}\rho_{g\mu,e\nu_{j}-1}^{\left(n\right)}$, $\tilde{\rho}_{f\mu,g\nu}^{\left(n\right)}\equiv e^{i\omega_{0}t}\rho_{f\mu,g\nu}^{\left(n\right)}$ and $\tilde{\rho}_{g\mu,f\nu}^{\left(n\right)}\equiv e^{-i\omega_{0}t}\rho_{g\mu,f\nu}^{\left(n\right)}$.

To proceed, we consider the steady-state equations for the coherence-like correlations (density matrix elements with different molecular states). From Eqs. (13) , (14), (15) and (16) we have

[TABLE]

[TABLE]

[TABLE]

[TABLE]

In order to reduce the dependence, we insert Eqs. (17) and (18) to Eqs. (19) and (20) to express $\tilde{\rho}_{e\mu_{j}-1,f\nu}^{\left(n\right)}$ and $\rho_{f\mu,e\nu_{j}-1}^{\left(n\right)}$ as functions of $\rho_{g\mu,f\nu}^{\left(n\right)}$, $\rho_{e\mu_{j}-1,e\nu_{j}-1}^{\left(n\right)}$ and $\rho_{g\mu,g\nu}$:

[TABLE]

[TABLE]

with the abbreviations:

[TABLE]

[TABLE]

Then, we consider the steady-state version of Eqs. (11) and (12):

[TABLE]

[TABLE]

Finally, we insert Eqs.(21) and (22) to Eqs.(27) and (28) to express $\tilde{\rho}_{g\mu,f\nu}^{\left(n\right)}$ and $\tilde{\rho}_{f\mu,g\nu}^{\left(n\right)}$ as functions of the population-like correlations:

[TABLE]

[TABLE]

with the abbreviations

[TABLE]

Since $\tilde{\rho}_{e\mu_{j}-1,f\nu}^{\left(n\right)}$ and $\rho_{f\mu,e\nu_{j}-1}^{\left(n\right)}$ depend on $\tilde{\rho}_{g\mu,f\nu}^{\left(n\right)}$ and $\tilde{\rho}_{f\mu,g\nu}^{\left(n\right)}$ through Eqs. (21) and (22), we can also express the former two correlations utilizing Eqs. (29) and (30) as functions of the population-like correlations:

[TABLE]

[TABLE]

with the abbreviations

[TABLE]

Finally, since $\rho_{g\mu,e\nu_{j}-1}^{\left(n\right)}$ and $\rho_{e\mu_{j}-1,g\nu}^{\left(n\right)}$ depend on $\tilde{\rho}_{e\mu_{j}-1,f\nu}^{\left(n\right)}$ and $\rho_{f\mu,e\nu_{j}-1}^{\left(n\right)}$ through Eqs. (17) and (18), we can express them also as functions of the population-like correlations:

[TABLE]

[TABLE]

Our next step is to obtain equations only for the population-like correlations. However, before doing so, it is helpful to consider the following combination of terms appearing in Eqs. (8) and (10):

[TABLE]

[TABLE]

with the abbreviations

[TABLE]

[TABLE]

We also consider the combination appearing in (9):

[TABLE]

with the abbreviations

[TABLE]

Now, we consider the steady-state version of Eqs. (8), (9) and (10):

[TABLE]

[TABLE]

[TABLE]

In Eq. (57), we have replaced $\rho_{f\mu_{j}-1,f\nu_{j}-1}^{\left(n\right)}$ by $\rho_{\mu_{j}-1\nu_{j}-1}-\rho_{g\mu_{j}-1,g\nu_{j}-1}^{\left(n\right)}-\rho_{e\mu_{j}-1,e\nu_{j}-1}^{\left(n\right)}$ and then approximated $\rho_{g\mu_{j}-1,g\nu_{j}-1}^{\left(n\right)}$ by $\rho_{g\mu,g\nu}^{\left(n\right)}$. In Eqs. (58) and (59), we have replaced $\rho_{e\mu,e\nu}^{\left(n\right)}$ by $\rho_{\mu\nu}-\rho_{g\mu,g\nu}^{\left(n\right)}-\rho_{f\mu,f\nu}^{\left(n\right)}$. Following this treatment, we get the dependence of the correlations and the plasmon RDM shown in Fig. 1 (c) in the main text. To proceed, we insert Eqs. (39) and (40) into Eqs. (57) and (58) and insert Eq. (53) into Eqs. (58) and (59) to get closed equations for the population-like correlations:

[TABLE]

[TABLE]

[TABLE]

where we have introduced

[TABLE]

in Eq. (60).

We notice that Eq. (60) can be rewritten in a matrix-form and the coefficients before $\rho_{e\mu_{j}-1,e\nu_{j}-1}^{\left(n\right)}$ form a coefficient matrix with the elements $M_{\mu\nu}^{\left(ijn\right)}=\delta_{ij}p_{\mu\nu}^{\left(in\right)}-\left(1-\delta_{ij}\right)i_{\mu\nu}^{\left(ijn\right)}.$ We assume the inverse matrix of the coefficient matrix is $q_{\mu\nu}^{\left(jkn\right)}$ and write the solution of Eq. (60) as

[TABLE]

We can also rewrite Eqs. (61) and (62) in a matrix form with the help of Eq. (64):

[TABLE]

with the abbreviations

[TABLE]

The solution of Eq. (74) is:

[TABLE]

[TABLE]

with

[TABLE]

In summary, we have expressed the coherence-like correlations as functions of the population-like correlations, cf. Eqs. (33), (34) , (37) and (38), and the population-like correlations as functions of the plasmon RDM through Eqs. (64), (84) and (85). By inserting those expressions back into Eq. (5), we get an explicit, linear equation for the reduced density matrix $\rho_{\mu\nu}$ of the plasmon modes:

[TABLE]

Here, we have summed the contribution from individual emitter and introduced the abbreviations $\beta_{\mu\nu}^{\left(jk\right)}\equiv\sum_{n=1}^{N_{\mathrm{e}}}\beta_{\mu\nu}^{\left(jkn\right)},\alpha{}_{\mu\nu}^{\left(j\right)}\equiv\sum_{n=1}^{N_{\mathrm{e}}}\alpha{}_{\mu\nu}^{\left(jn\right)}$ as well as $\tilde{\beta}_{\mu_{j}+1\nu_{j}+1}^{\left(jk\right)}\equiv\sum_{n=1}^{N_{\mathrm{e}}}\tilde{\beta}_{\mu_{j}+1\nu_{j}+1}^{\left(jkn\right)}$, $\tilde{\alpha}_{\mu_{j}+1\nu_{j}+1}^{\left(j\right)}\equiv\sum_{n=1}^{N_{\mathrm{e}}}\tilde{\alpha}_{\mu_{j}+1\nu_{j}+1}^{\left(jn\right)}$. The abbreviations are defined as follows:

[TABLE]

[TABLE]

and

[TABLE]

[TABLE]

with

[TABLE]

Appendix C Equation for Population of Plasmon Number State and Molecular States

The diagonal elements of the plasmon RDM can be interpreted as the population $P_{\mu}\equiv\rho_{\mu\mu}$ of plasmon number states and can be obtained by solving Eq. (3) in the main text. There, $\kappa_{\mu}^{\left(j\right)}\equiv-\sum_{n=1}^{N_{\mathrm{e}}}\alpha{}_{\mu\mu}^{\left(jn\right)}$ and $\eta_{\mu}^{\left(jk\right)}\equiv-\sum_{n=1}^{N_{\mathrm{e}}}\beta_{\mu\mu}^{\left(jkn\right)}$ are the molecule-induced plasmon damping and pumping rate respectively. Since they only depend on the diagonal elements of $\alpha{}_{\mu\nu}^{\left(jn\right)}$ and $\beta_{\mu\nu}^{\left(jkn\right)}$, they can be given explicitly as:

[TABLE]

In the above and also following expressions, all the quantities are the diagonal elements of the corresponding quantities appearing in Sec. B, for example $w_{\mu}^{\left(n\right)}\equiv w_{\mu\mu}^{\left(n\right)}$. The quantities in Eqs. (94) and (95) have the following expressions:

[TABLE]

[TABLE]

where we have introduced

[TABLE]

and

[TABLE]

[TABLE]

as well as

[TABLE]

In the steady-state of the systems, the time-derivative is zero in Eq. (3) in the main text and the resulting equation leads to a recursion relation for the population, cf. Eq. (4) in the main text. In the following, we explain the procedure to calculate the population with the recursion relation for three modes (one and two modes follow as special cases). First, we assume a fixed value for $P_{000}$ and use it to calculate the edge elements $P_{\mu_{1}00}$, $P_{0\mu_{2}0}$ and $P_{00\mu_{3}}$ with simplified versions of Eq. (4), cf. Fig. 8 (a) :

[TABLE]

Secondly, we calculate the surface elements $P_{\mu_{1}\mu_{2}0}$ with the known $P_{\mu_{1}00}$ and $P_{0\mu_{2}0}$ according to a simplified version of Eq. (4), cf. Fig. 8 (b):

[TABLE]

Similarly, we can also calculate other surface elements $P_{0\mu_{2}\mu_{3}}$ and $P_{\mu_{1}0\mu_{3}}$ according to simplified versions of Eq. (4):

[TABLE]

Finally, we calculate the body elements $P_{\mu_{1}\mu_{2}\mu_{3}}$ with the known $P_{\mu_{1}\mu_{2}0}$, $P_{0\mu_{2}\mu_{3}}$ and $P_{\mu_{1}0\mu_{3}}$ by applying repeatedly Eq. (4), cf. Fig. 8 (c). The reason why we can achieve the above simplified versions of Eq. (4) is that the rates with $\mu_{j}<0$ will vanish.

To calculate the population of molecular states $P_{g}^{\left(n\right)}=\sum_{\mu}\rho_{g\mu,g\mu}^{\left(n\right)}$, $P_{f}^{\left(n\right)}=\sum_{\mu}\rho_{f\mu,f\mu}^{\left(n\right)}$ and $P_{e}^{\left(n\right)}=\sum_{\mu}\rho_{e\mu,e\mu}^{\left(n\right)}$, we extract $\rho_{g\mu,g\mu}^{\left(n\right)}$,$\rho_{f\mu,f\mu}^{\left(n\right)}$,$\rho_{e\mu,e\mu}^{\left(n\right)}$ from Eqs. (64), (84) and (85) and express them as functions of the plasmon state population $P_{\mu}=\rho_{\mu\mu}$.

[TABLE]

[TABLE]

[TABLE]