Optimal power generation using dark states in dimers strongly coupled to their environment

TL;DR

This paper investigates how dark states in strongly vibrationally coupled dimers can enhance power output in light harvesting, revealing differences between homodimers and heterodimers and their robustness to environmental perturbations.

Contribution

It extends dark state analysis to strong vibrational coupling regimes without geometric restrictions, highlighting the stability of heterodimers for power enhancement.

Findings

Homodimers achieve high power but are sensitive to perturbations.

Heterodimers offer intermediate power and are more robust.

Strong vibrational coupling favors heterodimer dark states.

Abstract

Dark state protection has been proposed as a mechanism to increase the power output of light harvesting devices by reducing the rate of radiative recombination. Indeed many theoretical studies have reported increased power outputs in dimer systems which use quantum interference to generate dark states. These models have typically been restricted to particular geometries and to weakly coupled vibrational baths. Here we consider the experimentally-relevant strong vibrational coupling regime with no geometric restrictions on the dimer. We analyze how dark states can be formed in the dimer by numerically minimizing the emission rate of the lowest energy excited eigenstate, and then calculate the power output of the molecules with these dark states. We find that there are two distinct types of dark states depending on whether the monomers form homodimers, where energy splittings and dipole strengths are identical, or heterodimers, where there is some difference. Homodimers, which exploit destructive quantum interference, produce high power outputs but strong phonon couplings and perturbations from ideal geometries are extremely detrimental. Heterodimers, which are closer to the classical picture of a distinct donor and acceptor molecule, produce an intermediate power output that is relatively stable to these changes. The strong vibrational couplings typically found in organic molecules will suppress destructive interference and thus favour the dark-state enhancement offered by heterodimers.

Equations

$$H_{\mathrm{S}}=\sum _{j=1}^2{\delta }_j|j⟩⟨j|+\frac{C}{2}\left(|1⟩⟨2|+|2⟩⟨1|\right),$$

$$\frac{C}{2}=\frac{1}{4\pi r_{12}^3}\left[\text{bi}{d}_1\cdot \text{bi}{d}_2-3\left(\text{bi}{d}_1\cdot \text{bi}{r}_{12}\right)\left(\text{bi}{d}_2\cdot \text{bi}{r}_{12}\right)\right],$$

$$H_{\mathrm{I}}^{\gamma ,j}=\text{rmi}|0⟩⟨j|\sum _{\text{bi}q,\lambda }\text{bi}{d}_j\cdot \text{bi}{u}_{\text{bi}q\lambda }^{\ast }(\text{bi}{r}_j)a_{\text{bi}q\lambda }^{†}+H.c.,$$

$$u_{\text{bi}q\lambda }(\text{bi}r)=\sqrt{\frac{{\nu }_{\text{bi}q}}{2ϵV}}\text{bi}{e}_{\text{bi}q\lambda }{e}^{i\text{bi}q\cdot \text{bi}r},$$

$$H_{\mathrm{I}}^{\mathrm{pn},j}=|j⟩⟨j|\sum _{\text{bi}k}(g_{\text{bi}k,j}{b}_{\text{bi}k,j}^{†}+g_{\text{bi}k,j}^{\ast }{b}_{\text{bi}k,j}),$$

$$H=H_{\mathrm{S}}+\sum _{j=1}^2(H_{\mathrm{I}}^{\gamma ,j}+H_{\mathrm{I}}^{\mathrm{pn},j}+H_{\mathrm{E}}^{\mathrm{pn},j})+H_{\mathrm{E}}^{\gamma }.$$

$$B_j^{\pm }=\exp [\pm \sum _{\text{bi}k}(g_{\text{bi}k,j}{b}_{\text{bi}k,j}^{†}-g_{\text{bi}k,j}^{\ast }{b}_{\text{bi}k,j})/{\omega }_{\text{bi}k,j}]\equiv \prod _{\text{bi}k}{D}_{\text{bi}k,j}\left(\pm \frac{g_{\text{bi}k,j}}{{\omega }_{\text{bi}k,j}}\right),$$

$$H_{\mathrm{SP}}=\sum _{j=1}^2{\delta }_j^{\prime }|j⟩⟨j|+\frac{C^{\prime }}{2}\left(|1⟩⟨2|+|2⟩⟨1|\right),$$

$$J_j(\omega )=\sum _{\text{bi}k}{\left|g_{\text{bi}k,j}\right|}^2\delta (\omega -{\omega }_{\text{bi}k,j}),$$

$${\lambda }_j={\int }_0^{\infty }\text{rmd}\omega \frac{J_j(\omega )}{\omega },$$

$${\kappa }_j\equiv ⟨B_j^{\pm }{⟩}_{\mathrm{pnj}}={\text{rme}}^{-\frac{1}{2}{\varphi }_j(0)},$$

$${\varphi }_j(t)={\int }_0^{\infty }\text{rmd}\omega \frac{J_j(\omega )}{{\omega }^2}\left[\cos (\omega t)\coth \left(\frac{{\beta }_{\mathrm{pn}}\omega }{2}\right)-\text{rmi}\sin (\omega t)\right],$$

$$H_{\mathrm{IP}}^{\gamma }=\sum _{j=1}^2{H}_{\mathrm{IP}}^{\gamma ,j}=\sum _{j=1}^2\text{rmi}|0⟩⟨j|B_j^{-}\sum _{\text{bi}q,\lambda }\text{bi}{d}_j\cdot \text{bi}{u}_{\text{bi}q\lambda }^{\ast }(\text{bi}{r}_j)a_{\text{bi}q\lambda }^{†}+H.c.,$$

$$H_{\mathrm{IP}}^{\mathrm{pn}}=\frac{C}{2}\left(\mathcal{B}|1⟩⟨2|+{\mathcal{B}}^{†}|2⟩⟨1|\right),$$

$${\delta }_{\pm }=\frac{1}{2}\left({\delta }_1^{\prime }+{\delta }_2^{\prime }\pm \eta \right),$$

$$|+⟩=\cos \frac{\chi }{2}|1⟩+\sin \frac{\chi }{2}|2⟩,$$

$$|-⟩=-\sin \frac{\chi }{2}|1⟩+\cos \frac{\chi }{2}|2⟩,$$

$$\cos \chi =\frac{\Delta }{\eta },$$

$$\sin \chi =\frac{C^{\prime }}{\eta }.$$

$$\tilde{H}_{\mathrm{IP}}^{\gamma}=\sum_{\sigma=\pm}\rmi\ket{0}\bra{\sigma}\sum_{\bi{q},\lambda}\big{(}D^{(\sigma)}_{\bi{q}\lambda}(\bi{r}_{1},\bi{r}_{2})\big{)}^{\dagger}a_{\bi{q}\lambda}^{\dagger}+\mathrm{H.c.},$$

$$D_{\text{bi}q\lambda }^{(+)}(\text{bi}{r}_1,\text{bi}{r}_2)=\cos \frac{\chi }{2}\text{bi}{d}_1\cdot \text{bi}{u}_{\text{bi}q\lambda }(\text{bi}{r}_1)B_1^{+}+\sin \frac{\chi }{2}\text{bi}{d}_2\cdot \text{bi}{u}_{\text{bi}q\lambda }(\text{bi}{r}_2)B_2^{+},$$

$$D_{\text{bi}q\lambda }^{(-)}(\text{bi}{r}_1,\text{bi}{r}_2)=-\sin \frac{\chi }{2}\text{bi}{d}_1\cdot \text{bi}{u}_{\text{bi}q\lambda }(\text{bi}{r}_1)B_1^{+}+\cos \frac{\chi }{2}\text{bi}{d}_2\cdot \text{bi}{u}_{\text{bi}q\lambda }(\text{bi}{r}_2)B_2^{+}.$$

$${\tilde{H}}_{\mathrm{IP}}^{\mathrm{pn}}=\sum _{\alpha =x,y,z}{\tau }_{\alpha }{B}_{\alpha },$$

$${\tau }_{\mathrm{x}}=|+⟩⟨-|+|-⟩⟨+|,$$

$${\tau }_{\mathrm{y}}=-\text{rmi}\left(|+⟩⟨-|-|-⟩⟨+|\right),$$

$${\tau }_{\mathrm{z}}=|+⟩⟨+|-|-⟩⟨-|,$$

$$B_{\mathrm{x}}=\frac{C^{\prime }}{2}({\mathcal{B}}^{†}+\mathcal{B})\cos \chi ,$$

$$B_{\mathrm{y}}=\frac{C^{\prime }}{2\text{rmi}}({\mathcal{B}}^{†}-\mathcal{B}),$$

$$B_{\mathrm{z}}=\frac{C^{\prime }}{2}({\mathcal{B}}^{†}+\mathcal{B})\sin \chi .$$

$${\tilde{H}}_{\mathrm{P}}={\tilde{H}}_{\mathrm{SP}}+{\tilde{H}}_{\mathrm{IP}}^{\gamma }+{\tilde{H}}_{\mathrm{IP}}^{\mathrm{pn}}+H_{\mathrm{E}}^{\gamma }+\sum _{j=1}^2{H}_{\mathrm{E}}^{\mathrm{pn},j}.$$

$${\rho }_{\mathrm{E}}^x=\frac{{\text{rme}}^{-{\beta }_x{H}_{\mathrm{E}}^x}}{{\text{tr}}_{\mathrm{E},x}\left({\text{rme}}^{-{\beta }_x{H}_{\mathrm{E}}^x}\right)},$$

$$\frac{\text{rmd}}{\text{rmd}t}{\rho }_{\mathrm{S}}(t)=-\text{rmi}[{\tilde{H}}_{\mathrm{SP}},{\rho }_{\mathrm{S}}(t)]+{\mathcal{D}}^{\gamma }\left({\rho }_{\mathrm{S}}(t)\right)+{\mathcal{D}}^{\mathrm{pn}}\left({\rho }_{\mathrm{S}}(t)\right),$$

