We explore the electroluminescence efficiency for a quantum mechanical model of a large number of molecular emitters embedded in an optical microcavity. We characterize the circumstances under which a microcavity enhances harvesting of triplet excitons via reverse intersystem-crossing (R-ISC) into singlet populations that can emit light. For that end, we develop a time-local master equation in a variationally optimized frame, which allows for the exploration of the population dynamics of chemically relevant species in different regimes of emitter coupling to the condensed phase vibrational bath and to the microcavity photonic mode. For a vibrational bath that equilibrates faster than R-ISC (in emitters with weak singlet-triplet mixing), our results reveal that significant improvements in efficiencies with respect to the cavity-free counterpart can be obtained for strong coupling of the singlet exciton to a photonic mode, as long as the singlet to triplet exciton transition is within the inverted Marcus regime; under these circumstances, the activation energy barrier from the triplet to the lower polariton can be greatly reduced with respect to that from the triplet to the singlet exciton, thus overcoming the detrimental delocalization of the polariton states across a macroscopic number of molecules. On the other hand, for a vibrational bath that equilibrates slower than R-ISC (i.e., emitters with strong singlet-triplet mixing), we find that while enhancements in photoluminescence can be obtained via vibrational relaxation into polaritons, this only occurs for a small number of emitters coupled to the photon mode, with delocalization of the polaritons across many emitters eventually being detrimental to electroluminescence efficiency. These findings provide insight into the tunability of optoelectronic processes in molecular materials due to weak and strong light-matter coupling.

## I. INTRODUCTION

The strong light-matter coupling (SC) or polaritonic regime reached with organic molecules embedded in confined electromagnetic environments, such as optical microcavities, has attracted considerable attention due to the possibilities it offers to the control of chemical processes at the molecular level.^{1–5} In this regime, the interaction energy between the molecular ensemble and the photonic modes surpasses their respective linewidths, resulting in the formation of hybrid photon-matter excitations termed polaritons.^{6} The novel tunability of chemical dynamics and related processes afforded in the SC regime with organic molecules has been proven in photoisomerization,^{7} excitation energy transfer,^{8–10} optical selection rules,^{11} exciton transport,^{12,13} and more recently triplet harvesting.^{14–16} The SC regime has not been exclusively harnessed at optical frequencies as unique signatures of chemical control^{17–19} and nonlinear responses^{20,21} have also been demonstrated for molecular vibrational transitions coupled to a microcavity field.

Polariton setups have also emerged as promising candidates to boost the efficiency and versatility of light-emitting diodes (LEDs)^{22} and organic photodiodes (OPDs).^{23} As a proof of concept, electrical injection of carriers into inorganic polariton architectures has been shown to modulate light-matter coupling,^{24} thus opening new avenues toward polariton-based optoelectronic switches. Similarly, organic LEDs have been demonstrated utilizing molecular dyes in the ultrastrong coupling regime^{25–27} and it has recently been shown that materials operating in the latter can feature a complete inversion of molecular dark and light-emissive states.^{28} These exciting applications prompt the development of theoretical models to account for experimental observations of the molecular SC regime,^{29–32} as well as to rationally design new polaritonic setups that could enhance chemical processes of contemporary interest.^{33,34}

In this article, we explore the triplet electroluminescence efficiency of a microcavity containing molecules, which feature a range of electronic parameters and couplings to the condensed phase vibrational bath. Our model aims to describe polaritonic OLEDs such as those reported in Refs. 27 and 35, where (optically dark) triplet excitons are generated upon electrical injection and the latter can transition into fluorescent singlet states that emit light. Our approach relies on a master equation operating in a variationally optimized polaron frame, originally introduced by Silbey and Harris^{36,37} with generalizations due to Pollock,^{38} Wu,^{39} and their respective co-workers. The foundation of this method lies on the application of a unitary transformation to the total Hamiltonian, which yields a renormalized system and system-bath interactions that are weak enough to be perturbatively treated.

The outline of this article is summarized as follows. In Sec. II, we introduce our quantum mechanical model and the variational approach to address its dynamics. Next, in Sec. III, we give a formal definition of the triplet electroluminescence efficiency and describe our approach to calculate the time-evolution of populations in the relevant chemical states of the system. Then, in Sec. IV, we apply our approach to systems in two regimes of coupling to vibrational degrees of freedom (DOFs). In Sec. V, we identify the main limitation that must be overcome to reach a polariton-enhanced triplet-harvesting regime and propose an approach to circumvent this drawback. Finally, in Sec. VI, we summarize our study and conclude with an outlook of how polaritons could enhance the optoelectronic properties of molecular materials.

## II. THEORY

We consider an ensemble of *N* identical molecules embedded in an optical microcavity and interacting with the electromagnetic modes supported by the latter; for simplicity, we describe a single photonic mode interacting with the molecular ensemble, a coarse-graining approximation which is based on the much larger density of states (DOS) of the molecular degrees of freedom compared to the photonic ones.^{29,40} Thus, *N* should not be interpreted as the total number of molecules in the cavity but rather as the average number of molecules that couple to a single photon mode.

