6 W$_{ens}$ = ℏ kc 2 Q$_{ens}$, Q$_{ens}$ = ∑ J Γ r J Γ$_{J}$ + Γ nr 0 , (37) dW ( ω ) dω = √ ϵ$_{0}$ ℏ ω$_{0}$ π ∑ jk v$_{j}$v ∗ $_{k}$A$_{jk}$ , (32) v$_{j}$ = − i 2 ∑ k [ 1 ω$_{0}$ − ω + ˆ Σ ] jk e − $^{iϕ$_{k}$}$, (33) dW ( ω ) dω = ℏ kc 4 π Tr [ 1 ω − ω$_{0}$ − ˆ Σ † ˆ Γ r 1 ω − ω$_{0}$ − ˆ Σ ] , (34) γ r µν = N Γ r 3 B$_{µν}$, (38) B$_{µν}$ = 3 N ∑ j e µ $_{j}$e ν j (39) Γ r = Γ r 0 ∣ ∣$^{1 +}$ a$_{s}$α$_{1}$ r 3 0 ∣ ∣$^{2}$$^{,}$ Γ nr = 3Γ r 0 2 k 3 ∑ l b ( l ) s α ′′ l r 2 l +4 0 , Δ = Γ r 0 ( a$_{s}$α ′′ 1 r 3 0 − 3 4 k 3 ∑ l b ( l ) s α ′ l r 2 l +4 0 ) , (35) Γ r µ = N Γ r 3 λ$_{µ}$, (40) W = ℏ kc 2 Γ r Γ + Γ nr 0 = ℏ kc 2 Q, (36) IV. DISCUSSION AND NUMERICAL RESULTS The energy density then takes the form where A$_{jk}$ = Γ r 0 [ ( e$_{j}$ · e$_{k}$ ) − ˜ α$_{1}$K (1) jk − ˜ α ∗ $_{1}$K (1) kj + | ˜ α$_{1}$ | $^{2}$T (1) jk ] . Matrix A$_{jk}$ is not symmetrical, however only its symmet- rical part, equal to Γ r $_{jk}$, contributes to the final expres- sion. The solution of Eq. (8) can be presented as ˆ Σ | J 〉 = (Δ$_{J}$ − i $_{2}$Γ$_{J}$ ) | J 〉 , where Δ$_{J}$ is frequency shift of collective eigenstate J relative to ω$_{0}$ and Γ$_{J}$ is its decay rate. The molecular relaxation can be accounted for by adding to Σ$_{jk}$ a diagonal term, − i $_{2}$δ$_{jk}$ Γ nr 0 . Then, after frequency integration of Eq. (34), the ensemble radiated energy takes the form where Γ r J = 〈 J | ˆ Γ $^{r}$| J 〉 is radiative decay rate of state | J 〉 . In the photonic Dicke effect, superradiant and subra- diant states are eigenstates of the radiative decay matrix Γ 0 jk obtained from the free space Green diadic. Similarly, in the plasmonic Dicke effect, superradiant states are eigenstates of plasmon-enhanced radiative decay matrix Γ r $_{jk}$. Let us illustrate the emergence of plasmon-mediated superradiance for a simple case when all dipoles are at the same distance from NP surface and are oriented normal or parallel to it. Then it is easy to see that the corre- sponding decay matrix, Γ r jk = Γ $^{r}$e$_{j}$ · e$_{k}$ with Γ r given by Eq. (35), has just three nonzero eigenvalues. Indeed, let us introduce new decay matrices as and after averaging out over the initial random phases ϕ$_{j}$ , we finally obtain where trace is taken over the indexes ( jk ). The analysis of this expression and the results of numerical calcula- tions will be presented in the next section. where Let us start with a single dipole located at r$_{0}$ near a metal NP. In this case, the self-energy is a complex num- ber, Σ = Δ − i $_{2}$Γ, where Δ = Δ$_{jj}$ and Γ = Γ $^{r}$+Γnr = Γ$_{jj}$ are single-dipole energy shift and decay rate, respectively. For normal ( s = ⊥ ) and parallel ( s = ‖ ) dipole orienta- tions, using Eqs. (20) and (21), these are given by 20 is 3 × 3 matrix in coordinate space with Tr ˆ B = 3. It is easy to see that Tr [ ( ˆ Γ $^{r}$)$^{n}$] = Tr [ (ˆ γ $^{r}$)$^{n}$] for any integer n , i.e., the N × N matrix Γ r jk has only three non-zero eigenvalues coinciding with those of matrix γ r µν where λ$_{µ}$ ∼ 1 are eigenvalues of B$_{µν}$ . Note that the decay rates of the remaining N − 3 subradiant states vanish in the long wave approximation used here; they acquire finite values in the next order in ( kr$_{0}$ ) $^{2}$. with a$_{⊥}$ = 2, b ( l ) ⊥ = ( l +1) 2 and a$_{‖}$ = − 1, b ( l ) ‖ = l ( l +1) / 2. Note that both terms in Δ are smaller than their coun- terparts Γ r and Γ nr due to plasmon pole in the imaginary part of NP polarizability α ′′ l . Radiated energy of single dipole-NP system, obtained by frequency integration of Eq. (34), is given by Let us turn to non-radiative coupling, described by matrix Γ nr $_{jk}$. Its diagonal elements, Γ nr jj , describe nonra- diative energy exchange between excited dipole and NP plasmon modes with all angular momenta, as indicated by polarizabilities α ′′ l in Eq. (23). The non-diagonal el- ements of Γ nr jk describe a process by which a plasmon nonradiatively excited in the NP by dipole k transfers its energy to another dipole j . In general, due to high- l plas- mons involved in nonradiative coupling, the eigenstates of Γ nr jk are different from those of plasmon-enhanced ra- diative coupling Γ r jk which contains only dipole ( l = 1) plasmon mode. Therefore, the eigenstates of full decay matrix Γ$_{jk}$ = Γ nr jk + Γ nr jk are not pure superradiant and subradiant states but their admixtures. However, the high- l plasmon contribution to Γ nr jk is significant only at where we included internal molecular relaxation rate, Γ nr 0 , into quantum efficiency Q . For N uncoupled dipoles, i.e., for purely diagonal Σ$_{jk}$ = δ$_{jk}$ ( Δ − i $_{2}$Γ ) , Eq. (34) decouples into sum of N independent terms, yielding W$_{ens}$ = NW . In the presence of inter-dipole coupling, the system eigenstates, | J 〉 , are those of the self-energy matrix, Eq. (12). The corresponding eigenvalues are complex,