$${\mathcal{D}}^x\left({\rho }_{\mathrm{S}}(t)\right)=\sum _{\alpha ,\beta }\sum _{\omega ,{\omega }^{\prime }}{\Xi }_{\alpha \beta }^x(\omega )\left[A_{\beta }^x(\omega ){\rho }_{\mathrm{S}}(t)A_{\alpha }^{x†}({\omega }^{\prime })-A_{\alpha }^{x†}({\omega }^{\prime })A_{\beta }^x(\omega ){\rho }_{\mathrm{S}}(t)\right]$$

$$+H.c.,$$

$${\Xi }_{\alpha \beta }^x(\omega )={\int }_0^{\infty }\text{rmd}t{\text{rme}}^{\text{rmi}\omega t}⟨C_{\alpha }^{x†}(t)C_{\beta }^x(0){⟩}_{\mathrm{E}}.$$

Full text

Keywords: dark state protection, light harvesting, polaron transform, organic solar cell, quantum heat engine

1 Introduction

Organic light harvesting molecules offer the possibility of cheap, stable, flexible and portable photovoltaic devices [1, 2]. However, the efficiency of these devices is much lower than those found in more conventional solar cells, with experimentally reported efficiencies reaching up to 13% after molecular optimization [3], but commonly much lower [1, 4]. One way to improve these figures could be to exploit the results of recent theoretical studies, which show that if quantum interference can be harnessed in coupled organic molecules, then the efficiency of organic materials could be much higher [5, 6, 7, 8].

Pairs of coupled monomers, dimers, are the building blocks of the proposed quantum mechanical light harvesting devices. Modelling has shown that careful optimization of the coupling of the monomers, through tuning their optical properties and relative position and orientation, leads to a hybridization and energy splitting of the single exciton eigenstates. The lower of these states can be optically inactive – the so-called dark state [5, 6, 9, 10], whereas the higher lying state has an enhanced transition dipole – this is the bright state. The dark state is accessed by non-radiative, vibrational relaxation from the bright state. So long as the energy splitting of the bright and dark states is large enough to suppress vibronic re-excitation, then the absorbed energy becomes trapped as exciton recombination is suppressed. This process breaks the detailed balance that leads to the Shockley-Queisser limit, and has been predicted to increase the power output and efficiency of a device built from such molecules [5, 6, 7, 9, 11, 12, 13, 14, 15].

Dark state protection has been shown to work effectively for specific geometries of both homodimers [5] and heterodimers [6], and more recently in linear chains of coupled homodimers [16]. However, these studies are restricted to the case of weak environmental coupling; the picture is more complicated in typical organic molecules and dye pigments, since these have strong coupling to their vibrational environments [17, 18, 19, 20]. In this paper, we therefore go beyond these earlier works by treating strong environmental coupling in unrestricted geometries, and including realistic inter-monomer dipole-dipole couplings. We treat the strong coupling by using a polaron transformation [21]. In the polaron frame, the monomer energy splittings and inter-monomer coupling become renormalized such that the vibronic environment is in equilibrium with the excited state, rather than displaced by it. Energy transfer in the polaron frame has been studied before [21, 22, 23], but not in the context of dark state protection. We will show that strong coupling has a profound effect on quantum interference processes such that it no longer necessarily improves the efficiency of devices. Rather we will propose a different strategy for designing energy harvesting dimer molecules, in which a dark state is chiefly localized on the lower energy monomer.

The paper is structured as follows: Section 2 begins by introducing the dimer model, before showing how the Hamiltonian is transformed in the polaron frame. In Section 3 we discuss the Born-Markov master equation formalism before, in Section 4, we derive the emission rate from the dark state and explore the conditions under which dark state protection is optimized. We identify different mechanisms for dark state formation for homogeneous and heterogeneous dimers. In Section 5 we show absorption and emission spectra for the dimers under conditions in which a dark state is expected to form and compare this to conditions which prevent dark state formation. In doing so we provide a way to determine experimentally if a dark state can form. In Section 6 we explain the theoretical model of power extraction and calculate the power output of dimers over a wide range of realistic parameters. We find that the optimal power output occurs when the dark state is formed in either the homogeneous or heterogeneous dimer and discuss in which regimes either is preferable. Finally, we conclude in Section 7.

2 Model

The dimer consists of two monomers with single dipole moments, $\bi{d}_{j}$ located at positions $\bi{r}_{j}$, where $j=1,2$ labels the dipoles. This system is illustrated in a cartoon in Figure 1(a) where the dipoles exist as part of a larger protein structure. By virtue of the dipole-dipole coupling, the monomers form symmetric and antisymmetric eigenstates (denoted $\ket{+}$ and $\ket{-}$ respectively), depicted as green clouds. In an ideal scenario quantum interference is constructive at the $\ket{+}$, leading to enhanced photon absorption, and is conversely suppressed at the $\ket{-}$, forming a dark state. After absorption at the brighter $\ket{+}$, the excitation is transferred non-radiatively to the darker $\ket{-}$ where it becomes trapped until extraction to an idealized load to produce power. We will now describe the dimer system mathematically.

2.1 Dimer Hamiltonian

We assume that each monomer can be treated as a two level system. The excited states are denoted $\ket{1}$ and $\ket{2}$, and the ground state as $\ket{0}$. This lab frame energy level diagram is shown in Figure 1(c). We use natural units throughout ($\hbar=c\equiv 1$). The Hamiltonian describing this system part is

[TABLE]

where $\delta_{j}$ is the energy splitting of monomer-$j$ and $C$ is the static dipole-dipole coupling:

[TABLE]

where $\bi{r}_{12}=\bi{r}_{1}-\bi{r}_{2}$ is the separation vector of the monomers. No restrictions are placed on any of the system parameters. We only consider the single excitation subspace since under illumination with solar temperature radiation the dimer will very rarely hold two excitations.

Each monomer couples to a global, multimode photon environment described by $H_{\mathrm{E}}^{\gamma}=\sum_{\bi{q},\lambda}\nu_{\bi{q}}a^{\dagger}_{\bi{q},\lambda}a_{\bi{q},\lambda}$ where $\nu_{\bi{q}}$ is the frequency of mode $\bi{q}$ and $a^{\dagger}_{\bi{q},\lambda}$ ($a_{\bi{q},\lambda}$) is the creation (annihilation) operator for a photon of mode $\bi{q}$ and polarization $\lambda$. Under the rotating wave approximation, the interaction of monomer-$j$ with this field results in the Hamiltonian term,

[TABLE]

where H.c. denotes the Hermitian conjugate.111We have directly verified that neglecting counter rotating terms in (3), and ignoring the double excited state, makes negligible difference to any results we present. We take the spatial mode functions $u_{\bi{q}\lambda}(\bi{r})$ of the field to be those of free space with volume $V$ and permittivity $\epsilon$,

[TABLE]

where $\bi{e}_{\bi{q}\lambda}$ are unit vectors describing the polarization state.

Finally, the monomers interact with their own local phonon environment $H_{\mathrm{E}}^{\mathrm{pn},j}=\sum_{\bi{k}}\omega_{\bi{k},j}b_{\bi{k},j}^{\dagger}b_{\bi{k},j}$ where $\omega_{\bi{k},j}$ is the frequency of a phonon of mode $\bi{k}$ within monomer-$j$ and $b^{\dagger}_{\bi{k},j}$ ($b_{\bi{k},j}$) is the creation (annihilation) operator for that phonon. These interactions are represented by the usual displacement of the excited states,

[TABLE]

with coupling strength $g_{\bi{k},j}$. The full Hamiltonian is then

[TABLE]

2.2 Polaron transformation

The presence of large Stokes’ shifts in candidate monomers for implementing dark state protection ideas, such as those found in [6], indicates that these molecules have strong coupling to vibrational modes. Therefore, we transform to the polaron frame which takes some of the phonon interaction Hamiltonian into the system, leaving a residual interaction that can be treated within a weak coupling theory.

The polaron transformation is generated by the unitary operator $U_{\mathrm{P}}=\rme^{G}$ where $G=\sum_{\bi{k},j}\ket{j}\bra{j}(g_{\bi{k},j}b^{\dagger}_{\bi{k},j}-g^{*}_{\bi{k},j}b_{\bi{k},j})/\omega_{\bi{k},j}$. This can be decomposed into the dipole basis as $\rme^{\pm G}=\ket{0}\bra{0}+\sum_{j}B^{\pm}_{j}\ket{j}\bra{j}$, where

[TABLE]

is a product of displacement operators, themselves defined by $D_{\mathrm{x}}(\alpha)=\rme^{\alpha b^{\dagger}_{\mathrm{x}}-\alpha^{*}b_{\mathrm{x}}}$ [21, 22, 23].

After applying the polaron transformation to the full Hamiltonian (6), we subsequently partition the terms into system, environment and interaction Hamiltonians in the usual way [21, 23, 24]. We find that both the photon and phonon environment Hamiltonians are unchanged. The system Hamiltonian becomes

[TABLE]

where the primed variables indicate renormalization by phonon interactions. Specifically, defining the spectral density of the phonon environment coupled to monomer-$j$ as

[TABLE]

we have $\delta_{j}^{\prime}=\delta_{j}-\lambda_{j}$ where

[TABLE]

is the reorganization energy of the phonon environment, and $C^{\prime}=\kappa_{1}\kappa_{2}C$ where

[TABLE]

and

[TABLE]

is the phonon propagator for the phonon environment with inverse temperature $\beta_{\mathrm{pn}}~{}=~{}\left(k_{\mathrm{B}}T_{\mathrm{pn}}\right)^{-1}$. The polaron-frame system Hamiltonian is depicted in Figure 1(c). In the polaron frame, the photon and phonon interaction Hamiltonians become

[TABLE]

where $\mathcal{B}=B_{1}^{+}B_{2}^{-}-\kappa_{1}\kappa_{2}$, and note that $(B_{j}^{+})^{\dagger}=B_{j}^{-}$.