The Hamiltonian of our model can be written as

Each molecule is modeled as a three-level electronic system, namely, a singlet electronic ground state and two excited states with singlet and triplet spin characters, respectively. The energetics of the latter are correspondingly described by *H*_{S} and *H*_{T}. Assuming electrical pumping in the linear regime, we can restrict the Hamiltonian to the single excitation manifold such that

Here, |*n*⟩ (|*T*_{n}⟩) denotes a localized singlet (triplet) exciton at the *n*th molecular site, in the absence of microcavity photons, and *ϵ*_{S} (*ϵ*_{T}) is the vertical singlet (triplet) electronic excitation energy, while (*ℏ* = 1)

accounts for the vibrational degrees of freedom in the condensed phase environment, where $bv,n\u2020$ (*b*_{$v$,n}) is the creation (annihilation) operator for an excitation in the $v$th harmonic mode with frequency *ω*_{$v$} on site *n*. The mode-dependent singlet (triplet) vibronic couplings $gS(v)$ ($gT(v)$) are included in *H*_{S–B} (*H*_{T–B}),

The singlet-triplet intersystem crossing electronic coupling is given by

Finally, the photonic microcavity degree of freedom is encoded in

where |*G*, 1_{ph}⟩ accounts for the state of all molecules in the electronic ground state and one excitation in the photonic mode whose interaction with singlet excitons (in the rotating wave approximation, RWA) is described by

*g* being a single-molecule dipolar coupling. Notice that the truncation of the electronic and photonic degrees of freedom introduced in Eq. (7) is aligned with our single-excitation manifold restriction, $\u27e8a\u2020a\u27e9+\u27e8IS+IT\u27e9=1$, where $IS=\u2211n|n\u2009\u2009n|$, $IT=\u2211n|Tn\u2009\u2009Tn|$, and *a*^{†} (*a*) is the cavity-photon creation (annihilation) operator. Furthermore, in the Hamiltonian in Eq. (1), we ignore intersite couplings, which significantly simplifies our model and at the same time allows for the exploration of the dynamics induced purely by the photonic field. Our approximation is justified if intersite hopping is the smallest energy scale of the problem, significantly weaker than the collective light-matter interaction.

To describe the emergent dynamics upon electrical pumping of the optically dark triplet states, we employ a variational polaron transformed master equation. The latter has proven to reproduce reliable quantum dynamics in a wide range of coupling strengths between electronic degrees of freedom and phononic reservoirs,^{41,42} while being computationally inexpensive compared to more accurate approaches such as path integral,^{43,44} multiconfiguration time-dependent Hartree,^{45,46} hierarchical equations of motion,^{47} linearized density matrix evolution,^{48} surface-hopping,^{49–51} and exact factorization^{52} formulations. In particular, we consider a similar multisite approach to the one developed by Pollock and co-workers,^{38} who considered a variational polaron transformation to compute population dynamics in exciton networks with local phonon reservoirs. Furthermore, we adapt a generalization introduced by Wu and co-workers^{39} who treated the photon-exciton and exciton-vibration on equal footing to describe the photoluminescence and vibrational dressing in polaritonic setups. As a first step toward the development of a master equation, the Hamiltonian in Eq. (1) is transformed to the so-called polaron frame, for which we introduce the unitary transformation *e*^{P}, where

is partially based on previous works.^{39,53} According to the ansatz in Eq. (8), the electronic and photonic degrees of freedom are dressed with vibrational bath excitations to an extent quantified by the set of parameters ${\u2009fl(v),fT(v),h(v)}$ which are variationally determined, as will be explained below. The summation over *l* in Eq. (8) is introduced to account for the photon-induced nuclear deformations across many molecules^{39,53} such that the vibrational disentanglement between the electronic and photonic degrees of freedom of the polariton modes is reproduced for large enough light-matter interaction.^{31} Since the electronic triplet states do not directly interact with the photonic mode, the vibrational deformation of the triplet excited manifold is expected to be more localized (small polaron limit). In the polaron frame, we have

where

and

In Eqs. (10) and (11), we introduced a delocalized Fourier basis for the singlet $|k\u2009=1N\u2211neikn|n\u2009$ and triplet excitons $|Tk\u2009=1N\u2211neikn|Tn\u2009$, with $k=2\pi N$, *m* = 0, 1, 2, …, *N* − 1. In this basis, the state |*k* = 0⟩ couples to the photonic mode with a superradiantly enhanced strength given by $Ng=\Omega /2$. The renormalized on-site energies in (10) are given by

and the renormalization constant *η*_{ST} (*η*_{S–ph}) is the thermal average of the relative displacement operator between the singlet and triplet (photonic) harmonic potential energy surfaces (PESs),

