  We calculate the quantum state of the plasmon field excited by an ensemble of
molecular emitters, which are driven by exchange of electrons with metallic
nano-particle electrodes. Assuming identical emitters that are coupled
collectively to the plasmon mode but are otherwise subject to independent
relaxation channels, we show that symmetry constraints on the total system
density matrix imply a drastic reduction in the numerical complexity. For
$N_{\text{m}}$ three-level molecules we may thus represent the density matrix
by a number of terms scaling as $(N_{\rm m}+8)!/(8!N_{\rm m}!)$ instead of
$9^{N_{\text{m}}}$, and this allows exact simulations of up to
$N_{\text{m}}=10$ molecules. Our simulations demonstrate that many emitters
compensate strong plasmon damping and lead to the population of high plasmon
number states and a narrowed linewidth of the plasmon field. For large
$N_{\text{m}}$, our exact results are reproduced by an approximate approach
based on the plasmon reduced density matrix. With this approach, we have
extended the simulations to more than $50$ molecules and shown that the plasmon
number state population follows a Poisson-like distribution. An alternative
approach based on nonlinear rate equations for the molecular state populations
and the mean plasmon number also reproduce the main lasing characteristics of
the system.