2.3 Diagonalization

The stronger the dipole-dipole coupling $C^{\prime}$ is compared to the detuning of the monomers $\Delta$ in the polaron frame, the more delocalized the excitons are over the monomers [25]. The phase relationship between the monomer components of the eigenstate wave functions then determines their optical dipole matrix element with the ground state, and fine tuning the degree of delocalization leads to the formation of dark states for various monomers [6].

The diagonalized system Hamiltonian (8) is written $\tilde{H}_{\mathrm{SP}}=\sum_{\sigma=\pm}\delta_{\sigma}\ket{\sigma}\bra{\sigma}$, where the eigenvalues are

[TABLE]

with eigenstate detuning $\eta=\sqrt{\Delta^{2}+C^{\prime 2}}$, and renormalized monomer detuning $\Delta~{}=~{}\delta_{1}^{\prime}-\delta_{2}^{\prime}$. The symmetric and antisymmetric eigenstates (see Figure 1(c) for a depiction) are

[TABLE]

where

[TABLE]

Although we have named the $\ket{+}$ and the $\ket{-}$ the symmetric and antisymmetric eigenstates respectively, this definition is strictly only correct for positive couplings, $C$. When the coupling becomes negative, the symmetry of the eigenstates swap. In the eigenbasis, the photon interaction Hamiltonian becomes

[TABLE]

where

[TABLE]

Finally, the phonon Hamiltonian becomes

[TABLE]

where we have introduced the Pauli operators

[TABLE]

and the following groupings of phonon operators

[TABLE]

Thus, the full Hamiltonian in the polaron frame reads:

[TABLE]

3 Master equation formalism

To calculate the energy transfer dynamics we use the standard second order Born-Markov master equation in the polaron frame. A full discussion on the master equation formalism can be found in [26]. The Born-Markov assumptions are that the density matrix $\rho(t)$ is separable at all times into a system and environment part, and that the environment retains no memory of its interactions with the system. Therefore we can write the total density matrix in the factorized form $\rho(t)=\rho_{\mathrm{S}}(t)\otimes\rho_{\mathrm{E}}^{\mathrm{pn}}\otimes\rho_{\mathrm{E}}^{\gamma}$ where $\rho_{\mathrm{S}}(t)$ is the system density matrix, and $\rho_{\mathrm{E}}^{x}$ are the density matrices for boson environments $x$, and are constant in time and in thermal equilibrium, i.e.

[TABLE]

where $\beta_{x}$ ($\tr_{\mathrm{E},x}$) is the inverse temperature of (trace over) environment $x$. We are interested in the dynamics of the system, and because of the density matrix partitioning we are able to trace out the environment degrees of freedom and derive a dynamical equation $\rho_{\mathrm{S}}(t)=\tr_{\mathrm{E}}\left[\rho(t)\right]$ where E denotes both environments.

The closed evolution of the system is then treated exactly, but the effects of the system-environment interactions on the dynamics are treated perturbatively to second order within the Born-Markov approximation. To summarize a standard result, one finds that in the Schrödinger picture,

[TABLE]

where $\mathcal{D}^{x}(\rho)$ are the dissipators for environment $x$. After writing the interaction Hamiltonians as $\tilde{H}_{\mathrm{IP}}^{x}=\sum_{\alpha}A_{\alpha}^{x}\otimes C_{\alpha}^{x}$ where the $A_{\alpha}^{x}$ and $C_{\alpha}^{x}$ operators are in the system and environment(s) Hilbert spaces respectively, then in the Schrödinger picture (and polaron frame) the dissipators have the form

[TABLE]

where the rates of each process are given by the environment correlation functions (ECFs)

[TABLE]

In (16qsuvz), the system operators are in eigenoperator form, defined as $A_{\alpha}^{x}(\omega)=\sum_{\epsilon}\mathcal{P}(\epsilon)A_{\alpha}^{x}\mathcal{P}(\epsilon+\omega)$ where $\mathcal{P}(\epsilon)=\ket{\epsilon}\bra{\epsilon}$, $\tilde{H}_{\mathrm{SP}}\ket{\epsilon}=\epsilon\ket{\epsilon}$ and the sum extends over all eigenvalues $\epsilon$. Similarly, the environment operators in (16qsuvaa) are evaluated in the interaction picture.

We summarize the contributions to the master equation from each term of (16qsuvy) in B. In the following we only explicitly discuss terms which influence the formation of the dark state.

4 Creating dark states

For an eigenstate to be dark, the emission of a photon must be suppressed. It is also necessary that the dominant role of phonon processes is to cause excitations to be transferred into the dark state from the bright state, rather than vice-versa. The latter is easily achieved by choosing our target dark state as the lower energy, antisymmetric eigenstate, $\ket{-}$. The former is only achieved by carefully tuning all parameters of the system, and the phonon coupling strength, such that the photon emission rate is suppressed from the $\ket{-}$.

4.1 Antisymmetric eigenstate photon rates.

We derive the photon emission and absorption rates in A but here we quote the results. We find that the rates between the $\ket{-}$ and the ground state can be written as

[TABLE]

where $\mu=\mathrm{A},\mathrm{E}$ for absorption and emission, and the angle $\chi$ is defined by (16qa) and (16qb). In (16qsuvab), $\Gamma_{j}^{\mu}(\omega)$ are the rates for the phonon-renormalized electronic transitions of monomer-$j$ and $\Gamma_{12}^{\mu}(\omega)$ are the contributions to the rates due to the collective effects of the monomers in the dimer. Therefore, recalling the definition of the $\ket{-}$ (16b), the first two terms in (16qsuvab) are seen to describe the contribution from each monomer independently. These include a trigonometric factor describing the localization on either monomer, as well as the electronic transition rate of the associated monomer. The final term describes the effect of the coupling between the monomers and the overall sign determines if this leads to constructive ($>0$) or destructive ($<0$) interference. A term due to destructive interference in the rate expression is beneficial for creating a dark state.

The terms due to the monomer coupling, $\Gamma^{\mu}_{12}(\omega)$ can be derived analytically (see A.2) and for absorption and emission are found to be

[TABLE]

where

[TABLE]

is the bare electronic transition rate, $N(\omega)=\left(\rme^{\beta_{\gamma}\omega}-1\right)^{-1}$ is the photon Bose-Einstein distribution and $\kappa_{j}$ are defined by (11). The cross function is

[TABLE]

with

[TABLE]

where the hat denotes the unit vector. The range of the cross function is $-{1\leq\mathcal{F}(\omega r_{12})\leq 1}$ and its value indicates the effect of the delocalization on the photon rates (16qsuvab). In realistic dimer systems, the eigenstate energy splittings $\delta_{\pm}\sim\mathrm{eV}$ and monomer separation $r_{12}\sim\mathrm{nm}$, in which case we can use the following limit

[TABLE]

We can immediately see that for the collective effects to cause destructive interference and suppress the $\ket{-}$ rates we require $\mathrm{sign}\left[C\right]=\mathrm{sign}\left[\mathcal{F}\right]$. If this is not true, then one instead finds that the interference of the $\ket{-}$ rates is constructive. The magnitude of the interference term depends on: the size of the renormalized monomer detuning $\Delta$ compared to the renormalized dipole-dipole coupling $C^{\prime}$; the geometry of the dimer though the cross function $\mathcal{F}$ and coupling $C^{\prime}$; the bare monomer rates $\gamma_{j}(\omega)$ and the phonon coupling strengths through the factor $\kappa_{1}\kappa_{2}$ and renormalization of the monomer transition rates.

The monomer electronic transition rates, $\Gamma_{j}^{\mu}(\omega)$ in (16qsuvab) are discussed in A.3. These cannot be derived analytically and the usual approximation would be to assume that the photon spectral density is flat for frequencies near to the eigenfrequencies of the system [21, 24]. This is the flat spectral density approximation (FSDA). If we used this approximation we would find that $\Gamma_{j}^{\mathrm{A}}(\omega)\rightarrow\Gamma_{j}^{\mathrm{A},\mathrm{FSDA}}(\omega)=\gamma_{j}(\omega)N(\omega)$ and $\Gamma_{j}^{\mathrm{E}}(\omega)\rightarrow\Gamma_{j}^{\mathrm{E},\mathrm{FSDA}}(\omega)=\gamma_{j}(\omega)\left[1+N(\omega)\right]$. Hence, under the FSDA the phonon coupling does not affect the electronic transition rates. It is therefore not surprising that the rates derived using the FSDA are equivalent to those derived in the limit of weak phonon coupling. For strong phonon couplings this approximation quickly breaks down [27], and we show this explicitly in Figure 6 in A.3. Instead of making the FSDA, and because the effect of increasing the phonon coupling is to increase the lifetime of the monomers, we write

[TABLE]

where the $\zeta^{\mu}_{j}(\omega)$ have values $0\leq\zeta_{j}^{\mu}(\omega)\leq 1$, tending towards zero for increasing phonon coupling strength. Deviations from unity measure the error in assuming that the photon spectral density can be regarded as being flat, i.e. in assuming that there is no influence on optical emission rates by the vibrational environments. A similar result has been recently found in the context of mapping the vibrational environment onto a collective coordinate [27].

The functional form of $\zeta^{\mu}_{j}(\omega)$ is dependent on the choice of phonon spectral density, $J(\omega)$. An exact solution can only be derived for simple spectral densities. However, an approximate but analytic solution can be found for many spectral densities. We give this approximate expression (16qsuvacaiatazbj) in A.3. In Figure 6 of A.3, we plot the approximated and exact expressions for $\zeta^{\mu}_{j}(\omega)$ for a simple spectral density (where the exact solution exists) for a range of phonon couplings to illustrate that the approximate solution works well in this case. Throughout this paper, we use the spectral density

[TABLE]