where $\u27e8A\u27e9=TrB{Ae\u2212\beta HB}ZHB$, Tr_{B}{·} denoting a trace over the vibrational degrees of freedom, *β* being the inverse temperature, and $ZHB=TrB{e\u2212\beta HB}$ being the vibrational partition function. With these definitions, note that $H\u03030$ = $ePHe\u2212P\u2212HB$ [see Eqs. (9) and (10)], that is, it corresponds to a vibrationally thermally-averaged polariton Hamiltonian. Notice that by taking $fl(v)=gS(v)\delta l,0$ in Eq. (12a), we recover the full-polaron transformation for the singlet excitation, and $\u03f5\u0303S\u2212\u03f5S=\u2212\lambda S=\u2212\u2211v(gS(v))2/\omega v$, where λ_{S} is the reorganization energy of the singlet excited state. In Eqs. (10) and (11), we also introduced a partition of the Hamiltonian into the molecular contribution that is totally symmetric under site permutations $H\u03030,k=0$ and the nontotally symmetric one $\u2211k\u22600H\u03030,k$. Since $[H\u03030,k=0,\u2211k\u22600H\u03030,k]=0$, the eigenstates of $H\u03030$ can be found by independent diagonalization of each contribution,

where we define the upper (UP), middle (MP), and lower (LP) polariton states $|j\u2009=c\u2009Sj|k=0\u2009+cTj|Tk=0\u2009+cphj|G,1ph\u2009$, *j* = UP, MP, LP. On the other hand,

describes mixed singlet-triplet eigenstates $|k(\xb1)\u2009=cS\xb1|k\u2009+cT\xb1|Tk\u2009$ with no photonic component, whose eigenenergies are given by $\u03f5\u0303\xb1=\u03f5\u0303S+\u03f5\u0303T2\xb112(\u03f5\u0303S\u2212\u03f5\u0303T)2+4VST2\eta ST2$. These states are commonly known as dark or exciton reservoirs.^{54,55} For clarity, we refer to the highest (lowest) energy eigenstates in Eq. (14) as the upper (lower) dark states. Notice that since we are ignoring intersite couplings, we neglect the dispersive character of the exciton in the calculation of the eigenenergies (14), a valid assumption in view of the flat exciton dispersion relation compared to the photonic one in the *k* range of interest. The residual interaction Hamiltonian in Eq. (9) is

which has been written as a sum of different contributions defined as follows:

The calculation of the parameters ${\u2009fl(v),fT(v),h(v)}$ is carried out by taking advantage of the Feynman-Bogoliubov inequality $F\u2264F0\u2032+\u27e8H\u0303I\u27e9H\u03030\u2032$,^{36,37} where *F* ($F0\u2032$) is the free energy of the system governed by *H* ($H\u03030\u2032=H\u03030+HB$) and $\u27e8A\u27e9H\u03030\u2032=Tr{Ae\u2212\beta H\u03030\u2032}[Tr{e\u2212\beta H\u03030\u2032}]\u22121$, where Tr{·} denotes the trace over the eigenstates of $H\u03030\u2032$. Notice that by construction, $\u27e8H\u0303I\u27e9H\u03030\u2032=0$, and the leading correction to the exact equilibrium reduced density matrix (RDM) is $O(H\u0303I2)$,^{42} which justifies the (second order) perturbative treatment in $H\u0303I$. It follows that the parameters ${\u2009fl(v),fT(v),h(v)}$ can be found by minimizing $F0\u2032=\u2212\beta \u22121Tr{e\u2212\beta H\u03030}$, which amounts to solving $\u2202F0\u2032\u2202fi(v)=\u2202F0\u2032\u2202h(v)=0$. In our calculations, we consider a continuum vibrational bath limit, which is described in terms of the spectral densities $Ji(\omega )=\u2211v|gvi|2\delta (\omega \u2212\omega v)=ai\omega 3\omega c2e\u2212\omega /\omega c$, *i* = *S*, *T*, where *ω*_{c} is the cutoff frequency and *a*_{i} is the dimensionless parameter that encodes the strength of coupling between the excited electronic state with character *i* to the vibrational bath. This choice of spectral density qualitatively describes the behavior of low-frequency phonon modes in three-dimensional solids although our results should be qualitatively the same for other choice of spectral densities as long as the cutoff frequency is finite and does not lead to the infrared divergence that occurs when *s* ≤ 1. The computational details of the solution for ${\u2009fl(v),fT(v),h(v)}$ are summarized in Appendix A of the supplementary material. In the discussion that follows, we focus on the calculation of the population dynamics of the chemically relevant species as well as of the electroluminescence efficiencies.

## III. DYNAMICS IN THE POLARON FRAME AND DEFINITION OF TRIPLET ELECTROLUMINESCENCE EFFICIENCY

The open-quantum system dynamics associated with the Hamiltonian in Eq. (1) can be described in terms of the time evolution of the reduced density matrix (RDM) *ρ*(*t*) = Tr_{B}{*ρ*_{tot}(*t*)}, [*ρ*_{tot}(*t*) is the total density matrix], governed by the Liouville equation in the Schrödinger picture,

where $L\rho (t)=\u2212iTrB{[H,\rho tot(t)]}$ and {*ρ*(*t*), Γ_{i}} = *ρ*(*t*)Γ_{i} + Γ_{i}*ρ*(*t*); Γ_{nrad} and Γ_{ph} phenomenologically account for the dissipative processes associated with the nonradiative decay of excitons^{56} and photonic leakage, respectively. We define the triplet electroluminescence efficiency as

where the trace is taken over all the degrees of freedom, namely, the electronic, vibrational, and photonic. Equation (17) is analogous to the integrated probability used to define the efficiency of energy trapping in chromophoric complexes^{56} but in the present context acquires the meaning of the efficiency of emission of a photon. Here, we define $\Gamma nrad=knradS2\u2211n|n\u2009\u27e8n|+knradT2\u2211n|Tn\u27e9\u2009Tn|$ and $\Gamma ph=\kappa 2|G,1ph\u2009\u2009G,1ph|$, $knradi$ being the nonradiative decay rate of electronic state *i* = *S*, *T*, and *κ* being the cavity photon leakage rate. Since

*ϵ* is an invariant quantity under the polaron transformation. This last condition permits the computation of Eq. (17) in the polaron frame, where the time evolution of the polaron transformed RDM $\rho \u0303(t)=TrB{\rho \u0303tot}$ can be described in terms of second-order perturbation theory on $H\u0303I$ within the secular Born-Markov approximation. This procedure guarantees that the long-time evolution of the polaritonic system properly thermalizes into the reduced equilibrium state due to $H\u03030$ which, due to the optimization of the variational parameters according to the Feynman-Bogoliubov bound (see Sec. II), provides a good description of the equilibrium state of the full system involving electronic states, vibrations, and the photonic mode.

For the calculation of *ϵ*, we assume that the initial RDM corresponds to an incoherent mixture of localized excitations in the triplet electronic manifold, i.e.*,* $\u2211n\u27e8Tn|\rho \u0303(0)|Tn\u27e9=1$, $\u2211n\u27e8Tn|\rho \u0303(0)|n\u27e9=\u2211n\u27e8Tn|\rho \u0303(0)|G,1ph\u27e9=\u2211n\u27e8n|\rho \u0303(0)|n\u27e9=\u27e8G,1ph|\rho \u0303(0)|G,1ph\u27e9=0$*.* Under these assumptions, the initial RDM in the polaron frame can be written as

where the approximation is the result of considering a localized initial state, which in the delocalized polaron picture translates into the population distributed uniformly among the dark and polariton states. Since the former feature a larger DOS than the latter, the population is predominantly concentrated in the dark states. We can rewrite Eq. (17) as

where the vector

accounts for the photonic character of the different chemical species that take part in the population dynamics. The time evolution of $\rho \u0303\u2192(\tau )$ is described in terms of a Pauli master equation in the secular Markovian approximation, where

For simplicity in the calculations and interpretation, we consider the evolution of the total population in the dark state manifolds $\u2211k\u22600\u27e8k\xb1|\rho \u0303(t)|k\xb1\u27e9$ rather than dissected among the individual populations ${\u27e8k\xb1|\rho \u0303(t)|k\xb1\u27e9}$. The matrix $M$ is given by

where *k*_{j←i} is the rate of transfer from the state/manifold |*i*⟩ to |*j*⟩, and we remind the reader that $IT=\u2211n|Tn\u2009\u2009Tn|$, $IS=\u2211n|n\u2009\u2009n|$. Details of the calculation of the population transfer rates are included in Appendix B of the supplementary material. Using Eq. (21) together with Eq. (20), we have that

To get further insight into the contributions of the different pathways of population transfer to *ϵ*, we consider a Green’s function partition approach,^{56} under which

where we let $M=Mpol+Mnpol$, $Mpol$ being the rate matrix that accounts for the feed and depletion of polariton states only, which amounts to considering the rates *k*_{i←j} where either *i* or *j* ∈ {|*UP*⟩, |*MP*⟩, |*LP*⟩} while setting the other entries of $M$ to 0. Using Eqs. (23) and (24),

where $\u03f5dir=\kappa PphMnpol\u22121MpolM\u22121\rho \u0303\u2192(0)$, $\u03f5indir=\u2212\kappa PphMnpol\u22121\rho \u0303\u2192(0)$, and *ϵ* can be written as a sum of a contribution due to polariton-participating processes and another one where they do not play a role.

## IV. RESULTS AND DISCUSSION

Based on the approach above, we are particularly interested in the dependence of *ϵ* with Rabi energy Ω and detuning between the cavity photon and the vertical singlet energy transition (Δ = *ω*_{ph} − *ϵ*_{S}), considering molecular emitters which feature different regimes of interaction with the vibrational bath degrees of freedom (DOFs). For that purpose, we introduce two cases based on parameters chosen to represent a wide range of organic molecular emitters.