for which there is no exact expression of $\zeta^{\mu}_{j}(\omega)$. Therefore, we assume the approximated expression (16qsuvacaiatazbj) is accurate for this spectral density.

We can finally summarize the $\ket{-}$ absorption and emission rates by writing

[TABLE]

where the rate coefficients are

[TABLE]

for $\mu=\mathrm{A},\mathrm{E}$.

4.2 The dark state

To make the $\ket{-}$ a dark state, we need to minimize the rate coefficient (16qsuvacaiam). To do so it is useful to define the dipole magnitude ratio, $z\equiv\left|\bi{d}_{2}/\bi{d}_{1}\right|$ which, under our definitions of monomer-1 and monomer-2, can take the values $0<z\leq 1$. We identify two parameter regimes where the emission rate is minimized by two distinct mechanisms.

The first mechanism occurs when the monomer detuning is much greater than the dipole-dipole coupling ($\Delta\gg C^{\prime}$). From (16b) we see that this causes the antisymmetric eigenstate to localize onto monomer-2 ($\ket{-}\approx\ket{2}$) which is reflected in the emission and absorption rates of the $\ket{-}$ because

[TABLE]

and these are the respective rates for monomer-2. For the large detuning and small dipole magnitude ratio necessary to achieve localization, the emission rate of monomer-2 is small because $\gamma_{2}(\delta_{2}^{\prime})\propto z^{2}\left(\delta_{1}^{\prime}-\Delta\right)^{3}$ (see (16qsuvacad)). Therefore, very heterogeneous dimers can lead to a dark $\ket{-}$. This mechanism does not make use of quantum effects, but simply localizes the eigenstate onto an already relatively optically inactive monomer. Therefore, these heterodimers cannot create perfect dark states without completely decoupling monomer-2 from the photon field. This is close to the classical picture of a dimer consisting of a donor (monomer-1) and acceptor (monomer-2) molecule.

The second mechanism maximizes the destructive interference term of (16qsuvacaiam). This therefore requires $C^{\prime}\gg\Delta$ (so that $\sin\chi\approx 1$) which is achieved for homodimers because these have $(z,\Delta)\rightarrow(1,0)$. This can be more useful than the first mechanism because we can produce a highly (and in some ideal situations, perfectly) decoupled eigenstate but still have large dipole-dipole couplings. Since the phonon transition rates between the $\ket{+}$ and the $\ket{-}$ scale as $\sim C^{\prime 2}$, this also enables easier access to the dark state. The drawback of this mechanism is that the interference is easily suppressed by strong phonon couplings and sub-optimal monomer orientations; this is in contrast to dark states in heterodimers, which are largely unaffected. This can can be inferred from the interference term in (16qsuvacaiam), which is proportional to $\kappa_{1}\kappa_{2}$ and $\mathcal{F}$, and so away from ideal geometries ($\mathcal{F}=1$) and very weak phonon couplings ($\kappa_{j}\approx 1$) the destructive interference is reduced.

4.3 Minimization of the emission rate

In Figure 2 we explore the dependence of the $\ket{-}$ emission rate on detuning, orientation and phonon coupling strength using the complete expression for the dipole-dipole coupling given by (2). We plot the logarithm of the $\ket{-}$ emission rate, $\Gamma_{-}^{\mathrm{E}}$, divided by the emission rate of monomer-2, $\Gamma_{2}^{\mathrm{E}}=\zeta_{2}^{\mathrm{E}}(\delta_{2}^{\prime})\gamma_{2}(\delta_{2}^{\prime})$, where $\delta_{2}^{\prime}=\delta_{1}^{\prime}-\Delta$ and we fix $\delta_{1}^{\prime}=2.8$ eV. We are therefore using the lower energy monomer as a benchmark to test the decoupling. The coordinate system of the dipoles is defined in Figure 1(b).

If the emission rate of the $\ket{-}$ is larger than that of monomer-2 then this does not necessarily mean that the power output of the dimer will also be lower. This is explored in Section 6. In Figure 2, we fix $\phi_{1}=\frac{\pi}{2}$ in all subfigures so that for all orientations $\mathrm{sign}\left[C\right]=\mathrm{sign}\left[\mathcal{F}\right]$ and therefore the interference is alwayy destructive. Additionally, we choose the value of $z$ such that the $\ket{-}$ emission rate (16qsuvacaiam) is minimized in the ideal orientation - i.e. in each subfigure we use the same value of $z$ for all orientations. The ideal orientation for maximizing destructive interference is dipoles parallel or anti-parallel to each other, and perpendicular to the interconnecting vector so that $\mathcal{F}=1$ (i.e. $(\phi_{1},\phi_{2},\theta_{12})=(\frac{\pi}{2},\frac{\pi}{2},n\pi)$ for $n=0,1,2$). Figure 2 then allows us to determine the effect of imperfect orientation.

We keep the renormalized coupling fixed at $C^{\prime}=\left(100\ \mathrm{meV}\right)z$ and study a wide range of detunings. Therefore, we explore both the homodimer parameter regime where destructive interference is important ($C^{\prime}>\Delta$), as well as the heterodimer parameter regime where instead localization is important ($C^{\prime}\lesssim\Delta$). For parameters where the interference term is strongly suppressed the only minimum in the emission rate occurs at $z=0$, even in the optimal orientation. In Figure 2, these situations are indicated with an asterisk ($\ast$) and we instead choose a $z$ value by using

[TABLE]

derived in C. This expression gives the value of $z$ for which the destructive interference minimizes the emission rate (16qsuvacaiam) in the simplified case where $C^{\prime}$ does not depend on dipole strength, an assumption made in [6].

Figure 2 clearly illustrates that the destructive interference is larger for detunings much smaller than the renormalized coupling, as evident in the bottom two rows. Additionally, we can see the extreme sensitivity of the emission rate to orientation and phonon coupling strength: As a guide, reorganization energies in organic molecules are typically of the order $100~{}\mathrm{meV}$. We also see in Figure 2 that to form dark states at strong phonon couplings it is more beneficial to have heterodimers ($\Delta>C^{\prime}$). As expected in these cases, the $\ket{-}$ emission is equal to the benchmark emission, which would be small due to the reduced $z$ and $\delta_{2}^{\prime}$ values. At strong coupling in homodimers ($\Delta=0$) we see that the emission rate of the $\ket{-}$ can exceed that of monomer-2. This is because with suppressed destructive interference, the partial localization of the eigenstate over the more optically active monomer-1 dominates. In fact, at very strong couplings and small detunings it is preferable to orientate the dipoles such that the coupling $C^{\prime}$ is reduced much smaller than $\Delta$ to localize onto the lower energy monomer, so forming a heterodimer. This is not possible in the bottom row where the detuning is zero, except for $\theta_{12}=\frac{\pi}{2},\frac{3\pi}{2}$ for which $C^{\prime}=0$.

5 Spectra

Producing homodimers that maximize destructive interference requires precise control of the monomer orientations. This is especially true at strong phonon couplings where the interference has already been suppressed. To highlight experimentally detectable signatures of dark state formation in a homodimer, we plot theoretically calculated absorption and emission spectra in Figure 3. In D.1 we plot spectra for the heterodimer.

We choose a realistic molecular phonon coupling strength for both environments, $\lambda_{j}=100$ meV for $j=1,2$ and keep the phonon renormalized dipole-dipole coupling at $C^{\prime}/z=100$ meV in the optimum orientation. In each row of Figure 3 we choose different orientations of the dipoles. The lowest row corresponds to the orientation which optimizes dark state formation (signified by $\mathcal{F}=1$) and progressively higher rows correspond to less desirable geometries (signified by $\mathcal{F}<1$). The specific angles are detailed in the caption. We use the same dipole magnitude ratio, $z=0.735$, in each row. This value minimizes the $\ket{-}$ emission rate, i.e. creates the dark state in the homodimer, for the optimum orientation. In this way, the lowest row corresponds to the dimer one should aim to produce in order to maximize decoupling from the photon field, and the higher two rows show realistic experimental deviations from this perfect scenario. We note that the minimizing value of $z$ is not 1 (i.e. the dimer is not a completely homogeneous) and this is because of the strong phonon coupling.

In each subfigure of Figure 3 we plot the spectra with and without the phonon sideband (as solid blue and dashed green lines respectively). Using these plots we are also able to calculate the fraction of emission or absorption which occurs through the phonon sideband numerically, which is the parameter $f^{\mathrm{pn}}$: The specifics of this calculation, and how we find the spectra, are discussed in detail in D.2.

From Figure 3 we can identify two different signatures that a dark state has formed in the homodimer. Firstly, the darker the $\ket{-}$, the more reliant the dimer is on phonons for emission. This is reflected by the trend in $f^{\mathrm{pn}}$. For these strong phonon and dipole-dipole couplings, excitations are mostly transferred to the $\ket{-}$ before they are immediately re-radiated from the $\ket{+}$. Therefore, emission can only occur after phonon transfer out of the dark state, which becomes increasingly more true the darker the $\ket{-}$ is. Crucially, the excitations are trapped without phonon transfer. This trend is not seen in the absorption spectra because the dimer can still be excited through the $\ket{+}$ even if the $\ket{-}$ is dark. In this case, there is no absolute blocking of absorption in the absence of phonon processes. As we show in D.3, one can derive expressions for $f^{\mathrm{pn}}$ which are exact for absorption and approximate for emission. We also demonstrate the accuracy of these expressions in Figure 8. For the homodimer, one finds that (for $\kappa_{1}=\kappa_{2}\equiv\kappa$) the fractional absorption and emission in the phonon sideband are

[TABLE]

For the parameters used in Figure 3 the value of $1-\kappa^{2}=0.33$, in agreement with the numerical values calculated from the spectra. Indeed, (16qsuvacaiaq) allows one to determine $\mathcal{F}$ experimentally for a homodimer by measuring $f^{\mathrm{pn}}_{\mathrm{E}}$.