First, we consider (large) strengths of coupling to the vibrational environment *a*_{S} = 10, *a*_{T} = 15, and a frequency cutoff *ω*_{c} = 0.01 eV. The corresponding reorganization energy of the singlet electronic state is given by $\lambda S=\u222b0\u221eJS(\omega )\omega d\omega =0.2$ eV, and the analogous quantity for the triplet is $\lambda T=\u222b0\u221eJT(\omega )\omega d\omega =0.3$ eV. The assumption of a single vibrational reservoir coupled to both the singlet and triplet states allows for the calculation of the reorganization energy of the singlet-triplet transition λ_{ST} (see Appendix A of the supplementary material for details), which for this specific case amounts to 0.01 eV. For simplicity, we do not consider high frequency vibrational bath modes (those which feature frequencies *ω* ≫ *β*^{−1}) in our calculations. Furthermore, we assume an energy gap between the equilibrium vibrational configurations of the triplet and singlet Δ*G* = (*ϵ*_{T} − λ_{T}) − (*ϵ*_{S} − λ_{S}) = −0.1 eV and a weak intersystem crossing coupling amplitude *V*_{ST} = 2 × 10^{−5} eV. These parameters are typical for organic molecules that undergo thermally activated delayed fluorescence (TADF),^{57–59} where the population in the dark triplet states transfers to the bright singlets, with a reverse intersystem-crossing (R-ISC) rate that depends on a thermal energy barrier between the latter. We also consider a relatively slow nonradiative decay of the triplet $knradT=2\xd710\u22127$ ps^{−1}.^{58} Since the characteristic time scale of relaxation of the vibrational environment $\omega c\u22121\u226a\u210f/|VST|$, we expect a fast relaxation of the latter into its equilibrium configuration before R-ISC ensues, a scenario termed “fast bath” in the literature.^{42} This is the behavior expected in the limit Ω → 0 (small polaron formation before electronic transition) which was confirmed in our numerical calculations.

The second case we considered is the opposite to the previous one, in which the vibrational DOFs of the environment are sluggish in comparison with the R-ISC transition time scale (*V*_{ST} ≫ *ω*_{c}), a scenario identified as “slow bath,”^{42} which corresponds to molecular emitters with sizable singlet-triplet mixing via spin-orbit coupling or hyperfine interaction. For this case, we also chose large strengths of coupling to the vibrational environment *a*_{S} = 10, *a*_{T} = 30, while setting a small cutoff frequency *ω*_{c} = 5 × 10^{−3} eV (which translate into λ_{S} = 0.01 eV, λ_{T} = 0.03 eV, and λ_{ST} = 0.1 eV). Furthermore, we set *V*_{ST} = 0.05 eV and Δ*G* = −0.01 eV, which can qualitatively describe organometallic complexes exhibiting ultrafast ISC rates {such as [Cu(dmp)_{2}]^{+}^{60}}, with a slow-relaxing vibrational environment. For this case, if $knradT$ is small as above, the cavity-free triplet electroluminescence efficiency *ϵ* ≈ 1 and placing it in the cavity is not productive; hence, we assume that the molecule is such that it features a fast nonradiative triplet decay $knradT=2\xd710\u22123$ ps^{−1} such that we explore the possibility to outcompete it by cavity-assisted processes. The temperature *T* = 300 K is assumed for both cases.

For reference, we computed electroluminescence efficiencies in the cavity-free (bare molecules) case which are *ϵ* = 0.027 for the fast bath case, and *ϵ* = 0.062 for the slow one. The latter are calculated using Eq. (17), substituting Γ_{ph} with $\Gamma rad=kradS2IS$, where $kradS=10\u22122$ ps^{−1}, a typical fluorescence rate for organic molecules.

Upon coupling to a cavity photon mode, the emergent dynamics are changed as a result of the modification of the bare molecular energy landscape and the reorganization energies between the polariton states. In the polaron frame, the former contribution is encoded in the energy spectrum introduced in Eqs. (13) and (14), whereas the latter is accounted for in the reorganization energies λ_{i−j}, *i* ≠ *j*, *i*, *j* ∈ *UP*, *MP*, *LP*, ±, where the subindices ± denote |*k*^{±}⟩ states (in our model, the reorganization energies λ_{i−(±)} are independent of the momentum of the dark states). The polaron frame picture of the dynamics in both the slow and the fast bath case is schematically illustrated in Fig. 1, and the variation of the electroluminescence efficiency with Δ and Ω for the two different molecular emitter cases is shown in Fig. 2.

For the discussion in this section, we limit ourselves to the *N* = 10^{2} case, postponing the important analysis of larger *N* values for the next paragraphs.