The second parameter which helps to identify dark state formation in a homodimer is the total amount of light that is absorbed and emitted by the dark state, $P_{\mathrm{d}}$. As expected, the total absorption and emission from the dark state dramatically increases as the destructive interference weakens due to poor geometry.

6 Power output

Finally, we have calculated the power output of the dimer, shown in Figure 4; we use six panels for different values of dipole coupling $C^{\prime}$ and phonon coupling $\lambda$, and within each panel vary $z$ and $\Delta$. To calculate the power output we formally introduce the idealized load model alluded to in Figure 1(c). Following [12, 13], we measure the power output of the dimer by coupling it to a distinct two level system called a trap, which serves as the idealized load. The extraction process is modelled by adding an incoherent process which moves excitations from a chosen eigenstate in the dimer (which we will choose to be the $\ket{-}$) to the ground state, whilst simultaneously exciting the trap, all at a rate $\gamma_{\mathrm{x}}$. There is also an incoherent process which causes transfer from the excited trap-state $\ket{\alpha}$ to the ground trap-state $\ket{\beta}$ at rate $\gamma_{\mathrm{t}}$, and this is always optimized to give maximal power output. The voltage and current through the trap at inverse temperature $\beta_{\mathrm{t}}=\left(k_{B}T_{\mathrm{t}}\right)^{-1}$ are then defined as

[TABLE]

where $\delta_{\mathrm{t}}$ is the trap energy splitting, $e$ is the electron charge and the $P_{i}$ for $i=\alpha,\beta$ are the trap populations in the steady state. The trap Hamiltonian is $H_{\mathrm{t}}=\delta_{\alpha}\ket{\alpha}\bra{\alpha}+\delta_{\beta}\ket{\beta}\bra{\beta}$, where $\delta_{\mathrm{t}}=\delta_{\alpha}-\delta_{\beta}$ is always set equal to the energy splitting of the state from which the extraction is occurring. The extraction and transfer processes are added phenomenologically with Lindblad operators [12]. For a full description see E. Other models for the extraction process and power generation are also possible [28].

Instead of plotting absolute power (which is shown in F), Figure 4 shows the power output of the dimer compared to a benchmark. The benchmark is the total power that would be extracted from the monomers in the dimer if they were completely independent and coupled to separate traps. In the benchmark, the extraction rate from either monomer is equal to the rate that we extract from the dark state in the dimer. Therefore, if the dimer produces more power than the benchmark, then we can say with certainty that using the dimer conveys a definite advantage because overall the dimer has only half the extraction rate of the benchmark. In Figure 4, the power output is plotted as a function of dipole magnitude ratio $z$ and detuning $\Delta$, over dipole-dipole couplings and phonon couplings ranging from the weak to strong regimes. The extraction rate is fixed at $\gamma_{x}=100$ neV. Indicated on the plots are four power maxima which correspond to forming dark states in the homogeneous and heterogeneous dimers.

The four power maxima are labelled ‘homodimer maximum-X’ (HOX) and ‘heterodimer maximum-X’ (HEX) for $\mathrm{X}=1,2$ which label the strong and weak dipole-dipole coupling regimes respectively. The label of the maximum refers to the detuning and relative dipole strengths of the monomers in the dimer and therefore, as discussed in Section 4, the mechanism by which the $\ket{-}$ is optically darker. We see that the HO1 has a higher power output than the HO2 because of stronger dipole-dipole coupling which increases the destructive interference and that, as expected, strong phonon couplings quench this. Since the monomer detunings are zero at the HOX, the eigenstate detuning $\eta=\sqrt{\Delta^{2}+C^{\prime 2}}=C^{\prime}$. To maximize overall phonon transfer from the $\ket{+}$ to the $\ket{-}$ the eigenstate detuning must optimized. The detuning must be sufficiently large to suppress excitation via phonon transfer back out of the dark state, but smaller than the cut-off frequency of the vibrational baths $\omega_{c}$ so that the modes that are most strongly coupled are used to complete the transfer. As before, we use a cut-off frequency found in realistic molecules $\omega_{c}=0.3$ eV [19] therefore at the HOX, $\eta<\omega_{c}$. Additionally, in the polaron frame the phonon rates scale with $C^{\prime 2}$ (see B.2). Therefore, in addition to being darker, phonon transfer is also much larger at the HO1 than the HO2.

As discussed in Section 4, at the HEX the $\ket{-}$ is localized onto monomer-2, which requires small $z$ and $\delta_{2}^{\prime}$. With these parameters monomer-2 is also relatively optically inactive. This localization also implies that $C^{\prime}\ll\Delta$ and therefore the eigenstate splitting is dominated by the monomer detuning, i.e. $~{}{\eta=\sqrt{\Delta^{2}+C^{\prime 2}}\approx\Delta}$. The sacrifice made with the small dipole strength means that phonon transfers, which scale with $C^{\prime 2}$ are suppressed. Therefore, the power maxima occur at the eigenstate detuning that will maximize phonon transfer from the $\ket{+}$ to the $\ket{-}$. For our chosen cut-off frequency, phonon transfers maximize for eigenstate detunings $\approx 0.2$ eV which explains the location of the HEX. We show important rates in understanding the power outputs in G.

In Figure 5 we plot both the benchmarked power and the absolute power at each of these maxima as a function of extraction rate. Although absolute power output is always greater for faster extraction to the trap, it does not necessarily mean that a real device would operate in these high extraction regimes. This is discussed in the supplementary material of [12]. There it is argued that if one assumes a simple linear relationship between the total extraction rate and the average number of traps adjacent to any light harvesting dimer. Natural photosynthetic systems commonly contain a few hundred antennae for each reaction centre [29], which leads to a relatively small transfer rate. It can therefore be important to improve power output at low extraction rate, which is the regime in which dark state protection typically works well.

From Figure 4 and Figure 5 we identify that the strengths of the dipole-dipole coupling and phonon couplings should dictate whether the target dimer should be homogeneous or heterogeneous. The power output achieved with homodimers can be large, however, it requires either extremely large dipole-dipole couplings ($C^{\prime}>100$ meV) or extremely weak phonon couplings $\lambda\approx 0.1$ meV. Such strong dipole-dipole couplings are only achievable with unrealistically short monomer lifetimes and/or separations. Organic monomers typically have reorganization energies close to $100$ meV and in this regime aiming for the HOX is not ideal. On the other hand, the HE1 produces a high power output which is relatively unaffected by increasing phonon coupling. This reflects the fact that the HE1 is not reliant on an ever-increasingly quenched destructive interference for low $\ket{-}$ emission. In the HE2, increasing phonon coupling actually increases power output because, at this small dipole-dipole coupling, the power bottleneck is the phonon transfer rate from the bright to the dark state.

In Figure 5 we see that for large extraction rate the advantage conferred by dark state protection is lost, for all of the HEX and HOX. This is because extraction must, beyond some threshold rate, always occur much more quickly then any exciton recombination. There is, however, also a low extraction rate regime where the power output of the benchmark slightly exceeds that of the HO2 dimer. This happens for strong phonon coupling ($\lambda=100$ meV) which destroys the destructive interference effect for this dark state. Since there is then no dark state advantage, and the benchmark has overall twice the extraction, the benchmark has a higher power output than the dimer.

7 Conclusion

Dimers are the relatively simple molecules that can form the building blocks of larger more sophisticated light harvesting complexes that could form the components of organic solar cells [5]. Therefore, a detailed understanding of the energy transfer processes in these systems subject to realistic environmental constraints is important. Here we have developed a theory for energy transfer in a dimer which significantly improves upon previous models. More importantly, we have tailored the model to better suit biologically inspired dimers where the coupling to vibrational modes can be strong and a polaron theory is needed to capture the details of the dynamics.

We have found that there can be an enhanced power output for the homodimers which exploit the quantum effect of delocalization, resulting in destructive interference of the optical emission rate. However, this is extremely sensitive to orientation, phonon coupling and the ratio of detuning with dipole-dipole coupling. For strong phonon coupling, which is close to experimental reality, it can be beneficial to instead design a heterodimer that localizes the eigenstates onto the monomers. The optimal dimer choice also depends on how quickly excitons are extracted into the electrical circuit and there are additional, distinct challenges in constructing the best homogeneous or heterogeneous dimer.

Engineering an heterodimer (HEX) requires careful pairing of monomers to localize the $\ket{-}$ onto a relatively dark monomer-2, whilst ensuring that there is sufficient phonon transfer into the $\ket{-}$ so that it is populated. However, this is mitigated by there being a vast number of monomers to choose from, as shown in the supplementary material of [6]. On the other hand, creating a homodimer (HOX) simply requires the use of identical monomers with weak vibrational environments. However, the dipoles of the monomers must be orientated with relatively high precision to maximize destructive interference, and we have identified experimentally detectable signatures of when this has been achieved. Since destructive interference is of little importance to heterodimers, their orientation can instead be utilized to fine-tune the magnitude of the dipole-dipole coupling. This enables one to find the optimal balance between localization and phonon transfer. In previous literature, there has been focus on producing dark states by maximizing the destructive interference using homodimers. Here we show that in organic solar cells, where the vibrational coupling is likely to be strong, this interference is significantly quenched, and it is instead better to produce heterodimers where the eigenstate localizes onto a relatively optically dark monomer.

Acknowledgements

We thank Jonathan Keeling, Alexandre Coates, William Brown, Aidan Strathearn and Dale Scerri for useful discussions. D.R. acknowledges studentship funding from EPSRC under grant no. EP/L015110/1. E.M.G. thanks the Royal Society of Edinburgh and the Scottish Government for support.

Appendices

Appendix A Emission and absorption rates from the antisymmetric eigenstate

To derive the absorption and emission rates of the $\ket{-}$ we start from the ECF definition (16qsuvaa). We first partition the photon interaction Hamiltonian (16qr) into system and environment operators as $\tilde{H}_{\mathrm{IP}}^{\gamma}=\sum_{\alpha}A^{\gamma}_{\alpha}\otimes C_{\alpha}^{\gamma}$ where $A_{\mathrm{a}}^{\gamma}=\ket{0}\bra{+}$, $A_{\mathrm{b}}^{\gamma}=\ket{0}\bra{-}$, $A_{\mathrm{c}}^{\gamma}=A_{\mathrm{a}}^{\gamma\dagger}$ and $A_{\mathrm{d}}^{\gamma}=A_{\mathrm{b}}^{\gamma\dagger}$. The corresponding environment operators are $C_{\mathrm{a}}^{\gamma}=\rmi\sum_{\bi{q}\lambda}D^{(+)\dagger}_{\bi{q}\lambda}a^{\dagger}_{\bi{q}\lambda}$, $C_{\mathrm{b}}^{\gamma}=\rmi\sum_{\bi{q}\lambda}D^{(-)\dagger}_{\bi{q}\lambda}a^{\dagger}_{\bi{q}\lambda}$, $C_{\mathrm{c}}^{\gamma}=C_{\mathrm{a}}^{\gamma\dagger}$ and $C_{\mathrm{d}}^{\gamma}=C_{\mathrm{b}}^{\gamma\dagger}$ where $D^{(\pm)}_{\bi{q}\lambda}$ are given by (16qsa) and (16qsb). Then, from the expression for the general dissipator (16qsuvz), one finds that the emission (E) and absorption (A) rates from the $\ket{-}$ to and from the ground state are given by

[TABLE]

In the following subsections we aim to calculate $\Xi_{\mathrm{bb}}^{\gamma}(\delta_{-})$ and $\Xi_{\mathrm{dd}}^{\gamma}(-\delta_{-})$. This first requires us to derive expressions for the two-time expectation functions over the environment space operators, which we do in the next subsection. Subsequently we must integrate these expressions over time.

A.1 Photon and phonon expectation values.

We first note that the environment space operators $C_{\alpha}^{\gamma}$ live in both the photon and phonon spaces, and that the expectation value in the ECF definition (16qsuvaa) is over both environments. The Born approximation allows us to separate the expectation value over both environments into a product of averages over single environments. Introducing the concise notations,

[TABLE]

where $i,j$ labels the monomers, we can write

[TABLE]

where in the interaction picture $a_{\bi{q}\lambda}(t)=a_{\bi{q}\lambda}\rme^{-\rmi\nu_{\bi{q}}t}$ and

[TABLE]

As is shown in [21], $\langle B^{+}_{i}(t)B^{-}_{j}(0)\rangle_{\mathrm{E},{\mathrm{pn}}}=\langle B^{-}_{i}(t)B^{+}_{j}(0)\rangle_{\mathrm{E},{\mathrm{pn}}}$, and therefore we can write that

[TABLE]

Aside from the averages over the photon operators, the expressions for the emission (16qsuvacaiataw) and absorption (16qsuvacaiatay) expectation values are largely identical. The thermal averages of the photonic operators are

[TABLE]

We next want to evaluate $\sum_{\bi{q}\lambda}G_{ij}$ for $i,j=1,2$. This is completed by first taking the continuum limit,

[TABLE]

where $\rmd\Omega_{\mathrm{q}}$ is an infinitesimal solid angle. Then, by defining a coordinate system one can show that in this limit

[TABLE]

where $\gamma_{j}(\omega)$ and $\mathcal{F}(x)$ are defined in (16qsuvacad) and (16qsuvacae) respectively. This part of a similar derivation is discussed in detail in [30].

We also need the expectation values of the phonon displacement operators. Using standard results derived in [21] we have that:

[TABLE]

where $\phi_{j}(t)$ is defined by (12).

We now have the ingredients to calculate the ECFs. These require us to evaluate double integrals over photon frequency and time, which contain functions originating from both photon and phonon environments. Using the expressions for the expectation values, we find that these integrals are of the form

[TABLE]

where we have defined $F_{ij}^{+}(\nu_{\mathrm{q}})=F_{ij}(\nu_{\mathrm{q}})N(\nu_{\mathrm{q}})$ and $F_{ij}^{-}(\nu_{\mathrm{q}})=F_{ij}(\nu_{\mathrm{q}})[1+N(\nu_{\mathrm{q}})]$. This enables us to succinctly write the emission (16qsuvacaiata) and absorption (16qsuvacaiatb) rates of the $\ket{-}$ as

[TABLE]

where we have used the fact that $I^{\pm}_{12}(\omega)=I^{\pm}_{21}(\omega)$.

The integrals $I^{\pm}_{12}(\mp\omega)$ lead to the collective effects induced by the monomer coupling on the $\ket{-}$ rates with energy $\omega$, whereas $I^{\pm}_{jj}(\mp\omega)$ lead to the independent monomer rates for a monomer-$j$ with splitting $\omega$.

A.2 Monomer collective effects.

$I^{\pm}_{12}(\mp\omega)$ is calculated using the standard identity

[TABLE]

where $\delta(\epsilon)$ is the Dirac delta function and $\mathcal{P}$ indicates that the principal value of the subsequent integral should be taken. Using this it is straightforward to show that

[TABLE]

The real part of this quantity characterizes the changes to the eigenstate photon absorption and emission rates, and the imaginary part describes the renormalization of the eigenfrequency $\delta_{-}$, which does not diverge and is given explicitly in B.1. Since these both originate from the dipole-dipole coupling between the monomers, they have a significant dependence on the orientations of the dipoles. The entire expression is renormalized by the phonon coupling strength of both monomers to their respective environments through $\kappa_{1}\kappa_{2}$.

A.3 Monomer phonon-renormalized electronic rates.

Finally, we need to calculate the independent rates for each monomer. The relevant integrals are

[TABLE]

Unlike in the collective case, we see from (16qsuvacaiatazbc) that the renormalization induced by either phonon environment has complicated time dependence. A thorough investigation into this integral will be published in an additional paper [31]. We find that

[TABLE]

which is given in the main text (16qsuvacaia). To reach an analytic solution for $\zeta_{j}^{\pm}(\mp\omega)$, we can approximate the continuum of vibrational modes, defined by the spectral density, as a single mode. Specifically, $J(\omega)=\sum_{k}\left|g_{k}\right|^{2}\delta(\omega-\omega_{k})\to\left|\tilde{g}\right|^{2}\delta(\omega-\tilde{\omega})$, where $\tilde{g}$ and $\tilde{\omega}$ are the single mode representation coupling strength and frequency. Estimations for these parameters can be found from the moments of the distribution, $\int_{0}^{\infty}\rmd\omega\left(J(\omega)/\omega^{2}\right)\omega^{n}\equiv m^{(n)}$. One can identify that $m^{(0)}=\phi(0)$ and $m^{(1)}=\lambda$ which are the reorganization energy (10) and phonon propagator (12) respectively. We can then show that $\left|\tilde{g}\right|^{2}=m^{(0)}=\phi(0)$ and $\tilde{\omega}=m^{(1)}/m^{(0)}=\lambda/\phi(0)$. Using this, we find that the phonon renormalization at zero temperature is

[TABLE]

This renormalization can be calculated exactly for certain simple phonon spectral densities. Here, for a single two level system, we test our approximate solution of $\zeta^{\pm}(\omega)$ compared to the exact version for $J(\omega)=A\omega^{3}\mathrm{exp}\left(-\omega/\Omega\right)$. In Figure 6 we plot the approximated and exact expressions $\zeta^{\pm}(\mp\delta)$ against renormalization energy, $\lambda$ for various cut-off frequencies, $\Omega$. The figure shows reasonable agreement between the two solutions, bearing in mind that under the usual FSDA or weak phonon coupling limit $\zeta^{\pm}(\mp\omega)=1$ for all reorganization energies. This justifies our using of the approximated form (16qsuvacaiatazbj) throughout the paper where we use the spectral density given by (16qsuvacaiaj) for which there is no exact expression for the phonon renormalization.

Appendix B Master equations

We write, \dot{\mbox{\boldmath\rho}}_{\mathbf{S}}=\left(\mathcal{L}^{\mathrm{S}}+\mathcal{L}^{\gamma}+\mathcal{L}^{\mathrm{pn}}\right)\mbox{\boldmath\rho}_{\mathbf{S}} where each Liouvillian $\mathcal{L}^{x}$ is a matrix, and the vector \mbox{\boldmath\rho}_{\mathbf{S}}=\left(\rho_{00},\rho_{++},\rho_{--},\rho_{+-},\rho_{-+}\right) where $\rho_{\alpha\beta}\equiv\bra{\alpha}\rho_{S}\ket{\beta}$. Explicitly, $\mathcal{L}^{\mathrm{S}}\rho_{\mathrm{S}}=-\rmi\left[H_{\mathrm{S}},\rho_{\mathrm{S}}\right]$ and $\mathcal{L}^{x}\rho_{\mathrm{S}}=\mathcal{D}^{x}(\rho_{\mathrm{S}})$ for $x=\gamma,\mathrm{pn}$. We ignore the other elements of the density matrix because the dimer level populations decouple from them.

B.1 Photon + system Liouvillian

The sum of the system and photon Liouvillians, $\mathcal{L}^{{\mathrm{S}}+\gamma}=\mathcal{L}^{\mathrm{S}}+\mathcal{L}^{\gamma}$ are found as

[TABLE]

where the absorption rates from the ground state to the $\ket{+}$ or the $\ket{-}$ are denoted $\gamma^{\mathrm{A}}_{\pm}\equiv\Gamma^{\mathrm{A}}_{\pm}N_{\pm}$ and the emission rates from the eigenstates to the ground state are $\gamma^{\mathrm{E}}_{\pm}\equiv\Gamma^{\mathrm{E}}_{\pm}\left(1+N_{\pm}\right)$. The form of $\Gamma^{\mu}_{-}$ is given in the main text for $\mu=\mathrm{A,E}$ (16qsuvacaiam), $N_{\pm}$ is the Bose-Einstein distribution evaluated at frequency $\delta_{\pm}$ respectively, and