*Fast bath scenario.* For the fast bath molecule, we find that for weak light-matter couplings (Rabi splittings Ω ≪ λ_{S}), *ϵ*(Δ, Ω) is below the cavity-free scenario [for the resonant case Δ = 0, see the inset in Fig. 3(a)]. In the polaron frame, this observation can be explained in terms of thermally activated transfer of population from the triplets to the photonic potential energy surface (PES) via the singlets [see Fig. 1(a-I)]. Since this process involves passage through two thermal energy barriers, the time needed to reach the photonic PES is longer than the triplet lifetime, and therefore, *ϵ* ≪ 1. On the other hand, when Ω/2λ_{S} ≈ 2/3, our calculations predict an abrupt increase in *ϵ* [compare the inset and main plot in Fig. 3(a)], which is due to a transition to a regime where the dynamics in the polaron frame are described as depicted in Fig. 1(a-II). Under this picture, polaron decoupling^{32,53} manifests, which is expected for sufficiently large Rabi splittings,^{31,53} under which the polariton PESs are essentially undisplaced with respect to the electronic ground PES. The intuition behind this decoupling is that the exchange of energy between photon and singlet excitons is much faster than that between vibrations and singlet excitons. Importantly, while this regime was previously predicted for Ω ≫ λ_{S}^{61} for a single high-frequency mode, we hereby numerically demonstrate that this stringent condition can be relaxed to smaller Rabi splittings in the case of a continuum of low-frequency modes. In this polaron decoupling regime, our calculations predict that λ_{(−)−UP} = λ_{(−)−MP} = λ_{(−)−LP} = λ_{(+)−(−)} > λ_{ST} [where the equality among the previous quantities is a result of the ansatz in Eq. (8), which for Ω/2λ_{S} > 2/3 corresponds to $fl(v)=f(v)=O(gS(v)N)$, for all *l*. In other words, the singlet exciton PESs are negligibly displaced with respect to the ground state PES, when *N* is large]. As a consequence, the R-ISC energy barrier can either decrease or increase with respect to the bare molecules by tuning the Rabi splitting. In fact, we notice that as Ω increases [see Figs. 2(a-i) and 1(a-II)], *ϵ* increases and then decreases with respect to the cavity-free value, given that the smallest R-ISC to polariton energy barrier decreases and then increases as this rate goes from normal to inverted Marcus regimes. Notice that Fig. 3(a) only shows an increasing behavior in *ϵ* because the plotted range of Rabi splittings is small.

*Slow bath scenario.* For the slow bath case, we observe a significantly different *ϵ*(Δ, Ω). In Fig. 2(b–i), we show that there is no need to invoke Rabi splittings beyond the reorganization energy to observe an enhancement in electroluminescence efficiency.

Similarly to the previous fast bath scenario, our calculations reveal two qualitatively different pictures of the population dynamics in the polaron frame depending on the magnitude of Ω. For sufficiently weak light-matter coupling, the initial state is comprised of approximately equally populated lower and upper dark state manifolds since $|cS+|2\u2248|cS\u2212|2$ [see Eq. (19) and Fig. 1(b-I)], a consequence of the prevalence of the electronic coupling *V*_{ST} over the weak coupling of the R-ISC transition to the phonon environment. In this regime, *ϵ* can be explained in terms of a Marcus picture where the photonic PES receives population from the upper and lower dark-state manifolds with a rate proportional to a diabatic coupling of order |*g*|^{2}, corresponding to single-molecule light-matter coupling [see Eq. (16e) and Fig. 1(b-I)]. High and small R-ISC energy barriers can be obtained by scanning across Δ, thus decreasing and increasing *ϵ*, respectively; for the computed range of Δ values, the dominant R-ISC pathway is the one starting from the upper dark states [they are closer in energy to the photonic PES, see Fig. 1(b-I)]. In fact, for a sufficiently small increase in Ω, the triplet electroluminescence efficiency *ϵ* also rises because |*g*|^{2} increases; see Fig. 3(b), which shows a cut of *ϵ* at light-matter resonance Δ = 0. In this plot, this mechanism is operative up to a first maximum in *ϵ* at Ω ≈ λ_{S}, upon which light-matter coupling changes from weak to moderate and *ϵ* decreases. This effect can be understood as follows: as Ω increases, it can compete with *V*_{ST}, thus reducing the magnitude of the dressed $V\u0303ST$, and concomitantly reducing the mixing between singlets and triplets. The result is that the lower and upper dark states become more tripletlike and singletlike, respectively. Therefore, the initial triplet population has more probability of residing in the lower dark states where they feature a much larger energy barrier to the photonic mode, thus leading to slower R-ISC and lower *ϵ*.