[TABLE]

The renormalized excited state detuning of the dimer is $\tilde{\Delta}\equiv\tilde{\delta}_{+}-\tilde{\delta}_{-}$, where

[TABLE]

and the Lamb shift of the excited state of monomer-$j$ is

[TABLE]

where $\mathcal{P}$ denotes the principal value. This integral diverges but it is assumed that the monomer Lamb shifts are negligible. The coherence generating function (CGF) is defined by the integral

[TABLE]

where in the second line we have made the approximation of zero photon temperature. In (16qsuvacaiatazbq) we have defined,

[TABLE]

[TABLE]

and we have made use of the standard functions

[TABLE]

Finally, the rates which measure the degree of communication of the excited states with the coherences between them have the forms

[TABLE]

where we have defined

[TABLE]

and

[TABLE]

and therefore note that $S_{+-}(-\omega)\approx 0$. The Liouvillian shows that due to the imaginary parts of the rates $\Theta_{\pm}$, there are oscillations between the excited state populations and coherences at frequencies $\nu_{\pm}=2\kappa_{1}\kappa_{2}\cos\chi\sqrt{\gamma_{1}(\delta_{\pm})\gamma_{2}(\delta_{\pm})}\mathcal{G}(\delta_{\pm}r_{12})$. This is the origin of the naming for the CGF. The coherences also have additional oscillations occuring at $\tilde{\Delta}$. If we set $\Theta_{+}=\Theta_{-}$ (so that $\nu_{+}=\nu_{-}\equiv\nu_{\pm}$), then it can be shown that the combination of these two effects manifests itself in the excited state populations as oscillations with frequency $\nu_{\mathrm{pop}}=\sqrt{\tilde{\Delta}^{2}+4\nu_{\pm}^{2}}$. There are corrections to this expression when $\nu_{+}\neq\nu_{-}$.

For means of comparison, here we show the Liouvillian for the case where the rotating wave approximation (RWA) is not made in the original photon interaction Hamiltonian. In this case it also becomes necessary to include the double excited state, and only at the very end approximating that it is never populated. This way we can still capture the virtual processes connecting the single excitation eigenstates that pass through the double excited state. We emphasize that we have checked that by making these approximations in the paper we have not lost any information. In this case, the full Liouvillian is

[TABLE]

where without the RWA the CGF is now given by

[TABLE]

where we have again made the approximation of zero photon temperature. Furthermore, the Lamb shifts become

[TABLE]

but are still assumed to be zero. The additional terms in (16qsuvacaiatazcl) that originate from virtual processes passing through the double excited state are of the form

[TABLE]

and the eigenstate splitting gets further renormalized to $\tilde{\Delta}_{f}=\tilde{\delta}_{+,f}-\tilde{\delta}_{-,f}$ where

[TABLE]

B.2 Phonon Liouvillian

We find that

[TABLE]

where the rates are given by

[TABLE]

with $\gamma^{\mathrm{pn}}_{ab}(\omega)\equiv 2\Re[\zeta_{ab}(\omega)]$ and $S^{\mathrm{pn}}_{ab}(\omega)\equiv\Im[\zeta_{ab}(\omega)]$ for $a,b=x,y,z$. The ECFs are

[TABLE]

where $B_{a}(t)$ are given by (16qsuva, 16qsuvb, 16qsuvc) for $a=x,y,z$. Using the properties of displacement operators derived in [21] we find that

[TABLE]

where $\phi(t)$ is defined by (12) and all other combinations equate to zero. The ECF integrals are computed numerically.

The Liouvillian (16qsuvacaiatazdd) describes phonon induced excitation from the $\ket{-}$ to the $\ket{+}$ at rate $\gamma_{+}^{\mathrm{pn}}(-\eta)$ and the reverse process at rate $\gamma_{+}^{\mathrm{pn}}(\eta)$. It also shows that the eigenstate detuning is increased by $\bar{\mu}_{+}^{\mathrm{pn}}$. Furthermore, the population eigenstates either diminish or gain through the coherences due to the $\gamma_{\mathrm{xz}}^{\mathrm{pn}}(0)$ term. The sign of this rate determines which population gains and which diminishes. However, unlike the photon interaction, the phonon interaction induces oscillations between the coherences at frequency $\bar{\mu}_{-}^{\mathrm{pn}}$. Finally, the phonon interaction provides additional decay of the coherences at rate $\bar{\gamma}_{+}^{\mathrm{pn}}+2\gamma_{\mathrm{zz}}^{\mathrm{pn}}(0)$.

Appendix C Approximate analytic dark state condition

Minimizing the rate coefficient (16qsuvacaiam) with respect to $z$ allows us to derive the critical dipole ratio for which the destructive interference maximally suppresses the emission rate of the $\ket{-}$. This has previously been calculated for weak vibrational coupling, for homogeneous [5] and heterogeneous [6] dimers. In both of these cases, however, it was also assumed that the coupling between the monomers $C$ was independent of dipole strength. In our calculations, $C\propto z$ and this means that minimizing (16qsuvacaiam) with respect to $z$ can only be done numerically (as we show in Figure 2 of the main text). To derive simple analytic formulae for the location of the dark state and its efficacy once formed we approximate $C^{\prime}$ to be independent of $z$. Really, though, this is just a necessity to find analytic formulae. Despite this, the formulae we derive are largely obeyed when the correct dependence on $z$ is reinstated. To do so, we write $\gamma_{2}(\omega)=z^{2}\gamma_{1}(\omega)$ in (16qsuvacaiam) and find the minimum of the rate with respect to $z$. We find that this occurs when $z$ is equal to

[TABLE]

In this approximation, when a dimer is formed with dipole magnitude ratio $z_{c}$ the $\ket{-}$ emission rate will be minimal, with the rate coefficient equal to

[TABLE]

Evidently in this simple case, maximal optical decoupling of the minus state occurs for weaker phonon coupling and geometries for which $\mathcal{F}$ is closer to $1$. From (16qsuvacaiatazdm) this also means that $z_{\mathrm{c}}\to 1$, leading to homodimers. If the geometry deviates from $\mathcal{F}=1$, or there is some degree of phonon coupling, it is not possible for the dimer to become fully decoupled by utilizing destructive interference. We note that under the same conditions proposed in [6] (that is for very weak phonon coupling and a geometry giving $\mathcal{F}=1$) we recover the same condition on $z_{\mathrm{c}}$ and the resulting $\Gamma_{\mathrm{c}-}^{\mathrm{E}}$ for perfect dark state formation.

When the dipole-dipole coupling $C$ takes its $z$-dependent form (2), the qualitative dependencies of the critical dipole ratio $z_{\mathrm{c}}$ (16qsuvacaiatazdm) and rate coefficient $\Gamma_{\mathrm{c}-}^{\mathrm{E}}$ (16qsuvacaiatazdn) on system parameters are still followed when the dark state forms by optimizing destructive interference. This can be seen by comparing these expressions to Figure 2. For example, we still find that maximal decoupling occurs when $\mathcal{F}=1$ and for weak phonon couplings. However, there are a couple of important consequences of using the full expression of $C$: (a) a minimum with respect to $z$ does not necessarily exist in the rate for a given set of parameters, forming only for fairly homogeneous dimers; and (b) the Bose-Einstein distributions $N(\delta_{-})$ now have $z$ dependence, however, we have found numerically that the $z_{\mathrm{c}}$ of the eigenstates have negligible dependence on photon temperature.

Appendix D Absorption and emission spectra

D.1 Heterogeneous dimer spectra

In Figure 7 we plot absorption and emission spectra for the heterodimer, specifically for parameters at the HE1 power maximum in Figure 4. So we can compare with the spectra in the main paper for the homodimer (Figure 3) we choose the same geometries in each row, indicated by $\mathcal{F}$, and the same phonon coupling strength and dipole-dipole coupling strength, with $\lambda=100$ meV and $C^{\prime}=\left(100\ \mathrm{meV}\right)z$. Compared to the homodimer, emission in the heterodimer has much smaller dependence on orientation. This is because the heterodimer is not reliant on destructive interference to create a dark state, and therefore is not sensitive to perturbations from ideal geometries. However, the destructive interference is still playing a role in reducing emission and absorption directly into the dark state which can be seen by the changing value of $P_{\mathrm{d}}$ with orientation.

D.2 Calculating spectra

The emission and absorption spectra are calculated as the Fourier transforms of the two-time correlation functions between the positive and negative frequency components of the electric field. It is well documented that, in the Markov approximation, the absorption ($\mu=\mathrm{A}$) and emission ($\mu=\mathrm{E}$) spectrum are given by

[TABLE]

where $\gamma_{ij}$ are the (unnormalized) monomer photon rates and $g^{\mu}_{ij}(t,t+\tau)$ are the two-time expectation functions discussed shortly [30, 32]. In (16qsuvacaiatazdo), the sum runs through $i,j=1,2$ for the monomers in the dimer. Introducing the notation $\gamma_{i}(\delta_{i}^{\prime})\equiv\gamma_{i}$, the photon rates are given by

[TABLE]

which are the unnormalized single monomer photon rate ($i=j$) and collective rates ($i\neq j$) derived in A. The two-time expectation functions for absorption and emission are

[TABLE]

where $\sigma^{\pm}(t)$ are the raising and lowering operators for monomer-$j$ in the Heisenberg picture [30, 32]. Note that to get the overall absorption of some weak probe-field one must subtract the emission spectrum from the absorption [30, 32]. However, here we are only interested in the raw absorption spectrum. To account for the strong phonon interaction we again move to the polaron frame using the transformation described in Section 2.2. It can be shown that after applying the polaron transform, the Heisenberg picture raising and lowering operators become

[TABLE]