Keeping our discussion around Fig. 3(b), our model predicts a sharp change in efficiencies under strong light-matter coupling starting at Ω ≈ 4λ_{S}, with *ϵ* increasing with Ω, featuring a second maximum at Ω ≈ 6λ_{S}, and then decreasing for larger Ω. In analogy to the fast bath case above, this abrupt change corresponds to a transition to a different regime of population dynamics, pictured in Fig. 1(b-II). An intriguing feature of this regime is that, for the chosen parameters, the upper and lower dark states have largely triplet and singlet characters, respectively; that is, the configuration of minimal energy of the vibrational DOFs coupled to the singlet and triplet states in the polaron frame is such that their energy ordering is inverted with respect to that in the original frame [Fig. 1(b-II)], where *ϵ*_{S} − λ_{S} > *ϵ*_{T} − λ_{T}. This effect is due to the slow relaxation of the bath that renders the vibrational dressing of the states at nonequilibrium nuclear configurations close to those accessed by vertical transitions from the ground state. Under these circumstances, population is primarily initialized in the upper dark states and transfers to the MP [Fig. 1(b-II)] which is closer in energy than the other polaritons or dark states. Notice that in contrast to the weak light-matter coupling case, MP has significant singlet, triplet, and photonic character, and the population transfer mechanism can be essentially dissected into sizable contributions from multiphonon (a Marcus-like contribution due to displacement between the upper dark states and MP harmonic potentials) and one-phonon (Redfield-like) processes. In our example, a steady boost in *ϵ* is given by the Redfield processes which dominate over the Marcus-like ones, given the large activation energy barriers to be surmounted in the latter. These Redfield processes are primarily due to the triplet character in both upper dark states and MP. After the second maximum in *ϵ*, a further increase in Ω leads to a phonon blockade and a concomitant decrease in *ϵ*: the MP energy keeps lowering, but no one-phonon process is available to mediate the transfer from the upper dark states (the vibrational DOS decreases exponentially as a function of frequency in our chosen spectral density).

To conclude this section, we further comment on the nonsmooth behavior of efficiencies displayed in Figs. 2 and 3. The presence of abrupt changes is due to the existence of more than one solution to the self-consistent equations during the variational calculation: for the fast bath case, one solution describes a scenario where the singlet state is localized [see Fig. 1(a-I)], while the second solution corresponds to a (delocalized) polariton [see Fig. 1(a-II)]. Similarly, for the slow-bath scenario, the two solutions correspond to mixed singlet-triplet states localized in one molecular site, and to a polariton picture, respectively [see Figs. 1(b-I) and 1(b-II), respectively]. Following the variational prescription, we choose the solution that features the lowest free energy $F0\u2032$. As it has been pointed out in previous studies, a slightly more accurate description of this crossover from a localized to a delocalized regime (which is devoid of these abrupt transitions) can be described by including higher order perturbative corrections^{62} or by utilizing an improved variational ansatz;^{63} these small corrections are beyond the scope of this article and do not affect the main conclusions of our work.

## V. THE LARGE N ISSUE

So far, we have addressed the *N* = 100 case in detail. In fact, strong light-matter coupling with small delocalization lengths as those addressed here has been reported using plasmonic nanoparticle environments.^{64,65} However, most microcavity systems have much larger *N* values. We thus proceed to comment on the *ϵ* behavior for increasing values of *N*, while keeping Rabi splittings Ω fixed (and hence energy barriers and energy differences between states, see Fig. 2). For weak light-matter coupling, both fast (a) and slow (b) baths yield insensitive changes of *ϵ* with increasing *N*. This is due to the fact that changing *N* does not alter the energy barriers to be surmounted in the R-ISC process.

For strong light-matter coupling, the overall electroluminescence efficiency is damped as a result of the ratio of triplet to polariton states increasing. For the fast bath case (a), since both dark and polariton states feature the same reorganization energies in the polaron frame, *ϵ* is determined by the competition between the nonradiative triplet decay and transfer to the polariton state closer in energy, which is because the latter channel features a lower activation energy barrier compared to the dark states [see Fig. 1(a-II)]. However, at light-matter resonance, the R-ISC rate to the final polariton state scales as $1N$, given that the polariton is delocalized across *N* singlets and only one of them can undergo coupling to a given triplet (see Fig. 6). Importantly, for *N* ≥ 10^{4}, the rate of transfer to the polaritons falls below the cavity-free R-ISC rate for the entire (Δ, *ω*) range displayed [see Fig. 2 top panels and Fig. 4(a)]. The conclusion is similar for the slow bath case (b). While for small *N* values, high efficiencies are obtained [see Figs. 2(b-i) and 2(b-ii)], the assistance of polariton modes decreases with increasing *N*. This is because the ratio of rates (dark states → polariton)/(polariton → dark states) scales as $1N$ given the delocalization of the polariton state; this implies that the backward process polariton → dark states becomes more relevant as *N* increases, an effect which is detrimental to *ϵ*. Finally, the decrease in photon character in the MP for Δ < 0 [see Fig. 1(b-II)] together with a fast back transfer to the denser manifold of upper dark states for larger *N* explains the drop of *ϵ* at Δ < 0, specially for *N* > 10^{4} [see Fig. 2 lower panels and Fig. 4(b)].

Based on the discussion above, a pressing question is the characterization of an optimal scenario that overcomes the deleterious effects of cavity R-ISC rates as *N* increases, the latter being as large as 10^{7}.^{40} Here, we propose a case where the energy tunability of the polaritons introduces an increase of many orders of magnitude relative to the cavity-free R-ISC rate, and this enhancement can outcompete the suppression due to the low polariton DOS. The scenario is built upon the fast bath case studied above, but assuming a sufficiently large singlet-triplet energy gap and weak couplings of the electronic states to the bath (small *a*_{S} and *a*_{T} parameters). Under these conditions, the reduced dynamics of transfer in the polaron frame in the cavity-free scenario corresponds to the inverted-Marcus regime (for the singlet to triplet transition), where the population of the triplet manifold needs to overcome an enormous energy barrier to reach the singlet state (see Fig. 5). Upon confinement to a microcavity, the emergent lower polariton could introduce a much smaller energy barrier from the triplet states. Our results and the precise parameters used in our calculations are shown in Fig. 5, which show that polariton-assistance introduces four orders of magnitude enhancement relative to the bare R-ISC rate.

On the other hand, we note that the polariton-assisted R-ISC rate is upper-bounded by $VST2N\pi \beta \lambda T-pol$, which must be higher than $knradT$ to lead to harvesting of triplets. If we take *N* = 10^{6} and use the parameters above for the fast bath scenario (*k*_{nrad} = 2 × 10^{−7} ps^{−1}), we need a very small λ_{T–pol} = 10^{−17} eV. Therefore, in order to obtain a sizable *ϵ*, emitters with negligible couplings to a vibrational environment are required (small reorganization energies λ_{ST} and λ_{S}); that is, poor R-ISC molecules (i.e., the opposite of TADF materials) will obtain most benefit from polaritonic effects. Whether this extreme possibility exists will be subject of future work.

For the slow bath case, there is no possibility to enhance electroluminescence for large *N* because the associated R-ISC processes via polaritons happen via one-phonon vibrational relaxation (Redfield) processes, rather than multiphonon (Marcus) processes that exhibit exponential sensitivity on energy barriers. Hence, the only option for enhancement of electroluminescence for that case is to use samples with small *N* values.

## VI. CONCLUSIONS

In this article, we have developed a variational model for the prediction of electroluminescence efficiency of molecular emitters interacting with a photonic mode of a microcavity. The latter can describe the emergent dynamics on a range of Rabi splittings that extend from the weak to the strong coupling regime.

Our calculations indicate that the main limitation to polariton-assist the electroluminescence process is the delocalized character of polaritons: the probability of a triplet state to R-ISC into the latter is strongly suppressed by a $1N$ factor, similarly to results discussed in previous works.^{9,28,34} There are two approaches that could overcome the detrimental delocalization effect and therefore introduce a polariton-enhanced electroluminescence with respect to the cavity-free scenario: (1) for fast and slow bath cases (i.e., for chromophores with weak and strong singlet-triplet mixing), as Figs. 2(a-i) and 2(b-i) show, strong coupling with a small number of molecules per photonic mode *N* introduces significant electroluminescence yield enhancements with respect to the bare molecular case. This regime should be feasible in the context of plasmonic nanoparticles coated with organic dyes, where *N* can be more than 3 orders of magnitude smaller than those needed in a microcavity to attain the strong coupling regime.^{64–66} (2) The second approach only works for the fast bath case and consists of tuning of polariton energies and taking advantage of the exponential sensitivity of R-ISC rates with respect to energy barriers. Therefore, a high thermal energy barrier for the R-ISC transition in the bare molecular case translates into the possibility of its partial or total reduction via a triplet to polariton R-ISC transition, thus potentially outcompeting the delocalization effect, and introducing a net electroluminescence enhancement. For a molecular emitter, the previous requirement is equivalent to an inverted Marcus regime for the singlet → triplet ISC transition. An inverted Marcus regime being the ideal scenario to introduce polariton assistance has also been found in Ref. 34 for singlet fission and in Ref. 67, the latter in the context of a single-molecule charge transfer process assisted by a photon mode.

On the other hand, for the scenario of a slow vibrational bath, the efficiencies are largely determined by the energetic proximity of the polariton and dark state resonances. The increase in Rabi splitting in this case diminishes the rate of transition of the triplet to the polariton state closer in energy because the vibrational DOS that can assist the transition decreases exponentially with the energy gap between them. Therefore, for the slow phonon bath scenario considered in this work, approach (1) is the only alternative to minimize the polariton-delocalization effect on efficiencies.

## SUPPLEMENTARY MATERIAL

See supplementary material for the appendices which include details on the variational optimization and computation of rate constants.

## ACKNOWLEDGMENTS

L.A.M.-M. is grateful for support of a UC-Mexus CONACyT scholarship for doctoral studies. L.A.M.-M. and J.Y.-Z. acknowledge support of NSF EAGER Grant No. CHE-1836599. S.K.-C. acknowledges support from the Canada Research Chairs program and NSERC Grant No. RGPIN-2014-06129. L.A.M.-M. acknowledges Raphael F. Ribeiro, Jorge Campos-González-Angulo, and Matthew Du for useful discussions.