where $B_{j}^{\pm}(s)$ are the displacement operators (7) in the Heisenberg picture. We then state the following two-time expectation values of the displacement operators in the Heisenberg picture (for a similar proof see [21]),

[TABLE]

where $\kappa_{j}$ and $\phi_{j}(s)$ are defined in (11) and (12) respectively, and

[TABLE]

is the conjugate phonon propagator. Therefore, in the polaron frame, the two-time expectation functions (16qsuvacaiatazdq) are

[TABLE]

and

[TABLE]

The two-time expectation values of the lowering and raising operators are evaluated with the quantum regression theorem [21, 30, 32, 33]. This maps these expectation values onto elements of the system density matrix. Therefore, these encode the information on the system dynamics (photon/phonon processes). The phonon sidebands in the spectra arise entirely because of the factors of $\rme^{\phi_{j}(s)}$ and $\rme^{\overline{\phi}_{j}(s)}$ in (16qsuvacaiatazdsdu). So, to produce the spectra in Figure 3 and Figure 7 that do not have phonon sidebands (whilst keeping vibrationally induced transfer processes) we have simply set these factors to equal unity.

D.3 Approximating $f^{\mathrm{pn}}$

The area under the total spectrum is

[TABLE]

where $S^{\mu}(\omega)$ is defined in (16qsuvacaiatazdo). We can immediately evaluate the integral over frequency using the identity $\Re\int_{0}^{\infty}\rmd\omega\rme^{\rmi\omega\tau}=\pi\delta(0)$. Then, the integral over $\tau$ in (16qsuvacaiatazdo) results in the area being given by the integrand of $\pi S^{\mu}(\omega)$ evaluated at $\tau=0$. The fractional emission or absorption from the phonon sideband is defined as

[TABLE]

where $A_{\gamma}^{\mu}$ is the integral (16qsuvacaiatazdsdw) with the the factors of $\rme^{\phi_{j}(s)}$ and $\rme^{\overline{\phi}_{j}(s)}$ set to unity in the spectrum (the part of the spectrum not generated through the phonon sideband). Before we can evaluate this we must briefly introduce the Markovian quantum regression theorem (QRT) which is used to evaluate the two time system correlation functions in (16qsuvacaiatazdsdu) and (16qsuvacaiatazdsdv).

The theory behind the Markovian quantum regression theorem (QRT) is well documented [21, 30, 33], so we will just state the relevant results. One finds from using the QRT, that

[TABLE]

where $\chi_{j0}^{i}(\tau)=\bra{j}\chi(\tau)^{i}\ket{0}$ and $\chi^{i}(\tau)$ is an effective density matrix described by

[TABLE]

for the same total Liouvillian $\mathcal{L}^{\mathrm{T}}$ that describes the evolution of the system density matrix, $\rho_{\mathrm{S}}(t)$ i.e. $\dot{\rho}_{\mathrm{S}}(t)=\mathcal{L}^{\mathrm{T}}\rho_{\mathrm{S}}(t)$. The superscript $i$ refers to the initial conditions of the density matrix which, through the QRT, are found as

[TABLE]

Similarly one can show that

[TABLE]

with initial condition

[TABLE]

We can therefore write general expressions for the areas of the spectrum with and without the phonon sidebands, which are

[TABLE]

where

[TABLE]

The functions $Q^{\mu}_{j}$ are the effective density matrices evaluated at zero time. Specifically,

[TABLE]

For algebraic ease, we set $\kappa_{1}=\kappa_{2}\equiv\kappa$ so that $\alpha^{\mu}=\kappa^{2}$ for $\mu=\mathrm{A},\mathrm{E}$. Therefore, we arrive at the general expression

[TABLE]

where $\nu^{\mu}=F^{\mu}/G^{\mu}$ is the ratio of collective and individual monomer electronic rates.

The zero time values of the effective density matrices concerned with absorption, $\Lambda_{j0}^{i}(0)$, do not depend on whether the dimer is homogeneous or heterogeneous and can be readily found as

[TABLE]

where we have introduced the notation $P_{ab}\equiv\bra{a}\rho_{\mathrm{S}}(t=\infty)\ket{b}$. Therefore, without approximation,

[TABLE]

and we arrive at (16qsuvacaiap) in the main text. Calculating the emission fraction is not so simple, because

[TABLE]

where the trigonometric functions describe the rotation between the dipole basis and the eigenbasis and the angle $\chi$ is defined by (16qa) and (16qb). For an exact solution, one can calculate the steady states numerically, however, here we make approximations to get simpler solutions. We will assume that $P_{--}\gg P_{++}$ which is justified if there is a sufficient bias for phonon transfer to move excitations from the $\ket{+}$ to the $\ket{-}$.

We can further simplify by assuming that the dimer is either homogeneous or heterogeneous. If the dimer is homogeneous, then $\Delta\ll C^{\prime}$ and so $\sin\frac{\chi}{2}\approx\cos\frac{\chi}{2}\approx\frac{1}{\sqrt{2}}$. Furthermore, because the eigenstate populations completely decouple from their coherences in the homodimer we can assume that $P_{+-}=P_{-+}=0$. Under these assumptions, $\chi_{j0}^{j}(0)\approx\frac{1}{2}P_{--}$ and $\chi_{10}^{2}(0)\approx\chi_{20}^{1}(0)\approx-\frac{1}{2}P_{--}\mathrm{sign}\left[C\right]$. Consequently,

[TABLE]

where in the second equality we have made the further approximation that $\gamma_{1}\approx\gamma_{2}$ (true for highly homogeneous dimers) to arrive at (16qsuvacaiaq) in the main text.

If the dimer is heterogeneous ($\Delta\gg C^{\prime}$), one cannot derive an expression that does not involve the steady states of the density matrix. Defining $\epsilon=C^{\prime}/\Delta\ll 1$, and noting that $\gamma_{1}\gg\gamma_{2}$, we find that

[TABLE]

where we have used that $\Im\left[P_{+-}\right]\approx 0$. In Figure 8 we plot these derived expressions for $f^{\mathrm{pn}}_{\mu}$ with $\mathcal{F}$ for $\lambda=100$ meV. Additionally, we plot the values calculated by numerically integrating the homodimer and heterodimer spectra over all frequencies with and without the phonon sidebands.

Appendix E Idealized load model

We model the idealized load illustrated in Figure 1(c) as a two level system called a trap. The evolution of the dimer and trap is described by the total density matrix $\rho(t)\otimes\rho_{\mathrm{t}}(t)$ where $\rho(t)$ and $\rho_{\mathrm{t}}(t)$ are the dimer and trap density matrices. As described in the main text we follow standard Born-Markov procedure and trace out the dimer environment degrees of freedom leaving the relevant ones, $\rho_{\mathrm{rel}}(t)=\rho_{\mathrm{S}}(t)\otimes\rho_{\mathrm{t}}(t)$. Adding the trap with a tensor product was first discussed in [12] in place of the original set-up which added the trap energy levels to the system density matrix [5, 13, 14]. This was to reduce the number of free parameters needed to describe the extraction process. The evolution of the system density matrix has been derived in the previous sections and is governed by $\dot{\rho}_{\mathrm{S}}=\mathcal{L}^{\mathrm{T}}\rho_{\mathrm{S}}$ where $\mathcal{L}^{\mathrm{T}}$ is the total Liouvillian. We can write down a phenomenological equation describing the evolution of the trap as \dot{\rho}_{\mathrm{t}}=-i\big{[}H_{\mathrm{t}},\rho_{\mathrm{t}}\big{]}+\mathcal{D}_{\mathrm{t}}(\rho_{\mathrm{t}}), where the trap decay dissipator is given by

[TABLE]

and the trap Hamiltonian is $H_{\mathrm{t}}=\delta_{\alpha}\ket{\alpha}\bra{\alpha}+\delta_{\beta}\ket{\beta}\bra{\beta}$ where $\delta_{\mathrm{t}}=\delta_{\alpha}-\delta_{\beta}$. As described in the main text, there is then a further phenomenological dissipator added which describes non-radiative extraction of excitons from the minus state to the trap. This must act in both spaces and has the form

[TABLE]

Therefore, the evolution of the combined density matrix is described by

[TABLE]

Appendix F Absolute power output

Figure 9 shows the power output of the dimer calculated in Figure 4 but without dividing by the power output of the benchmark. The significant difference between the two power plots is that there is now a new power maximum for $z=1$ and large detunings, $\Delta$. This new maximum does not appear in the benchmarked power output because the benchmark similarly maximizes power output here. The reason that the power maximizes here is simply that reducing the energy of monomer-2 will increase the steady state Bose-Einstein population of this excited state. However, there is a trade-off because this decrease in energy of monomer-2 will also mean that the trap voltage has decreased by the same amount.

Appendix G Important rates in determining power output

In Figure 10 we show the important rates in understanding the positions of the maxima in Figure 4 and Figure 9. Each set of four contour plots in Figure 10 corresponds to the power plot in Figure 4 and Figure 9 which has the same $C^{\prime}/z$ and $\lambda$ values. In each set of four, we show the $\ket{-}$ emission rate (top left); $\ket{+}$ absorption rate (top right); overall phonon transfer rate from the $\ket{+}$ to $\ket{-}$ (bottom left) and the photon non-secular oscillation frequency from the $\ket{+}$ to the excited state coherences (bottom right).

Interestingly, by looking at how the emission rate from the $\ket{-}$ changes as phonon coupling increases for a given dipole-dipole coupling, one can see the dark state in the homodimer, at $(z,\Delta)=(1,0)$, being destroyed. The optical absorption into the $\ket{+}$ is also larger for the HOX than the HEX which prevents bottle-necking and, though less important than a reduced $\ket{-}$ emission, will increase the power output. This occurs because constructive interference of the $\ket{+}$ absorption is equal in magnitude to the destructive interference of the $\ket{-}$ emission. We also show a photon non-secular rate as an example of what these look like.