Design principles for long-range energy transfer at room temperature

TL;DR

This paper reveals how quantum effects in small molecular clusters can be harnessed to enable efficient, long-range energy transfer at room temperature, challenging traditional assumptions about noise suppression.

Contribution

It introduces design principles leveraging quantum coherence and dark states to enhance macroscopic energy transfer in molecular systems, supported by an exactly solvable model and numerical simulations.

Findings

Quantum coherence can extend energy transfer distances.

Dark states can control dissipation pathways.

Design guidelines for artificial light-harvesting systems.

Abstract

Under physiological conditions, ballistic long-range transfer of electronic excitations in molecular aggregates is generally expected to be suppressed by noise and dissipative processes. Hence, quantum phenomena are not considered to be relevant for the design of efficient and controllable energy transfer over significant length and time scales. Contrary to this conventional wisdom, here we show that the robust quantum properties of small configurations of repeating clusters of molecules can be used to tune energy transfer mechanism that take place on much larger scales. With the support of an exactly solvable model, we demonstrate that coherent exciton delocalization and dark states within unit cells can be used to harness dissipative phenomena of varying nature (thermalization, fluorescence, non-radiative decay and weak inter-site correlations) to support classical propagation over macroscopic distances. In particular, we argue that coherent delocalization of electronic excitations over just a few pigments can drastically alter the relevant dissipation pathways which influence the energy transfer mechanism, and thus serve as a molecular control tool for large-scale properties of molecular materials. Building on these principles, we use extensive numerical simulations to demonstrate that they can explain currently not understood measurements of micron-scale exciton diffusion in nano-fabricated arrays of bacterial photosynthetic complexes. Based on these results we provide quantum design guidelines at the molecular scale to optimize both energy transfer speed and range over macroscopic distances in artificial light-harvesting architectures.

Figures

Tables

Table 1. Table 1: Single pigment vs. aggregate fluorescence. Quantum yield of fluorescence ( ϕ italic-ϕ \phi ) and decay lifetime ( τ subscript 𝜏 \tau_{\text{}} ) of isolated Bchl- a molecules and LH2 complexes, together with their radiative ( Γ rad subscript Γ rad \Gamma_{\text{rad}} ) and nonradiative ( Γ non-rad subscript Γ non-rad \Gamma_{\text{non-rad}} ) decay rates. Adapted from Monshouwer  et al. ( 1997 ) .

	$\varphi$ (%)	${\tau }_{}$ (ns)	${\Gamma }_{\text{rad}}$ (ns-1)	${\Gamma }_{\text{non-rad}}$ (ns-1)
Bchl-a	18	3.140	0.057	0.26
LH2	9.86	0.986	0.10	0.91

Table 2. Table 2: Model parameters for a B850 ring.

	Value	Parameter description
$\theta$	${10.06}^{\circ }$	intra-dimer angle Caycedo-Soler et al. (2017)
$\Delta {\theta }_1$, $\Delta {\theta }_2$	${19.9}^{\circ }$, ${17.6}^{\circ }$	dipole tilt away from the tangent Caycedo-Soler et al. (2017)
${\phi }_1$, ${\phi }_2$	${6.24}^{\circ }$, ${4.81}^{\circ }$	dipole tilt away from ring plane Caycedo-Soler et al. (2017)
$r$	3.0 nm	ring radius Walz et al. (1998)
$\Delta z$	0.063 nm	vertical displacement between $\alpha$ and $\beta$ Bchl Caycedo-Soler et al. (2017)
$d$	6.4 D	$Q_y$ dipole moment of Bchl Pullerits et al. (1997)
$\epsilon$	12330 cm-1	average site energy Caycedo-Soler et al. (2018)
$J_1$, $J_2$	320, 255 cm-1	intra- and inter-dimer nearest neighbor coupling Scholes et al. (1999)
$\gamma$	(20 fs)-1	single pigment optical dephasing rate Caycedo-Soler et al. (2018)
$R_{01}$	(200 fs)-1	relaxation rate between the two lowest excitons Cohen Stuart et al. (2011)
${\Gamma }_{\text{th}}$	(1 ns)-1	decay rate from thermalized exciton manifold Monshouwer et al. (1997)
$T$	298 K	temperature of the vibrational bath
$\varphi$	0.1	quantum yield of fluorescence Monshouwer et al. (1997)
${\sigma }_{\text{p}}$, ${\sigma }_{\text{0}}$	265 cm-1, 33 cm-1	std. dev. of uncorrelated/ground state disorder Caycedo-Soler et al. (2018)111These values have been rescaled from Caycedo-Soler et al. (2018) to match the B850 absorption linewidth.
$r_{\text{c}}$	0.5 nm	correlation radius of local pigment environments Renger and Marcus (2002)

Equations

$${\hat{H}}^{(n)}=\sum _i{ϵ}_i|i_n⟩⟨i_n|+\sum _{i\ne j}{J}_{ij}|i_n⟩⟨j_n|,$$

$$V_{ij}<\hslash \gamma <J_{12},$$

$${\partial }_t{p}_{\alpha }^n=$$

$$-\sum _{{\alpha }^{\prime }(\ne \alpha )}{R}_{{\alpha }^{\prime }\alpha }{p}_{\alpha }^n+\sum _{{\alpha }^{\prime }(\ne \alpha )}{R}_{\alpha {\alpha }^{\prime }}{p}_{{\alpha }^{\prime }}^n+I_{\alpha }^n,$$

$$W_{\alpha \beta }=\{\begin{array}{ll}\frac{2|V_{\alpha \beta }{|}^2}{{\hslash }^2}\frac{1}{{\gamma }_{\alpha \beta }}& \text{ for α=β, }\\ \frac{2|V_{\alpha \beta }{|}^2}{{\hslash }^2}\frac{{\gamma }_{\alpha \beta }}{{\gamma }_{\alpha \beta }^2+4J_{12}^2/{\hslash }^2}& \text{ otherwise},\end{array}$$

$${\Gamma }_{\alpha }={\chi }_{\alpha }{\Gamma }_{\text{rad}}+{\kappa }_{\alpha }{\Gamma }_{\text{non-rad}}.$$

$${\Gamma }_{\text{th}}=\sum _{\alpha }{p}_{\alpha }^{\text{th}}{\Gamma }_{\alpha },$$

$$\varphi =\frac{{\sum }_{\alpha }{\chi }_{\alpha }{p}_{\alpha }^{\text{th}}{\Gamma }_{\text{rad}}}{{\Gamma }_{\text{th}}}.$$

$$\overline{W}$$

$$\overline{\Gamma }$$

$$\frac{{\bar{p}}_d}{{\bar{p}}_b}=\frac{R_{db}+2W_{db}}{R_{bd}+2W_{bd}},$$

$${\partial }_t\rho =-\frac{i}{\hslash }[{\hat{H}}^{(1)}+{\hat{H}}^{(2)}+{\hat{V}}^{(1,2)},\rho ]+({\mathcal{D}}^{(1)}+{\mathcal{D}}^{(2)})\rho ,$$

$$\overline{W}$$

$$\overline{\Gamma }$$

$${\partial }_t{p}_b^n$$

$$-2W_{db}{p}_b^n-R_{db}{p}_b^n+R_{bd}{p}_d^n-{\Gamma }_b{p}_b^n+I_b{\delta }_{n,0}$$

$${\partial }_t{p}_d^n$$

$$-2W_{bd}{p}_d^n+R_{db}{p}_b^n-R_{bd}{p}_d^n-{\Gamma }_d{p}_d^n+I_d{\delta }_{n,0},$$

$${\partial }_t{p}_b(x,t)$$

$$-(2w+R_{db}+{\Gamma }_b)p_b(x,t)+(2w+R_{bd})p_d(x,t)$$

$${\partial }_t{p}_d(x,t)$$

$$+(2w+R_{db})p_b(x,t)-(2w+R_{bd}+{\Gamma }_d)p_d(x,t).$$

$${\mathsf{G}}_2$$

$${\mathsf{G}}_0$$

$$\partial_{t}\mathsf{p}(x,t)=\big{(}\mathsf{G}_{2}l^{2}\partial_{x}^{2}-\mathsf{G}_{0}\big{)}\mathsf{p}(x,t)+\mathsf{I}\delta(x).$$

$$\hat{\mathsf{p}}(q,\infty )=\frac{1}{{\mathsf{G}}_0+q^2{l}^2{\mathsf{G}}_2}\mathsf{I}=\sum _{k=0}^{\infty }(-q^2{l}^2{\mathsf{G}}_0^{-1}{\mathsf{G}}_2{)}^k{\mathsf{G}}_0^{-1}\mathsf{I}.$$

$$\mathsf{G}_{0}^{-1}=\frac{1}{\overline{\Gamma}}\left(\begin{array}[]{cc}\bar{p}_{b}&\bar{p}_{b}\\ \bar{p}_{d}&\bar{p}_{d}\end{array}\right)$$

$$\hat{\mathsf{p}}(q,\infty)=\frac{I_{b}+I_{d}}{\overline{\Gamma}}\frac{1}{1+q^{2}l^{2}\overline{W}/\overline{\Gamma}}\left(\begin{array}[]{c}\bar{p}_{b}\\ \bar{p}_{d}\end{array}\right),$$

$$\mathsf{p}(x,\infty)=\frac{I_{b}+I_{d}}{\overline{\Gamma}}\frac{1}{2l_{\text{diff}}}e^{-|x|/l_{\text{diff}}}\left(\begin{array}[]{c}\bar{p}_{b}\\ \bar{p}_{d}\end{array}\right),$$

$${\hat{\mathsf{p}}}^{\prime }(q,\infty )=\frac{\hat{g}(q)}{{\mathsf{G}}_0+q^2{l}^2{\mathsf{G}}_2}\mathsf{I}=\hat{\mathsf{p}}(q,\infty )\hat{g}(q).$$

$${\mathsf{p}}^{\prime }(x,\infty )=\int \text{d}{x}^{\prime }\mathsf{p}(x^{\prime },\infty )g(x-x^{\prime }).$$

$$H(\mathbf{r},\mathbf{p};\mathbf{R},\mathbf{P})$$

$$+V_{\text{nuc-nuc}}(\mathbf{R})+T_{\text{nuc}}(\mathbf{P}),$$

$$\begin{gathered} [T_{\text{el}}(\hat{\mathbf{p}})+V_{\text{el-el}}(\hat{\mathbf{r}})+V_{\text{el-nuc}}(\hat{\mathbf{r}};\mathbf{R}) \\ +V_{\text{nuc-nuc}}(\mathbf{R})]|{\varphi }_{\mathbf{R}}^i⟩={\epsilon }_{\mathbf{R}}^i|{\varphi }_{\mathbf{R}}^i⟩. \end{gathered}$$

$$⟨{\varphi }_{\mathbf{R}}^i|H(\hat{\mathbf{r}},\hat{\mathbf{p}};\mathbf{R},\mathbf{P})|{\varphi }_{\mathbf{R}}^j⟩={\delta }^{ij}{\epsilon }_{\mathbf{R}}^i+⟨{\varphi }_{\mathbf{R}}^i|T_{\text{nuc}}(\mathbf{P})|{\varphi }_{\mathbf{R}}^j⟩.$$

$${\hat{H}}_{\text{BO}}(\mathbf{R},\mathbf{P})=\sum _i|{\varphi }_{\mathbf{R}}^i⟩⟨{\varphi }_{\mathbf{R}}^i|\left(T_{\text{nuc}}(\mathbf{P})+{\epsilon }_{\mathbf{R}}^i\right).$$

Full text

I Introduction

Over the last decades research into the role of coherent excitonic delocalization in the dynamics in photosynthetic membranes has shown that strong coherent coupling in subunits of tightly coupled pigments can result in short-ranged excitonic delocalization in the steady state Adolphs and Renger (2006); Müh et al. (2007); Novoderezhkin et al. (2006); Strümpfer and Schulten (2009, 2012); Pajusalu et al. (2015). Delocalization within these domains, typically restricted to individual proteins termed antenna complexes, is essential for modeling transient and steady state optical spectra of the full light-harvesting ensemble. Additionally, excitonic delocalization within antenna complexes is a crucial ingredient for a modular description of dynamics over longer distances and timescales Hu et al. (2002); Jang et al. (2004); Caycedo-Soler et al. (2010), which, as observed experimentally Salverda et al. (2003); Herek et al. (2000); Ma et al. (2016); Agarwal et al. (2002), can be interpreted as a series of incoherent transfer steps, described by simple rate processes. These energy-transfer rates depend on the properties of the states involved, and thereby rely heavily on the steady state excitonic delocalization within antenna complexes, as has been shown by numerically exact calculations Strümpfer and Schulten (2009, 2012); Bennett et al. (2018, 2018).

The modular architecture found in biological light-harvesing membranes, whereby antenna complexes containing a few pigments self-aggregate into larger structures, offers the potential for artificial solar energy conversion and molecular electronics based on such a modular design Scholes et al. (2011); Brédas et al. (2004); Baldo and Stojanovic (2009); Aida et al. (2012). This is made possible by the high degree of experimental control that is today available for the integration of synthetic and biological structures and for the directed assembly of photosynthetic antenna complexes isolated from living organisms Escalante et al. (2008, 2010); Eisenberg et al. (2014); Vasilev et al. (2014); Harris et al. (2013); Boulais et al. (2017) or of supramolecular dye arrays Eisele et al. (2012, 2014); Akselrod et al. (2014); Clark et al. (2014); Haedler et al. (2015); Caram et al. (2016); Wan et al. (2017). To facilitate the realization of these devices, however, a theoretical understanding of the mechanisms involved in to-date unexplained observations of the large diffusion lengths in several of these light-harvesting architectures is needed. For instance, the observed micron-scale diffusion of excitations in nano-fabricated arrays of purple bacteria antenna complexes and phycobilisome proteins Escalante et al. (2008, 2010); Eisenberg et al. (2014) exceeds by more than one order of magnitude the theoretical expectation for diffusion based on experimentally determined parameters under physiological conditions. On the one hand, this shows that an important enhancement of the diffusion can be achieved with hybrid technologies. On the other hand, it stresses the necessity for theory to establish and verify physical design principles by which delocalization of electronic excitations over a few pigments can enable the observed long-range energy transfer.

One viable strategy for improving the diffusion lengths of excitons is to extend their lifetime, thus allowing them to propagate across longer distances. Recent studies on the role of dark states in solar energy conversion have provided valuable insight into the potential advantages of protecting excitations from losses due to fluorescence in order to promote charge separation Creatore et al. (2013); Fruchtman et al. (2016); Wertnik et al. (2018); Hu et al. (2018); Rouse et al. (2019). These models, though, do not consider the microscopic origin of the interactions, thereby omitting the conditions that enable or inhibit the active participation of dark states in the dynamics. Because dark states cannot be excited directly by light and typically do not couple efficiently to the propagating bright states, their participation in the exciton propagation is often overlooked and requires a careful reconsideration of the typical models used to describe resonant energy transfer Sissa et al. (2011); Liu et al. (2015). A further challenge is presented by the need of accounting for non-radiative decay due to the interaction between electronic excitations and vibrational motion of pigment molecules, which is responsible for the usually low fluorescence yield of light-harvesting complexes Monshouwer et al. (1997); Helenius et al. (1997). In fact, the competition between radiative and non-radiative decay channels can drastically modify the dissipation landscape, leading to scenarios in which optically dark states have a much shorter lifetime than bright states. Moreover, the dynamics of excitons can be influenced by the presence of correlations between the local vibrational environments of each pigment. Most research on the topic has been focused on clarifying the extent to which these correlations influence the ultrafast spectroscopy of excitons, often leading to contrasting predictions Cao et al. (2020). However, more relevant for our case, weak inter-site correlations can be expected to play a role on much longer timescales too and thus influence non-radiative decay.

In the present work, we provide a theoretical model that identifies the desirable features of spectral structures and exciton delocalization within unit cells (subunits of tightly coupled pigments) that support efficient energy transfer. This model builds on the formation and participation of dark states in the dynamics in order to achieve long-range diffusion across arrays of these subunits. We show how this is possible by a combination of excitonic delocalization within unit cells and close proximity between unit cells. Identifying each LH2 complex as a unit cell of the transfer chain, we show that a model that can explain the long-range energy propagation reported in LH2 arrays does not need to resort to previously hypothesized Lloyd and Mohseni (2010); Abasto et al. (2012); Yang et al. (2012) long-range quantum coherence involving several LH2 complexes Pullerits et al. (1996); Chachisvilis et al. (1997); Monshouwer et al. (1997); Trinkunas et al. (2001); Cheng and Silbey (2006). Indeed, we demonstrate that a theoretical description that is consistent with available experimental data concerning structure and optical response of the LH2 antenna Pullerits et al. (1996); Chachisvilis et al. (1997); Monshouwer et al. (1997); Trinkunas et al. (2001); Cherezov et al. (2006) can be developed, which reproduces the experimentally observed exciton diffusion length Escalante et al. (2010).

This theoretical description therefore presents desirable features for transport across photosynthetic membranes. The low-energy part of the excitonic spectrum of a unit cell comprises states that are protected against dissipation, while high-energy excitons offer fast pathways for energy transfer to neighboring units by virtue of their delocalization. A combination of these two features allows for robust long-range energy transfer.

The remainder of this work is organized as follows. In Section II.A we present the exactly solvable model which contains the necessary features to discuss long-range exciton propagation across a modular array of excitonic unit cells. Here we discuss several energy-transfer regimes and their relation to the underlying excitonic properties. In Section II.B we introduce the bacterial antenna complex LH2 and provide a thorough characterization of their excitons relevant to energy transfer. In Section II.C we present the results of our simulations of a linear array of LH2 complexes and we discuss how the mechanisms presented in Section II.A apply to a real-world application, backed by existing experimental findings. Section III summarizes the main results of this work and sets a broader context in which these can provide useful guidelins for energy transfer design in molecular materials.

II Model, Results and Major Insights

II.1 Minimal model of a light-harvesting array

The beneficial role provided by exciton delocalization and the formation of states protected against dissipation in excitonic energy transfer can be readily understood in a model system made up of dimeric unit cells (Fig. 1a). Let us consider a homogeneous linear array of unit cells, each of them described by a Hamiltonian

[TABLE]

where the state $|i_{n}\rangle$ describes an electronic excitation of the $i$-th pigment belonging to the $n$-th unit cell. Energy transfer across the full array is made possible by the interaction between pigments of different unit cells, described by $\hat{V}^{(n,n+1)}=\sum_{i,j}V_{ij}|i_{n}\rangle\langle j_{n+1}|+\mathrm{H.c.}$, where $V_{ij}$ are inter-dimer dipole-dipole couplings.

If the single-pigment dephasing rate $\gamma$, the intra- and inter-dimer couplings $J_{12}$ and the inter-dimer dipole-dipole couplings$V_{ij}$ follow

[TABLE]

then delocalized intra-dimer excitons defined by $\hat{H}^{(n)}|\alpha_{n}\rangle=E_{\alpha}|\alpha_{n}\rangle$ form. These may be used to describe the (incoherent) energy transfer between adjacent dimers Jang et al. (2004). We refer to “incoherent” or “classical” transfer interchangeably, meaning that inter-unit cell coherences of the density operator of the full chain are negligible ($\langle i_{n}|\rho|j_{m}\rangle\approx 0$ for $n\neq m$), which then results in classical diffusion across subunits exhibiting local coherent dynamics. The hierarchy of interactions in Eq. (2) is commonly fulfilled in photosynthetic membranes and nano-engineered arrays, where pigments aggregate in antenna complexes within sub-nanometer distances, while the inter-complex distances can span several nanometers, as we will discuss later in detail. For the configuration of anti-parallel transition dipoles (i.e. $\mathbf{d}_{1}=-\mathbf{d}_{2}=d\mathbf{e}_{1}$) (indicated by red arrows in Fig. 1a), the $n$-th dimeric unit supports a dark (bright) exciton $|d_{n}\rangle$ ($|b_{n}\rangle$) given by the symmetric (antisymmetric) coherent superposition $|\alpha_{n}\rangle$ of single pigment states, namely $|\alpha_{n}\rangle=|\smash{{}^{d_{n}}_{b_{n}}}\rangle=(|1_{n}\rangle\pm|2_{n}\rangle)/\sqrt{2}$ when $\varepsilon_{1}=\varepsilon_{2}$.

Within this framework, a quantum master equation description of the full chain Hamiltonian in the presence of dephasing and relaxation mechanisms May and Kühn (2004) can be replaced by the classical rate equations

[TABLE]

for the population $p^{n}_{\alpha}$ of the $\alpha$-th exciton on the $n$-th dimeric unit cell. Here pairs of subindices $\alpha\beta$ ($\alpha\alpha^{\prime}$) label excitons on different (the same) unit cells. The rates that describe the transfer of excitations between unit cells $W_{\alpha\beta}$, their overall decay rate $\Gamma_{\alpha}$, and their thermalization rate $R_{\alpha\alpha^{\prime}}$, depend on the characteristics of the quantum states and their environments within these cells. Injection of excitations into the array can be spatially dependent and is taken to occur with rates $I^{n}_{\alpha}$ on site $n$. The rates $W_{\alpha\beta}$ can be obtained from the overlap between homogeneous lineshape functions and depend on the coupling matrix elements $V_{\alpha\beta}=\langle\alpha_{n}|\hat{V}^{(n,n+1)}|\beta_{n+1}\rangle$ between unit cell eigenvectors and the relative dephasing rate $\gamma_{\alpha\beta}$ between these states via

[TABLE]

as explained in more detail in Appendix D. The dephasing rate between non-overlapping excitons $\gamma_{\alpha\beta}=\gamma_{\alpha}+\gamma_{\beta}$ is the sum of the linewidths $\gamma_{\alpha}=\sum_{\alpha^{\prime}}R_{\alpha^{\prime}\alpha}/2$, which are typically dominated by pure dephasing $R_{\alpha\alpha}$ over the intra-unit cell thermalization rates $R_{\alpha^{\prime}\neq\alpha}$. The thermalization rate $R_{\alpha\neq\alpha^{\prime}}\propto\mathcal{J}(|E_{\alpha}-E_{\alpha^{\prime}}|)\ |n(E_{\alpha}-E_{\alpha^{\prime}})|$ is proportional to the phonon spectral density $\mathcal{J}(\hbar\omega)$ and to the thermal boson occupation number $n(\hbar\omega)$ across the excitonic manifold. This leads to a Boltzmann distribution of the excitons when injection, loss and inter unit-cell transfer are much slower than thermalization. Excitonic pure dephasing $R_{\alpha\alpha}/2=\sum_{i}|\langle i_{n}|\alpha_{n}\rangle|^{4}\gamma\equiv\mathcal{P}^{-1}_{\alpha}\gamma$ is slower than the pigment’s pure dephasing $\gamma$ by a factor given by the inverse participation ratio $\mathcal{P}^{-1}_{\alpha}$ Caycedo-Soler et al. (2018). Consequently, excitonic delocalization ($\mathcal{P}_{\alpha}>1$) results in slower dephasing rates for unit cell excitons in comparison to individual pigments, which, based on Eq. (4), implies an enhancement of the transfer rates between excitons in neighboring units with increasing unit cell delocalization. For this dimeric system where $\mathcal{P}_{\alpha}=2$, delocalization over two pigments of states $|b_{n}\rangle$ and $|d_{n}\rangle$ results in a two-fold speed up on the inter-dimer transfer $W_{dd}$ and $W_{bb}$, as compared to unit cells where strong dephasing $\gamma$ or mild coupling $J_{12}$ prevents the formation of delocalized excitons $|\alpha_{n}\rangle$. Excitonic delocalization within unit cells also redistributes the optical transition dipole strength of individual pigments, resulting in a fluorescence rate of excitons $\chi_{\alpha}\Gamma_{\text{rad}}$, where $\Gamma_{\text{rad}}$ is the single pigment fluorescence rate and $\chi_{\alpha}=\sum_{ij}\langle\alpha_{n}|i_{n}\rangle\chi_{ij}\langle j_{n}|\alpha_{n}\rangle$ characterizes the optical brightness of an exciton. The brightness quantifies the number of sites participating in the fluorescence from exciton $\alpha$ Meier et al. (1997), and is determined by \chi_{ij}=(\mathbf{e}_{i}\cdot\mathbf{e}_{j})\ j_{0}(2\pi r_{ij}/\lambda)-\frac{1}{2}\big{(}\mathbf{e}_{i}\cdot\mathbf{e}_{j}-3(\mathbf{e}_{i}\cdot\mathbf{n}_{ij})(\mathbf{e}_{j}\cdot\mathbf{n}_{ij})\big{)}\ j_{2}(2\pi r_{ij}/\lambda), where $j_{\nu}$ are spherical Bessel functions of the first kind, $\lambda$ is the wavelength associated to the pigment’s $Q_{y}$ transition and $\mathbf{r}_{ij}={r}_{ij}\mathbf{n}_{ij}$ is the relative position of pigments $i$ and $j$. In the limit $r_{ij}\ll\lambda$, the brightness reduces to the usual measure of superradiance $|\mathbf{D}_{\alpha}|^{2}/d^{2}$, i.e. the relative dipole strength, where $\mathbf{D}_{\alpha}=\sum_{i}\langle i|\alpha\rangle\mathbf{d}_{i}$, leading to a superradiant $|b\rangle$ ($\chi_{b}\approx 2$) and a subradiant $|d\rangle$ ($\chi_{d}\approx 0$) exciton for the dimeric unit cells of Fig. 1a.

Two ideas that will play a major role later are worth stressing at this point. First, we note that, while dark states do not couple significantly to electromagnetic fields, they still play a central role in the energy transfer dynamics across different unit cells. Intuitively, one expects that, if the distance between unit cells is much larger than their internal size ($l\gg l_{0}$), two neighboring dimers couple via their global dipoles $\mathbf{D}_{\alpha}$. However, if the dimers are placed sufficiently close such that their distance becomes comparable to their spatial extent ($l\gtrsim l_{0}$) the dipole approximation for their mutual interaction breaks down and states with vanishing dipole strength can start to interact via their higher moments, thus gaining some coupling strength $V_{dd}$. This interaction can even exceed the one of bright states $V_{bb}$ in certain configurations as shown in Fig. 1b, where the coupling strengths $V_{\alpha\beta}$ are plotted against the inter-dimer distance, for the arrangement shown in Fig. 1a. For distances $l\lesssim 3.2l_{0}$, dark states couple more strongly than bright states. As a result, energy transfer through the dark manifold can be achieved much faster than through its bright counterpart. Secondly, we note that typical light-harvesting complexes show a rather low quantum yield of fluorescence $\phi$, i.e. most of the optically generated excitations are lost through non-radiative decay channels. This means that dark states, although protected against fluorescence, could easily be more dissipative than bright states, as we show below. As discussed in Appendix B and C, non-radiative decay rates are influenced by the presence of correlations in the vibrational environments of single pigments. When these correlations are taken into account, the non-radiative decay rate is distributed to different excitons analogously to what happens to the radiative rates (i.e. superradiance). In fact, static correlations between local vibrational environments can be interpreted as arising from the presence of delocalized vibrations coupling to different sites (see Appendix C for a detailed explanation). Under these circumstances, the excitonic decay rates gain a prefactor $\kappa_{\alpha}$ which depends on the phase with which each site contributes to the excitonic wavefunction. Depending on the specific exciton, localized excitations can interfere constructively (destructively) to yield non-radiative rates which can be larger (smaller) than their single-pigment counterpart. In our case, the non-radiative decay rate of the symmetric (antisymmetric) state $|d_{n}\rangle$ ($|b_{n}\rangle$) is enhanced (reduced) by a factor $\kappa_{d(b)}=1\pm e^{-l_{0}/r_{\text{c}}}$ with respect to its single-pigment value $\Gamma_{\text{non-rad}}$, where $r_{\text{c}}$ is the correlation length of the pigment’s vibrational environment. Thus, the decay rate from exciton $|\alpha_{n}\rangle$ is given by

[TABLE]

Note that, for large correlation lengths, the non-radiative decay from the bright state can be significantly reduced compared to the one from the dark state, i.e. $\kappa_{b}\ll\kappa_{d}$. Since both radiative and non-radiative decay in light-harvesting complexes are typically much slower than exciton thermalization, the decay of photoexcitations in a single light-harvesting unit occurs from a thermalized exciton distribution with the average rate

[TABLE]

where $p_{\alpha}^{\text{th}}\propto e^{-E_{\alpha}/k_{\text{B}}T}$. This average decay rate is an experimentally accessible quantity, as is the quantum yield of fluorescence

[TABLE]

The latter only determines the overall relative contribution of radiative and non-radiative decay, whereas the factors $\chi_{\alpha}$ and $\kappa_{\alpha}$ determine how these rates are distributed across the excitonic manifold.

So far, we have discussed how delocalization, transition dipole geometry and environmental correlations determine the dissipation properties of the excitonic states within unit cells, radiative and non-radiative alike, and how the finite size of the unit cell in densely packed arrays opens up the possibility to engage dark states into the propagation dynamics. We now have all the necessary ingredients to determine how all these properties influence the diffusion length of excitons across the light-harvesting array. As we are interested in diffusion over macroscopic distances, it is useful to describe the position of a dimeric unit cell in terms of a continuous variable $x=nl$ as $l\to 0$. Thus, bright and dark state populations at discrete sites $p^{n}_{\alpha}(t)$ are replaced by the respective densities $p_{\alpha}(x,t)$ and Eq. (II.1) takes the form of two coupled continuous diffusion equations, which allow for an analytical solution of the steady-state density $p_{\alpha}(x,\infty)$, as detailed in Appendix A. If we assume that the injection and decay happen on a much slower timescale than transfer and relaxation, the solutions for local driving at $x=0$ take the simple form $p_{\alpha}(x,\infty)=(\bar{p}_{\alpha}/(2l_{\text{diff}}))e^{-|x|/l_{\text{diff}}}$, with $\bar{p}_{b}+\bar{p}_{d}=1$. The diffusion length is given by $l_{\text{diff}}=l\sqrt{\overline{W}/\overline{\Gamma}}$, introducing the effective transfer and decay rates

[TABLE]

as weighted averages on the populations $\bar{p}_{\alpha}$. These are determined by the ratio

[TABLE]

which only depends on the rates that provide mixing between bright and dark manifolds. In particular, we observe that the presence of symmetric inter-dimer bright-to-dark transfer $W_{db}=W_{bd}$ causes a deviation of the steady state populations from thermal equilibrium $R_{db}/R_{bd}=e^{-(E_{d}-E_{b})/k_{\text{B}}T}$, resulting in an increase of the stationary dark state population. Thus, the onset of this non-equilibrium state, which can be observed at short inter-dimer separations $l$ when the rates $W_{bd}=W_{db}$ become comparable to $R_{bd}$ and $R_{db}$ (Fig. 1b, gray dotted line), leads to a modification of effective transfer and decay rates with respect to their equilibrium values.

In order to study the effect of dark states on energy transfer, we focus on the dependence of the diffusion length $l_{\text{diff}}$ on the average intra-dimer relaxation rate $R=\sqrt{R_{bd}R_{db}}$, quantifiyng the strength of local electron-phonon coupling, at a fixed lattice spacing $l$ such that bright-to-dark transfer and relaxation take place on a similar timescale and exciton transfer through dark states is faster than through bright states (Fig. 1b, yellow vertical line).

First, we consider the ideal situation in which the main decay channel is radiative ($\phi=0.9$): Most of the dissipation takes place at the low-energy bright state, which propagates more slowly than the higher-energy dark exciton. In this scenario, long range energy transfer clearly benefits from the establishment of a non-equilibrium state with more population in the dark state. In fact, when reducing $R$, e.g. by decreasing the local exciton-phonon coupling, we are moving population from the slow, dissipative bright state to the fast, less dissipative dark state. Furthermore, we are at the same time reducing the effective decay rate $\overline{\Gamma}$ and increasing the effective transfer rate $\overline{W}$, resulting in larger diffusion length $l_{\text{diff}}$ (Fig. 1c, red line).

In a more realistic scenario, we expect non-radiative decay to play a much bigger role. Let us then set $\phi=0.1$, closer to values observed for biological light-harvesting complexes. In doing so, we also adjust the single pigment decay rates $\Gamma_{\text{rad}}$ and $\Gamma_{\text{non-rad}}$ to ensure that the effectve decay rate at thermal equilibrium $\Gamma_{\text{th}}$ from Eq. (6) remains the same as in the case just considered. Since now the main decay pathway is the non-radiative decay from the fast high-energy exciton, an increase in the relaxation rate $R$ would lead to a longer lived excitation, but also to a slower propagation. However, the former prevails: the reduction of effective decay $\overline{\Gamma}$ is large enough to ensure a longer-ranged propagation (Fig. 1c, blue line). In other words, by increasing the coupling to local phonons, we counter-intuitively increase the spatial extent of energy transfer, at the price of making it slower. This energy transfer regime can be seen as a natural extension to systems with realistic non-radiative decay channels of the dark state protection scheme proposed in context of quantum photocells Creatore et al. (2013); Fruchtman et al. (2016); Wertnik et al. (2018); Hu et al. (2018); Rouse et al. (2019) and recently applied to energy transfer Davidson et al. (2020). Our generalized “dark” state protection mechanism makes use of states that are protected against non-radiative decay (and are therefore “dark”), to extend the lifetime of the excitons, which can diffuse across longer distances.

Lastly, we observe that, while non-radiative decay pathways seem to fundamentally limit the transport efficiency of light-harvesting architectures, one can work around this constraint by expoiting the symmetry of the excitonic wavefunction. Let us consider an arrangement in which now the dipoles within each dimeric unit cell are parallel. In this situation, the intra-dimer coupling $J_{12}$ becomes negative and the low energy exciton is now the symmetric state $(|1_{n}\rangle+|2_{n}\rangle)\sqrt{2}$. While most of the optical dipole still resides in the low energy level, now this state also becomes the most sensitive to non-radiative decay. Thus, we are left with a high-energy exciton which shows little dissipation and fast transfer. As shown in Fig. 1c (green line), this situation is similar to the one that we considered with $\phi=0.9$. The important difference, however, is that in this case non-radiative decay is fully taken into account, but its effect is mitigated by exploiting excitonic delocalization within a unit cell and the ensuing redistribution of decay pathways.

II.2 Bacterial light-harvesting units

The simple dynamical model considered so far allows for analytical expressions that facilitate the identification of the different mechanisms on which we base our explanation of the to-date unexplained experimental observations of long-range energy transfer in nano-engineered arrays of LH2 photosynthetic complexes of purple bacteria Rb. sphaeroides Escalante et al. (2010). These complexes consist of a protein holding two concentric bacteriochlorophyll (Bchl) rings (Fig. 2a) with the inner B850 ring consisting of 18 strongly interacting Bchl pigments, which at room temperature exhibit an excitonic delocalization across about 3–6 pigments, determined by superradiance measurements Pullerits et al. (1996); Chachisvilis et al. (1997); Monshouwer et al. (1997); Trinkunas et al. (2001); Cheng and Silbey (2006). This subunit mediates the transfer between LH2 complexes under physiological conditions, whereas the B800 pigments supports localized excitations which extend the absorption range of the LH2 complex and regulate oxidation Polívka et al. (2007).

We consider the full B850 Hamiltonian $\hat{H}^{(n)}$ describing the interactions among $Q_{y}$ transitions of its $18$ BChls, sketched with red arrows in Fig. 2a, and study the delocalization properties of the single-ring excitons $|\alpha\rangle$, with $\alpha=1,\dots,18$. For realistic LH2 complexes, we need to consider different excitonic energies $E_{\alpha}$ from realizations of pigment energies $\epsilon_{i}$ in order to describe the inhomogeneities (static disorder) arising from local protein configurations. Our choice of spectral density, nearest neighbor couplings, static disorder, geometry and magnitude of the transition dipoles are justified by previous independent analysis of experimental observations Monshouwer et al. (1997); Cohen Stuart et al. (2011); Koolhaas et al. (1998); Hofmann et al. (2004); Jang et al. (2011); Pajusalu et al. (2015); Caycedo-Soler et al. (2018) and, when incorporated into our model, they reproduce observed absorption spectra and superradiance enhancement as shown in Fig. 2b–c (for details of calculation and parameters, see Appendix D and E).

This parametrization of the B850 ring results in cooperative fluorescence from $\langle\sum_{\alpha}\chi_{\alpha}p_{\alpha}^{\text{th}}\rangle\approx 3$ pigments on average (where $\langle\cdot\rangle$ represents the average over static disorder), which is consistent with the experimental observations of superradiance in LH2 Monshouwer et al. (1997). Moreover, fixing the decay rate for an isolated LH2 complex and the quantum yield of fluorescence to the experimentally observed values of $\Gamma_{\text{th}}=(1\text{ ns})^{-1}$ and $\phi=0.1$, allows us to obtain the distribution of decay rates $\Gamma_{\alpha}$ across the exciton manifold. A modest value of $r_{\text{c}}=5\text{ \AA}$ for the correlation radius of the pigment’s vibrational environments Renger and Marcus (2002); Adolphs and Renger (2006) is sufficient to generate a significant redistribution of the dissipation rates to the higher-energy part of the excitonic manifold (Fig. 2c), as we would expect for an anti-parallel arrangement of neighboring dipoles. The participation ratio shown in Fig. 2c exhibits a maximum $\langle\mathcal{P}_{\alpha}\rangle\approx 8$, also consistent with previous estimates Monshouwer et al. (1997), which underlie excitonic delocalization constrained to small portions of the B850 ring (Fig. 2a).

We should note that for realistic LH2 center-to-center distances $l$, the aggregation into arrays does not disrupt the excitonic manifold of single rings, as can be expected from the similarity of optical spectra of diluted and densely packed arrays Escalante et al. (2010). Typical physiological conditions and lipid-reconstituted membranes exhibit a center-to-center distance $l\approx 8$ nm, with increasing inter-complex distances for larger lipid concentration Pflock et al. (2008, 2011a, 2011b). On the other hand, the process of nano-fabrication of LH2 arrays exploits host-guest interactions on a nano-imprinted substrate and does not involve lipids Escalante et al. (2010), which allows us to assume the $6.2$ nm diameter of LH2 $\beta$-helices of Rb. sphaeroides Walz et al. (1998) as the absolute minimum for $l$. Hence, it is reasonable to consider center-to-center separations of $l\gtrsim 6.5$ nm in the nano-engineered arrays.

To assess the robustness of the single ring excitonic manifold against the coherent interaction between neighboring LH2s, we proceed to diagonalize the full two-ring Hamiltonian $(\hat{H}^{(1)}+\hat{H}^{(2)}+\hat{V}^{(1,2)})|\tilde{\mu}\rangle=E_{\tilde{\mu}}|\tilde{\mu}\rangle$ and present in Fig. 2c the average over static disorder of their participation ratios, relative dipole strengths and dissipation rates. Notice that even though the participation ratio in the LH2 pair $\mathcal{P}_{\tilde{\mu}}$ is slightly larger, very minor changes occur in the distributions of optical brightness $\chi_{\tilde{\mu}}$ and dissipation rates $\Gamma_{\tilde{\mu}}$ with respect to the single ring eigenstates. As a consequence, optical absorption spectra are only slightly affected by the coherent interaction between LH2 rings (Fig. 2b). This result confirms that the coherent electronic interaction for realistic values of $l$ does not perturb significantly the excitonic structure of isolated rings. The robustness of the single ring excitonic manifold can be understood by noticing that the maximum coupling between any two LH2 excitons residing on different rings (averaged over static noise and relative rotations on coplanar rings) is below 40 cm*-1* even for $l=6.5\text{ nm}$, which is much smaller than the nearest neighbor interactions within each ring $\approx 250\textrm{--}350$ cm*-1* Koolhaas et al. (1998); Caycedo-Soler et al. (2018). One last argument in favor of the robustness of the single ring excitonic manifold to the coherent coupling between rings is provided by looking at the residual inter-ring coherence after exciton thermalization. To do so, we set up a Lindblad equation describing the coherent interaction between two rings and local thermalization and dephasing,

[TABLE]

and solve for the steady state $\rho^{\text{ss}}$. Since both radiative and non-radiative decay take place on a much slower timescale, we neglect them here. The constant of proportionality between the thermalization rates $R_{\alpha\alpha^{\prime}}$ and the spectral density $\mathcal{J}(|E_{\alpha}-E_{\alpha^{\prime}}|)$ estimated via fluorescence line narrowing experiments Renger and Marcus (2002) is such that the $\approx$ 200 fs timescale of equilibration in LH2 Cohen Stuart et al. (2011) is reproduced. As we are interested in demonstrating that the stationary inter-ring coherence is typically negligible, we average the absolute value of the matrix elements $|\rho^{\text{ss}}_{\alpha\beta}|$ over static disorder to avoid ensemble dephasing. The results are shown in Fig. 2d for two inter-ring distances. Even for the shortest distance ($l=6.5\text{ nm}$), the average steady state coherence $\langle|\rho^{\text{ss}}_{\alpha\beta}|\rangle$ between any two single-ring excitons $|\alpha_{1}\rangle$ and $|\beta_{2}\rangle$ is much smaller than $1/2$, which is the value it would take for a maximally coherent superposition of two states $(|\alpha_{1}\rangle+|\beta_{2}\rangle)/\sqrt{2}$.

At this point, we have established that for center-to-center distances $l\geq 6.5$ nm, the incoherent energy transfer between neighboring B850 rings can be treated based upon the single ring eigenstates. Now we proceed to analyze the mechanisms that underlie this incoherent energy transfer between neighboring LH2s. In Fig. 2e we show how different single-ring excitons residing on two neighboring LH2 rings are connected via the average transfer rates $\langle W_{\alpha\beta}\rangle$. As we saw in the previous section, a shorter separation between unit cells leads to a substantial participation of optically dark states in the energy process, which can surpass the bright states in terms of transfer speed. In fact, also in this case, when reducing the distance $l$ from 8.5 nm to 6.5 nm, the fastest transfer pathways (red dashed lines) shift from the low energy part of the spectrum, where most of the dipole strength resides, to higher-energy excitons, which show larger delocalization.

This finding underlines that the interaction in densely packed arrays does not depend on the exciton transition dipoles but rather on their delocalization, as quantified by the participation ratio.

II.3 Nano-engineered LH2 arrays

Armed with these facts, we now proceed to discuss the origin of the long-range diffusion observed in Escalante et al. (2010). In this experiment, simultaneous excitation with a continuous-wave diffraction-limited laser beam and imaging of the spatial profile of emission through confocal fluorescence detection enabled the read-out of exciton propagation lengths of up to 2 $\upmu$m in quasi-1D assemblies of LH2 complexes. To the best of our knowledge, theoretical models could only explain such a diffusion by ignoring static disorder and underestimating dephasing Jang et al. (2001); Pajusalu et al. (2015); Cupellini et al. (2018), resulting in long-range excitonic delocalization across approximately 40 pigments Lloyd and Mohseni (2010); Abasto et al. (2012); Yang et al. (2012), a value that is in conflict with the experimental observations Pullerits et al. (1996); Chachisvilis et al. (1997); Monshouwer et al. (1997); Trinkunas et al. (2001); Cheng and Silbey (2006).

In order to examine this experiment, we determine the rates of the Pauli master equations Eq. (II.1) for stochastic realizations of the pigment energies $\epsilon_{i}$ and relative orientations of a 1D array of $1001$ coplanar LH2 complexes as shown schematically in Fig. 3a, and study the stationary exciton distribution $\bar{p}^{n}_{\alpha}$. For moderate excitation power which allows us to remain in the single excitaton sector, driving can be modelled to take place on a single ring via incoherent B800-to-B850 energy transfer (Appendix F). Exciton distributions arising from a Gaussian laser profile can then be obtained by convolution of the driving profile with the solution for localized driving (Appendix A).

Despite the presence of static disorder and the multi-level structure of each unit cell, the stationary population distribution across the LH2 array is still characterized by an exponential distribution around the ring at which driving takes place (Fig. 3a). This allows us to characterize the population profile by a single parameter $l_{\text{diff}}$. When considering an initial Gaussian beam 400 nm wide (full width at half maximum), our model is able to reproduce the experimentally observed micron-range exciton propagation lengths if $l=6.5\text{ nm}$ (Fig. 3b), while natural distances of 8.0–8.5 nm yield a barely noticeable spread of the exciton density. This result suggests that LH2 packing density is a key factor in determining the spatial extent of energy transfer. With a distance of $l=6.5$ nm, a competition between thermalization and transfer leads to the establishment of a non-equilibrium steady state exciton population within antenna units $\bar{p}_{\alpha}$, that has a larger weight on high energy dark states than the thermal distribution, as shown in Fig. 3c. This is a clear signature that the non-equilibrium transfer across these arrays partially proceeds via high energy dark states, which, as explained above, rely on excitonic delocalization within each ring unit cell.

Despite being optically dark, these high-energy excitons are more sensitive to non-radiative decay than their low-energy bright counterpart, as shown in Fig. 2c. Therefore, energy transfer and decay dynamics have competing effects on the diffusion length: Local exciton-phonon interactions leading to exciton thermalization slow down the transfer of excitons by moving population away from fast propagating high-energy states, while at the same time granting them more time to propagate, since low-energy excitons are less sensitive to non-radiative decay. Whether this leads to a larger diffusion length or not, depends on the specifics of the system under consideration. In order to test this possibility, we artificially slow down the exciton thermalization timescale from 200 fs to 2 ps. This corresponds in practice to forcing the steady state exciton distribution in a given ring $\bar{p}_{\alpha}^{n}$ to be further away from thermal equilibrium $p_{\alpha}^{\text{th}}$. The shift of population towards the higher excitons results in a shorter ranged transfer (Fig. 3d, orange dots), revealing that these light-harvesting arrays operate in the generalized “dark” state protection regime. The relevance of such shelving mechanism is also confirmed by the behavior of the effective transfer and decay rates,

[TABLE]

(Note that the presence of static disorder in the array forces us to keep track of the unit cell index $n$ and to average over both forward and backwards transfer rates.) For faster relaxation, i.e. larger coupling to local phonons, both effective transfer and decay rates are significantly decreased, as captured by Fig. 3d.

As noticed recently Jang (2018), the transfer rate $\overline{W}$ benefits from the engagement of dark states at small inter-ring separations. Their participation allows for an increase in $\overline{W}$ that exceeds the one predicted by interactions only mediated by the collective dipoles $\mathbf{D}_{\alpha}^{(n)}$ and $\mathbf{D}_{\beta}^{(n+1)}$ of two neighboring B850 rings. The relevance of the finite size of unit cells naturally makes the transfer more dependent on their geometrical details and relative arrangement. Indeed, we notice for example that a systematic out-of-plane angle of just 5*∘* as observed in lipid-reconstituted membranes Walz et al. (1998); Scheuring et al. (2003), which slightly increases the distance between the closest pigments in neighboring rings, slows down the effective transfer rate and therefore decreases diffusion length (Appendix F). However, the apparently beneficial involvement of fast-propagating dark states in the dynamics is countered by their sensitivity to non-radiative decay. Thus, in any situation in which we are interested in increasing the spatial extent of energy transfer rather than its speed, we need to consider the presence of realistic (radiative and non-radiative) decay channels and design the energy transfer process accordingly. Closely packed LH2 arrays seem to naturally operate in a parameter regime that sustains a generalized “dark” state protection mechanism, where fast exciton thermalization causes a shielding against non-radiative decay.

Finally, we notice from Fig. 3d that our model produces energy transfer rates that are in excellent agreement with existing theoretical estimates obtained for larger inter-ring separations ($l\gtrsim 7.5$ nm) Jang (2018); Jang et al. (2015), matching the conditions typically observed in reconstituted and biological light-harvesting membranes Scheuring et al. (2003, 2004). In this regime, inter-complex energy transfer proceeds from a completely equilibrated excitonic manifold (Fig. 3c), far from the non-equilibrium regime in which “dark” state shelving becomes relevant. This suggests that the design guidelines discussed in this work, while relevant for tightly packed nano-engineered systems, might be of secondary importance for more sparsely assembled biological LH2 membranes.

III Conclusions

In conclusion, we have shown with the help of an analytically solvable model that room temperature excitation energy transfer can benefit from quantum dynamics within modular unit cells and demonstrated that the resulting design principles apply in photosynthetic membranes with realistic physiological parameters, as well as in nano-fabricated architectures. The resulting hybrid “quantum-classical” design can increase both speed and propagation range thanks to the participation of the dark states due to excitonic delocalization within unit cells. On the one hand, improved speed can be achieved when the packing density of unit cells is made sufficiently large for the coupling of individual pigments on different unit cells to benefit from contributions of high-lying dark states. Crucially, close packing allows for exciton populations to depart from a thermal distribution, increasing the overall diffusion rate via non-equilibrium energy transfer. On the other hand, the exciton propagation length is extended by intra-unit cell relaxation, biasing electronic populations towards low-energy excitons, which are less sensitive to non-radiative decay and therefore increase the overall time window over which energy transfer can take place. This does not only exemplify the beneficial interplay of quantum coherent dynamics and environmental noise Plenio and Huelga (2008); Caruso et al. (2009); Huelga and Plenio (2013) but also provides basic mechanisms that underpin the micron-range propagation of excitations observed in artificial arrays of LH2 photosynthetic complexes. These can be explained by the speed up of inter-complex transfer rates induced by dense packing of light-harvesting units and the protection from non-radiative decay provided by low-energy excitons. Although we do not expect this transfer mechanism to be at play in biological LH2 membranes due to their large inter-complex separations, it could be tested on nano-engineered platforms in experiments where the exciton diffusion length is measured for different light-harvesting arrays, prepared with different LH2 packing densities. Further corroborating evidence for the mechanism that we propose is derived from experimental observation of reduced exciton lifetimes in reconstituted LH2 membranes compared to isolated complexes Pflock et al. (2008). Our model already captures this trend correctly in terms of an increased population of the high-energy excitons, and thus could serve as a solid starting point for a more thorough quantitative analysis of this effect. While the multiscale model used in this work contains all the necessary elements to discuss general energy transfer strategies in light-harvesting arrays, further improvements could be achieved by employing more refined theoretical descriptions of the light-harvesting units Jang and Mennucci (2018). We plan to do so in the future by applying some recently developed numerical techniques, allowing for an exact treatment of the non-Markovian exciton dynamics Somoza et al. (2019).

The fact that the speed and spatial extent of energy transfer can be directly related to excitonic delocalization suggests the possibility of using partial delocalization restricted to single unit cells (due to the magnitude of noise in real room-temperature scenarios) as a resource to optimize the range of propagation of electronic excitations in technological applications with the goal of outperforming the already extremely high efficiency of natural photosynthesis. The general nature of these design principles hints that this energy transfer scheme might find applications in a broad class of excitonic materials, not limited to the specific architecture excplicitly discussed in this work, although it will be task of future research to probe its technological feasibility.

Acknowledgments

The authors are grateful to A. Block and N. van Hulst for useful discussions. This work was supported by the ERC Synergy grants BioQ and HyperQ, an IQST PhD fellowship and was made possible through the support of a grant from the John Templeton Foundation. The opinions expressed in this publication are those of the authors and do not necessarily reflect the views of the John Templeton Foundation. The authors acknowledge support by the state of Baden-Württemberg through bwHPC and the German Research Foundation (DFG) through grant no INST 40/575-1 FUGG (JUSTUS 2 cluster).

APPENDIX A EXACTLY SOLVABLE MODEL

In this appendix, we present the main steps leading to the analytical solution of the minimal model describing incoherent energy transfer across a linear array of dimerized unit cells, which locally support quantum delocalization. We start from Eq. (II.1) in the main text, i.e. the discrete diffusion equation of a linear array composed by dimerized unit cells, each hosting two levels, $|b_{n}\rangle$ and $|d_{n}\rangle$ (bright and dark), which can hop to the neighboring cells and are subject to intra-cell relaxation and fluorescence. Writing explicitly both the equations for the bright and dark components for local injection of excitations at site $n=0$, we obtain

[TABLE]

where $p_{b}^{n}$ ($p_{d}^{n}$) denotes the population of the bright (dark) state of the $n$-th dimer. The continuum limit is achieved by identifying $p_{\alpha}^{n}(t)/l=p_{\alpha}(x,t)$ ($\alpha=b,d$) with $x=nl$ and taking the inter-dimer separation $l$ to be vanishingly small, so that $x$ becomes a continuous variable. To simplify the discussion, we assume that the cross-rates are equal, i.e. $W_{db}=W_{bd}=w$, which is the case for the configuration assumed in the main text. If $W_{bd}\neq W_{db}$, a drift term in the diffusion equation is introduced, and the final exciton distribution is not going to be symmetric around the injection point $x=0$. Thus, we obtain two coupled diffusion equations

[TABLE]

Introducing the vectors $\mathsf{p}=(p_{b},p_{d})^{\mathsf{T}}$, $\mathsf{I}=(I_{b},I_{d})^{\mathsf{T}}$ and the matrices

[TABLE]

we can rewrite Eq. 16–17 more compactly as

[TABLE]

Rewriting Eq. (24) in terms of the Fourier transform $\hat{\mathsf{p}}(q,t)=\int\textrm{d}x\ e^{-iqx}\mathsf{p}(x,t)$ and considering the stationary state at $t\to\infty$, we obtain the solution as

[TABLE]

The solution can be brought to a much more transparent form if we assume that the decay described by $\Gamma_{\alpha}$ takes place on a much slower timescale than inter- and intra-dimer energy transfer. Within this approximation, which is typically satisfied in light-harvesting complexes, we have

[TABLE]

to leading order in $\Gamma_{\alpha}$, where $\overline{\Gamma}$ and $\bar{p}_{\alpha}$ were defined in Eq. (9) and Eq. (10). This simple form allows us to sum the series in Eq. (25) to

[TABLE]

where $\overline{W}$ was defined in Eq. (8). Taking the inverse Fourier transform leads to the final form of the solution

[TABLE]

which is the exciton distribution discussed in the main text. As one would expect, if inter-exciton conversion due to the cross-rates $R_{bd}$, $R_{db}$ and $w$ is much faster than the decay timescale, both bright and dark manifold have the same final distribution, differing only by the a normalization constant $\bar{p}_{\alpha}$. It is easy to check that we would get the same exponential distribution if we considered a diffusion process involving unit cells containing a single level rather than two, with transfer and dissipation rates $\overline{W}$ and $\overline{\Gamma}$. Therefore, a dimeric unit cell results in the same type of diffusion as a monomeric one, where the large-scale diffusion properties can be tuned by changing the local dimer parameters.

A completely analogous procedure leads to the solution of the case in which $W_{bd}\neq W_{db}$, which we do not present in detail here. The imbalance between these two rates does not change the definitions of $\bar{p}_{\alpha}$, $\overline{W}$ and $\overline{\Gamma}$ (which only depend on the average cross-rate $w=(W_{bd}+W_{db})/2$), but leads to a drift coefficient proportional to $\Delta=(W_{bd}-W_{bd})(\bar{p}_{b}-\bar{p}_{d})$. As a final result, the exciton distribution is still exponentially localized around the injection point $x=0$, but the diffusion lengths for $x<0$ and $x>0$ are different. Thus, the diffusion is not symmetric. As a possible application, this property might be used to design excitonic wires that are able to switch the dominant direction of diffusion by small changes of their unit-cell properties, for example implementing artificial photoprotection.

Although in the case of the exactly solvable model we consider a completely homogeneous system, it is possible to show that the presence of moderate static disorder does not modify the exponential shape of the stationary exciton distribution, but only leads to a reduction of the diffusion length. This explains why, also in the case of the LH2 array presented in the main text, we still obtain an exponential distribution. To grasp of the effects of static disorder, we consider a linear array of monomeric unit cells with some slight random inhomogeneity in the transfer rates. The transfer rate from site $n$ to site $n+1$ (and vice versa) is given by $W_{n}=\overline{W}+\delta W_{n}$, where we assume $\langle\delta W_{n}\rangle=0$ and $\langle\delta W_{n}\delta W_{m}\rangle=\sigma^{2}\delta_{nm}$. The average population distribution $\langle p(x)\rangle$ for small disorder can be obtained by expanding the solution in the form of Eq. (25) to second order in $\delta W_{n}$ and taking the ensemble average. This leads to an exponential average exciton distribution with diffusion length $l_{\text{diff}}^{\prime}=l_{\text{diff}}(1-\sigma^{2}/2\overline{W}^{2})$, independent of the specific distribution of the rate fluctuations $\delta W_{n}$.

To conclude this section, we note that the stationary solution $\mathsf{p}(x,\infty)$ for local exciton injection $\delta(x)$ is sufficient to determine the exciton profile $\mathsf{p}^{\prime}(x,\infty)$ for any other (normalized) injection profile $g(x)$. In fact, defining $\hat{g}(q)$ the Fourier transform of $g(x)$, we immediately obtain that the new solution satisfies

[TABLE]

Transforming back to real space, we obtain the new solution for generic driving as a convolution between the solution for local driving and the generic driving profile, namely

[TABLE]

The same principle allows us to draw conclusions for exciton propagation in a LH2 array upon driving with a Gaussian laser profile, using only the results of simulations for local driving.

APPENDIX B MICROSCOPIC ORIGIN OF NON-RADIATIVE DECAY

In this appendix, we derive a microscopic Hamiltonian describing vibrationally-induced non-radiative decay of a photoexcited molecular state. We start from the ab initio Hamiltonian of a single chromophore. We determine the form of the non-adiabatic coupling at the basis of nonradiative decay in the usual harmonic approximation for intra-molecular vibrations. We use this result to compute the internal conversion rate for delocalized excitonic states of a molecular aggregate.

B.1 Molecular Hamiltonian

The full Hamiltonian of a molecule is given by May and Kühn (2004)

[TABLE]

where $T$ denotes the kinetic energy and $V$ the Coulomb interactions. The sets of electronic coordinates and momenta are denoted by ${\bf r}$ and ${\bf p}$, whereas the nuclear degrees of freedom are described by ${\bf R}$ and ${\bf P}$. The typical challenge in condensed matter physics is to find approximate eigenstates and eigenvalues of this Hamiltonian. The huge difference between electronic and nuclear masses justifies the typical Born-Oppenheimer (BO) approach, where the electronic degrees of freedom are treated on a fully quantum level (i.e. $\mathbf{r}\to\hat{\mathbf{r}}$ and $\mathbf{p}\to\hat{\mathbf{p}}$), whereas the nuclear coordinates are kept fixed and their momenta are initially neglected. This allows to find the adiabatic electronic eigenstates by diagonalization of the electronic Hamiltonian at fixed $\bf R$:

[TABLE]

The eigenstates $|\phi^{i}_{\mathbf{R}}\rangle$ and eigenenergies $\varepsilon^{i}_{\bf R}$—which take the name of potential energy surfaces (PES)—depend parametrically on the nuclear coordinates. Now we express the Hamiltonian (31) in terms of these adiabatic electronic states. Although their interpretation as physically relevant states depends on whether the nuclear kinetic energy is negligible or not, they always form a perfectly legitimate orthonormal basis of the electronic Hilbert space at fixed $\bf R$. Making use of (32), we see that

[TABLE]

Note that the second term on the right-hand side does not necessarily reduce to $\delta^{ij}T_{\text{nuc}}({\bf P})$, due to the parametric dependence on $\bf R$ of the electronic states. However, if nuclei are considered as classical particles, the nuclear kinetic energy does not affect the electronic adiabatic wavefunction in any way, and the full molecular Hamiltonian (31) takes the well-known BO form

[TABLE]

The adiabaticity of this Hamiltonian is represented by the fact that the nuclear motion does not couple different electronic states: when the electrons are in the state $|\phi^{i}_{\mathbf{R}}\rangle$, the nuclei evolve according to the effective Hamiltonian $h_{i}(\mathbf{R},\mathbf{P})=T_{\text{nuc}}({\bf P})+\varepsilon^{i}_{\bf R}$.

Introducing a diagonal nuclear mass tensor $\mathsf{M}$ and expanding the $i$-th PES for small deviations from its stable equilibrium configuration ${\bf R}^{i}$ (i.e. close to its minimum $\varepsilon^{i}_{{\bf R}^{i}}$), the nuclear Hamiltonian $h_{i}$ becomes

[TABLE]

where $\mathsf{H}\geq 0$ is the Hessian matrix of the PES $\varepsilon^{i}_{\bf R}$ at ${\bf R}^{i}$. On the second line, we introduced the mass-rescaled normal coordinates and momenta ${\bm{\xi}}=\mathsf{U}\mathsf{M}^{1/2}{\bf R}$ and ${\bm{\pi}}=\mathsf{U}\mathsf{M}^{-1/2}{\bf P}$, where $\mathsf{U}$ is the unitary transformation that diagonalizes the mass-rescaled Hessian, i.e. $\mathsf{U}^{\mathsf{T}}{}\mathsf{D}\mathsf{U}=\mathsf{M}^{-1/2}\mathsf{H}\mathsf{M}^{-1/2}$. $\mathsf{D}$ has only diagonal entries, corresponding to the square of the normal frequencies $\omega_{k}^{2}$. As usual, we have neglected the Duschinsky mixing, that is, we have assumed that the Hessian $\mathsf{H}$ at the nuclear equilibrium configuration is independent of the specific PES, so that it is diagonalized by the same transformation and yields the same eigenvalues (i.e. vibrational frequencies) for the PESs in which we are interested.

At this point, we can easily quantize the nuclear normal modes by redefining them in terms of a set of bosonic annihilation (creation) operators $\hat{b}_{k}$ ($\hat{b}_{k}^{\dagger}$). Assuming from now on $\hbar=1$, we have

[TABLE]

The quantized vibrational Hamiltonian (B.1) becomes

[TABLE]

Note that we have introduced the Huang-Rhys factors $s_{k}^{i}={\xi_{k}^{i}}^{2}\omega_{k}/2$, which in general depend on the PES under consideration.

Since our ultimate goal is to describe the dynamics of photo-excitations, we focus on two electronic states $|\phi^{g}_{\bf R}\rangle$ and $|\phi^{e}_{\bf R}\rangle$, describing the electronic ground state and the first optically excited state. We neglect any further dependence of these states on the nuclear coordinates and refer to them simply as $|\phi_{g}\rangle$ and $|\phi_{e}\rangle$. If we refer all energies and nuclear coordinates to the minimum of the PES of the electronic ground state, the BO Hamiltonian (34) takes the usual spin-boson form with pure dephasing interaction

[TABLE]

where $\varepsilon=\varepsilon^{e}_{{\bf R}^{e}}+\sum_{k}s_{k}\omega_{k}$ and $g_{k}=-\omega_{k}\sqrt{s^{e}_{k}}$. The adiabatic electron-phonon coupling can then be described by the spectral density

[TABLE]

B.2 Non-adiabatic electron-phonon coupling

Now, we turn back to (33) and consider the effect that the nuclear kinetic energy has on the adiabatic electronic wavefunction. For this, we take into account the quantum nature of nuclei as well, i.e. considering ${\bf R}\rightarrow\hat{\bf R}$ and ${\bf P}\rightarrow\hat{\bf P}$ as quantum mechanical operators. Note that the second term on the right-hand side of (33) acts non-trivially on both electronic and nuclear degrees of freedom. To better understand this, we consider how the nuclear kinetic energy operator acts on a product state formed by $|\phi^{j}_{\mathbf{R}}\rangle$ and a nuclear state $|\chi_{\nu}\rangle$. In coordinate representation, $\hat{\bf P}$ becomes the differential operator $-i\partial_{\bf R}$, therefore we have

[TABLE]

Using (41), we can now compute the matrix element of the nuclear kinetic energy operator between two arbitrary electronic-vibrational states as

[TABLE]

where the three terms appearing on the second line of (42) directly originate from those appearing in (41). It is easy to spot the first of them as the one giving rise to the BO Hamiltonian (34). On the other hand, the terms involving derivatives of the adiabatic electronic wavefunction with respect to the nucrear coordinates are not necessarily proportional to $\delta^{ij}$, and therefore describe the non-adiabatic coupling between different PESs. Within the BO framework, they are neglected by merit of the argument that the electronic wavefunctions depend very weakly on the nuclear coordinates, therefore their derivatives are negligible. In fact, electronic wavefunctions usually vary significantly over distances of 1–10 Å, corresponding to typical inter-nuclear separations. On the other hand, typical nuclei in organic molecules jiggle by about 0.01–0.1 Å around their equilibrium position. On this length-scale, therefore, electronic wavefunctions are expected to be fairly smooth and almost constant. Thus, the last two terms in (42) can be regarded as higher order contributions of a perturbative expansion whose zeroth order term corresponds to the BO Hamiltonian.

For simplicity, let us focus only on the lowest non-zero order of this non-adiabatic perturbation, i.e. the one involving the first derivative of the adiabatic electronic wavefunction. Let us define the quantity $A^{ij}_{k}({\bf R}):=-i\langle\phi^{i}_{\bf R}|\partial_{R_{k}}\phi^{j}_{\bf R}\rangle$, which is closely related to the Berry connection. Note that we can always assume the phase of the adiabatic electronic wavefunctions to be independent of ${\bf R}$, corresponding to a specific gauge choice for the Berry connection. This, together with the normalization condition $\langle\phi^{i}_{\bf R}|\phi^{i}_{\bf R}\rangle=1$, immediately implies that $F^{ii}_{k}({\bf R})=0$, meaning that this term only connects different electronic states and does not modify the vibrational Hamiltonian within a given PES. On top of that, the orthogonality condition $\langle\phi^{i}_{\bf R}|\phi^{j}_{\bf R}\rangle=0$ further imposes a hermiticity constraint ${A^{ij}_{k}({\bf R})^{*}}=A^{ji}_{k}({\bf R})$. As usually assumed for optical transitions, we neglect the dependence of the coupling between different PESs on the nuclear configuration. This means that $A^{ij}_{k}({\bf R})\approx A^{ij}_{k}({\bf R}^{i})$ is practically constant in the integral over $\bf R$ in (42). We also neglect any further dependence of the adiabatic states on the nuclear coordinates, as we did in the previous section in order to recast the BO Hamiltonian (34) into spin-boson form (39). Thus, we can finally write the initial molecular Hamiltonian (31) including the lowest order non-adiabatic correction as $\hat{H}_{\text{BO}}+\hat{H}_{\text{non-ad}}$, where

[TABLE]

where in the last step we have rewritten the nuclear momenta in terms of the normal modes defined above and introduced the transformation ${\bm{\alpha}}^{ij}=\mathsf{U}\mathsf{M}^{-1/2}\mathbf{A}^{ij}({\bf R}^{i})$ on the coefficients $A^{ij}_{k}({\bf R}^{i})$.

Thus, the non-adiabatic electron-phonon coupling can lead to transitions between different electronic states $|\phi_{i}\rangle$ and $|\phi_{j}\rangle$ initiated by the action of nuclear momenta. If we consider only two electronic states, namely the ground state $g$ and one optically excited state $e$, we can see this as causing non-radiative transitions between $g$ and $e$ mediated by the coupling constant $f_{k}:=\alpha^{ge}_{k}\sqrt{\omega_{k}/2}$. This process goes under the name of internal conversion (IC), and is therefore determined by the following spectral density

[TABLE]

analogous to (40).

B.3 Internal conversion in molecular aggregates

While it is generally true that the non-adiabatic couplings are much weaker than the pure-dephasing electron-phonon coupling, they still play a role on longer time-scales, where other phenomena such as fluorescence enter into play. Therefore we should devise a way to treat IC on the same footing as fluorescence in molecular aggregates. We consider an aggregate of $N$ interacting two-level chromophores, each one described by the Hamiltonian (43). We focus on the subspace spanned by the global ground state $|g\rangle=|\phi_{g}\rangle_{1}\dots|\phi_{g}\rangle_{N}$ and the single excitations $|i\rangle=|\phi_{g}\rangle_{1}\dots|\phi_{e}\rangle_{i}\dots|\phi_{g}\rangle_{N}$. Assuming that the only site dependent parameters of our model are the electronic excitation energies $\varepsilon_{i}$ and couplings $J_{ij}$, the Hamiltonian of the aggregate is

[TABLE]

On the first line, we recognize the free Hamiltonian of the excitonic system (diagonalized by excitons $|\alpha\rangle=\sum_{i}|i\rangle\langle i|\alpha\rangle$ with energy $E_{\alpha}$) and the free Hamiltonian of environmental vibrations. On the second and third line, we find the system-bath interactions, mediated by the two bath operators

[TABLE]

which respectively cause pure dephasing in the site basis and non-radiative transitions between site excitations and global ground state. Postponing a more detailed justification for later, we allow for correlations between different vibrational environments to be present when the bath is at equilibrium,

[TABLE]

where the time evolution is computed with respect to the free bath Hamiltonian and the expectation value $\langle\cdot\rangle_{\text{eq}}$ is taken on the stationary state of the bath. The correlation functions $\mathcal{G}(t)$ and $\mathcal{F}(t)$ are assumed to be site-independent, and only depend on the temperature of the vibrational bath and on the spectral densities $\mathcal{J}(\omega)$ and $\mathcal{J}_{\text{non-ad}}(\omega)$. We assume that the internal timescale of the bath is sufficiently fast and the system-bath coupling is sufficiently weak so that the dynamics of electronic excitations can be described by a Lindblad equation in the exciton basis.

We focus on internal conversion for now. Following the microscopic derivation outlined in Rivas and Huelga (2011), we obtain a Lindblad equation with jump operators $|g\rangle\langle\alpha|$ and rates

[TABLE]

appearing in Eq. (5). The rate $\Gamma_{\text{non-rad}}$ is the single-pigment non-radiative rate, determined by

[TABLE]

In principle, it depends on the frequency at which the specific transition to the ground state takes place (i.e. $\omega=E_{\alpha}$ for $|g\rangle\langle\alpha|$). However, since optical frequencies are much higher than vibrational and thermal energies, we can assume $\mathcal{J}_{\text{non-ad}}(\omega)$ to be fairly small and constant across the excitonic energies $E_{\alpha}$.

The correlations introduced in (48) and (49) are crucial for the determination of the distribution of the IC rate across the excitonic manifold, as shown in Fig. 4. For example, if the correlations are absent, i.e. $\kappa_{ij}=\delta_{ij}$, Eq. (51) predicts that the IC rate is essentially the same for all excitons (blue line). If, instead, we allow for some positive correlations decaying exponentially as a function of the distance between chromophores with some characteristic correlation length $r_{\text{c}}$ (i,e, $\kappa_{ij}=e^{-r_{ij}/r_{\text{c}}}$), excitons which delocalize over neighboring pigments show an increased decay rate (green and yellow solid lines). Correlations of such form are commonly used when modelling various optical spectra of pigment-protein complexes Renger and Marcus (2002). If the correlations are negative for neighbouring pigments, the opposite behavior is observed (dashed lines). Optical correlations between transition dipole moments give the usual fluorescence profile (red line), with few bright excitons at the bottom of the band. We will argue later in favour of the correlations between local vibrational environments that we have postulated, by considering some experimental results.

B.4 Other non-radiative decay pathways

So far, we have seen how IC can be described microscopically and how it influences the non-radiative decay of different optical excitations in a molecular aggregate. However, there could be other processes that compete with fluorescence and IC on the same timescale, which could influence the effective decay rate from a given exciton. The other decay channel typically present in molecular systems is inter-system crossing (ISC) Marian (2012). During this process, a singlet excitation $S_{1}$ generated by optical absorption from the singlet ground state $S_{0}$ can be turned into an excited triplet $T_{1}$ by means of spin-orbit interactions. Since the radiative transition from $T_{1}$ to $S_{0}$ is strongly suppressed, the excitation is lost non-radiatively either by IC or by triplet-triplet energy transfer to other pigments present in the molecular aggregates (i.e. carotenoids), before it can flip the spin of a neighboring oxygen molecule, thus forming singlet oxygen, which can cause oxidative damage to the photosynthetic apparatus.

If our focus is to follow the dynamics of excitons in a narrow band (i.e. the B850 of the LH2 complex), we do not want to take into account all these processes: we are only interested in the rate at which spin-orbit interactions transfer excitations from one adiabatic PES to another. As we explicitly worked out for IC, it has been shown that also the spin-orbit interaction results in coupling between adiabatic PESs mediated by the nuclear momentum Englman and Jortner (1970). This tells us that the effective Hamiltonian causing ISC will exhibit couplings of the same form as (43), with properly redefined coupling constants. The resulting ISC rate, therefore, has the same form as (51), i.e. it depends on the density of vibrational states at excitonic energies. Assuming that this density of states is sufficiently flat, as done before, we have also in this case that the dependence of the rate on the excitonic state is dominated by the pattern of correlations between local baths. ISC can be therefore absorbed together with IC into a single rate $\kappa_{\alpha}\Gamma_{\text{non-rad}}$, describing non-radiative decay of population from exciton $\alpha$.

APPENDIX C PHYSICAL ORIGIN OF CORRELATED ENVIRONMENTS

In this appendix we argue in favor of the presence of correlations between the vibrational environments of single pigments belonging to the same pigment-protein complex, which take place on the timescale of non-radiative decay. First, we propose an estimate based on experimental results, and later move on to discuss a physical mechanism which can support inter-site vibrational correlations, relying on inter-molecular vibrational modes.

C.1 Phenomenological view

Let us step back for a moment from the microscopic derivation of internal conversion, and follow another approach. We start from experimental evidence, and use this knowledge to set up a phenomenological model of nonradiative decay of a molecular aggregate. The most valuable insight comes from comparing excited state lifetimes ($\tau$) and fluorescence quantum yields ($\phi$) of isolated pigments with those of the aggregate. Table 1 shows theses quantities in the case of bacteriochlorophyll a molecules (Bchl-a) and LH2 complexes Monshouwer et al. (1997). Since the fluorescence quantum yield of LH2 is lower than the one of Bchl-a, there must be some additional decay channels in LH2 that are not present in the monomers and which can therefore be interpreted as an effect of aggregation.

To be more quantitative, let us define the quantum yield of fluorescence

[TABLE]

and express the radiative and nonradiative decay rates as $\Gamma_{\text{rad}}=\phi/\tau$ and $\Gamma_{\text{non-rad}}=(1-\phi)/\tau$. The values reported in Monshouwer et al. (1997) allow us to determine radiative and nonradiative decay rates for Bchl-a and LH2 (Table 1). While the difference in radiative lifetimes between the monomer and the aggregate is easily explained by exciton delocalization and superradiance, a clear molecular mechanism underpinning the mismatch between nonradiative decay rates is not known with certainty. Nonetheless, we can put forward a simple argument. Since the nonradiative decay in LH2 is about 3.5 times faster than in Bchl, we can imagine the existence of additional dissipation channels. In fact, we expect the same intra-molecular decay channels present in Bchl to be at play also in LH2. Thus, the additional dissipation present in LH2 must come from some other decay pathway that has no analogue in Bchl, which must account for $(\Gamma_{\text{non-rad}}^{\text{LH2}}-\Gamma_{\text{non-rad}}^{\text{Bchl}})/\Gamma_{\text{non-rad}}^{\text{LH2}}=71\%$ of the nonradiative dissipation. The first reasonable candidate is the protein environment, which can offer other IC pathways through vibrations and conformational changes. Since the B850 ring of LH2 is composed by dimeric subunits bound to the same protein, it makes sense to assume that the IC channels offered by the protein can in principle be correlated. Following this argument, we allow for correlations between the vibrational environments that couple to IC and ISC transitions. As a consequence, excitons that delocalize on neighboring pigments can experience a modified decay rate, as discussed above and shown in Fig. 4.

C.2 Microscopic mechanism

We have seen how, within a BO framework, a coupling arises between (optical) electronic transitions and vibrations of a molecule. A sudden change in the electronic state leads to a different electrostatic potential experienced by the nuclei, which therefore initiate their dynamics. For this reason, this coupling involves only vibrations of the nuclei over which the electronic states are delocalized. It is therefore clear that intra-molecular modes will experience this type of direct coupling to the electronic dynamics.

However, if the chromophore is embedded in a protein environment, the vibrational modes of the protein will also influence electronic energies, leading to a direct coupling between electronic excitations of the chromophore and longer wavelength vibrational modes. This observation lies at the heart of theoretical descriptions of electronic resonances coupled to a shared phonon environment Soules and Duke (1971). Thus, the pure dephasing coupling in (B.3) can be written explicitly in terms of these protein vibrational modes $\hat{c}_{q}$ with frequency $\nu_{q}$ as

[TABLE]

By redefining the operators $\hat{G}_{i}$ (46) to include also the protein modes $\hat{c}_{q}$, we obtain the following bath correlation functions

[TABLE]

While the term on the first line (corresponding to intra-molecular modes) vanishes for for $i\neq j$, the second term (arising from protein modes) is able to generate correlations between sites. Moreover, since protein motion takes place on a larger and slower scale than intra-molecular vibrations, it is reasonable to assume that two neighboring pigments $i$ and $j$ couple to the protein motion with the same phase ($\xi_{iq}\xi_{jq}>0$). This results in positive inter-site correlations, and thus can lead to the redistribution of non-radiative decay rates discussed in Appendix B.

Another possible mechanism able to generate positive inter-site vibrational correlations is based on the idea that the vibrational motion of a chromophore can mechanically couple to the slow vibrations of the protein environment. To formalize better this idea, let us consider a slightly modified version of the model described by the Hamiltonion (B.3). To simplify the discussion, we consider a molecular aggregate composed only of two pigments, i.e. $i=1,2$, and we neglect the direct coupling between electronic excitations and protein modes (52). If the two pigments are bound to the same protein, it makes sense to assume that their intra-molecular modes $\hat{b}_{1,k}$ and $\hat{b}_{2,k}$ couple to a set of common modes $\hat{c}_{q}$ with frequency $\nu_{q}$. Since these shared modes are associated to the protein structure, it is reasonable to assume that they are much slower and effectively classical, meaning that their energy $\nu_{q}$ is much smaller than thermal energy $k_{\text{B}}T$. We assume that the intra-molecular modes of the two pigments couple the protein vibrations with the same strength and phase. This choice makes sense if we think that the two pigments have opposite orientation and are embedded in an elastic medium (see Fig. 5). The Hamiltonian of the vibrational bath then reads

[TABLE]

where we have neglected the counter-rotating coupling terms to simplify the following treatment, although they could be included. Defining the symmetric and antisymmetric combinations of intra-molecular modes $\hat{b}_{k}=(\hat{b}_{1,k}+\hat{b}_{2,k})/\sqrt{2}$ and ${\hat{b}^{-}}_{k}=(\hat{b}_{1,k}-\hat{b}_{2,k})/\sqrt{2}$, we can rewrite the coupling as $\sum_{k,q}\sqrt{2}\eta_{kq}\hat{b}_{k}\hat{c}_{q}^{\dagger}+\text{ H.c.}$, i.e. the antisymmetric modes ${\hat{b}^{-}}_{k}$ are insensitive to the coupling to the protein. The symmetric modes, on the other hand, mix with the protein vibrations, yielding the new set of normal modes

[TABLE]

with frequencies $\tilde{\omega}_{k}$ and $\tilde{\nu}_{q}$. Thanks to its unitarity, the transformation can be easily inverted and the bath operators $\hat{G}_{i}$ defined in (46) can be reexpressed as

[TABLE]

where we have defined the transformed system-bath couplings

[TABLE]

Analogous expressions can be found for the operators $\hat{F}_{i}$ defined in (47).

Let us now focus on the bath correlation functions for the pure dephasing environment introduced in Eq. (48), evaluated on the thermal state $\rho_{\text{eq}}\propto e^{-\hat{H}_{\text{B}}/k_{\text{B}}T}$. (The system-bath coupling leading to internal conversion can be treated in a completely analogous way.) Using the new normal mode decomposition of $\hat{G}_{i}$ (57), we get

[TABLE]

To lowest order in the inter-site vibrational coupling $\eta_{kq}$, we find the following perturbative expressions

[TABLE]

Note that, in order to be perfectly consistent, we should keep also the second-order correction to $\mathsf{G}_{k}$. However, since the protein modes generally have much smaller energies $\nu_{q}$ than both $k_{\text{B}}T$ and intra-molecular modes $\omega_{k}$, their thermal occupation number is much higher, i.e. $n(\tilde{\nu}_{q})\gg n(\tilde{\omega}_{k})$. Therefore, the correction arising from the first line in (59) is negligible with respect to the one originating from the second line. Taking into account all these assumptions, we obtain

[TABLE]

Note that both these equations can be rewritten in terms of positive spectral densities: In the case of (62) it is the one defined in (40), whereas for (63) we can define

[TABLE]

Taking the Fourier transform of these correlation functions, we obtain the rates that allow us to write down a Lindblad equation for the reduced dynamics of the excitons, which are proportional to the spectral densities $\mathcal{J}(\omega)$ and $\Delta\mathcal{J}(\omega)$.

At this point, we can see that our particular choice of coupling the local intra-molecular modes to the common protein modes with the same phase leads to the establishment of positive inter-site correlations ($\Delta\mathcal{J}(\omega)>0$). If we chose to couple local modes with opposite phase instead, we would end up with negative correlations between the vibrational environments.

APPENDIX D LINESHAPES, LINEAR SPECTRA AND AND TRANSFER RATES

In this section we determine the lineshapes of an excitonic system (i.e. our unit cells), which allow for the calculation of both linear optical spectra and excitonic transfer rates. Consider an excitonic system with Hamiltonian $\hat{H}=\sum_{\alpha}E_{\alpha}|\alpha\rangle\langle\alpha|$, which we identify as our unit cell, where Greek letters denote exciton states of the unit cell $|\alpha\rangle=\sum_{i}|i\rangle\langle i|\alpha\rangle$. (Here we relax the convention used in the main text, according to which primed indices refer to the same subunit.) Within a single unit cell, excited state population thermalizes due to the interaction between electronic and vibrational degrees of freedom of both intra- and inter-molecular origin. The main electron-phonon coupling mechanism is described by the pure-dephasing interaction presented in the aggregate Hamiltonian (B.3). Under the conditions of weak electron-phonon coupling and fast vibrational relaxation, this process can be described by a Lindblad equation. In the presence of correlations $\kappa_{ij}$ between different local environments, we can write down the resulting Lindblad dissipator as

[TABLE]

where the first line describes population transfer across different excitons and the second describes pure dephasing processes. The rates can be obtained respectively as

[TABLE]

where $\omega_{\alpha\beta}=E_{\alpha}-E_{\beta}$, and $\gamma$ is the single pigment optical dephasing rate. Note that, when the local vibrational environments are not correlated, i.e. $\kappa_{ij}=\delta_{ij}$, the factors involving excitonic amplitudes in (66) and (67) reduce to the spatial overlap between excitons and their participation ratio. Lindblad dynamics in the exciton basis results in an exponential decay of the optical coherences between exciton $\alpha$ and the ground state $g$ with a rate $\gamma_{\alpha}=\gamma\zeta_{\alpha\alpha}+\sum_{\beta(\neq\alpha)}R_{\beta\alpha}/2$, whereas the inter-exciton coherence decays with a rate $\gamma_{\alpha\beta}=\gamma_{\alpha}+\gamma_{\beta}-2\gamma\zeta_{\alpha\beta}$. This leads to simple expressions when computing optical absorption and emission spectra. These are related to the absorption and emission tensors, $\mathsf{A}(\omega)$ and $\mathsf{E}(\omega)$ respectively, defined by their matrix elements

[TABLE]

Here $\hat{\sigma}_{\mu}(t)$ denotes the Heisenberg time evolution of the annihilation operator of exciton $\mu$, $\hat{\sigma}_{\mu}=|g\rangle\langle\mu|$, and $\langle\cdot\rangle_{g(e)}$ denotes the average over the equilibrium ground (excited) state. Evolving the transition dipole operator through the dual of $\mathcal{D}$, results in the simple expression $\hat{\sigma}_{\mu}(t)=|g\rangle\langle\mu|e^{-(\gamma_{\mu}+iE_{\mu})t}$, which leads to diagonal absorption and emission tensors, with each exciton having a Lorentzian lineshape, i.e. $\mathsf{A}_{\alpha\beta}(\omega)=\delta_{\alpha\beta}f_{\alpha}(\omega)$ and $\mathsf{E}_{\alpha\beta}(\omega)=\delta_{\alpha\beta}f_{\alpha}(\omega)p^{\textrm{th}}_{\alpha}$, where $p^{\textrm{th}}_{\alpha}$ is the thermal population of exciton $\alpha$, and

[TABLE]

Once we have the lineshape for each exciton, we weight each individual lineshape by its associated brightness $\chi_{\alpha}\approx|\mathbf{D}_{\alpha}|^{2}/d^{2}$, and calculate straightforwardly the absorption spectrum of the unit cell as $\omega\sum_{\alpha}\chi_{\alpha}f_{\alpha}(\omega)$. In the case of the B850 subunit, once averaged over static disorder, this expression gives excellent agreement with experimental results, as shown in the main text (Fig. 2b).

If we introduce a second excitonic system, weakly interacting (with respect to the timescales associated to dephasing) with the first one via dipole-dipole couplings

[TABLE]

we can calculate the incoherent energy transfer rate between excitons $\alpha$ and $\alpha^{\prime}$ on the two subunits via generalized Förster theory, according to which, we have

[TABLE]

where $V_{\alpha\alpha^{\prime}}=\sum_{i,i^{\prime}}\langle\alpha|i\rangle V_{ii^{\prime}}\langle i^{\prime}|\alpha^{\prime}\rangle$, which reduces to Eq. (4) from the main text in the case of dimeric unit cells, once $\hbar$ is reintroduced.

APPENDIX E PARAMETRIZATION OF THE B850 RING

Fig. 6 and Table 2 summarize the geometry and the parameters that have been used in the simulations of the B850 rings. The ring structure is dimerized, meaning that each pigment $i$ is identified by two indices, one specifying the dimer $(n=1,\dots,9)$ and the other the position within the dimer $(\nu=1,2)$. The two bacteriochlorophyll (Bchl) molecules belonging to the same dimer are usually labelled $\alpha$ and $\beta$. A single ring is described by the Hamiltonian

[TABLE]

where the primed sum indicates summation over all couples of non-adjacent pigments and $J_{ij}$ stands in this case for the dipole-dipole interaction between pigments. The fluctuations of the site energies $\delta\varepsilon_{i}$ are given by the sum of two Gaussian random variables: one with standard deviation $\sigma_{\text{p}}$, completely uncorrelated for different pigments, describing local energy shifts due to the slightly different protein environments; one with standard deviation $\sigma_{0}$, which is the same for all pigments, describing global shifts of the ground state energy.

Throughout this work, the B800 ring enters the dynamics only during the excitation process as we describe more in detail in the next section. However, we explicitly included B800 rings in the calculation of the absorption spectra shown in Fig. 2b, for a better comparison with the experimental data. We model the B800 ring as a set of $9$ transition dipoles arranged on a circle of radius $3.1$ nm, concentric to the B850 ring and vertically displaced from it by $1.7$ nm. We consider the B800 dipoles to be perfectly coplanar and tangent to the circle. In accordance with previous works Caycedo-Soler et al. (2018), we set the ratio between B850 to B800 pigment dipole strength to $1.1$. The site energy of B800 pigments $\varepsilon_{\text{B800}}$ has an average value of $1.25\ 10^{4}$ cm*-1*, with Gaussian static disorder with standard deviation $\sigma_{\text{B800}}=100$ cm*-1*, whereas the single pigment optical dephasing rate is $\gamma_{\text{B800}}=70$ cm*-1* Caycedo-Soler et al. (2018).

The B800 pigments couple to each other and to the B850 pigments via dipole-dipole interactions. As a result, one may expect that excitons can in principle delocalize across both rings. However, the presence of static disorder is sufficient to destroy any inter-ring delocalization Caycedo-Soler et al. (2018). This can be seen in Fig. 7, where we plot the average population of a LH2 exciton $|\tilde{\alpha}\rangle$ (diagonalizing the total B800 + B850 Hamiltonian) on site $|i\rangle$. The absence of significant B800-B850 coherent mixing thus justifies the approach adopted throughout the paper, where we only consider the indirect effect of incoherent B800-to-B850 energy transfer to populate the B850 manifold after initial laser excitation.

APPENDIX F SIMULATIONS OF A LH2 ARRAY

In this section we give more details on the simulations of the full LH2 linear array. In analogy with the experiment of Escalante et al. Escalante et al. (2010), we look at the stationary exciton probability profile along the linear array upon continuous wave driving, i.e. we look numerically for the steady state of Eq. (II.1) (with all rates dependent also on the specific subunit $n$, since every subunit exhibits a different realization of static disorder, i.e. different spectra and pigment positions, affecting relaxation, fluorescence, transfer and injection rates). In the experiment the driving is provided by a 800 nm laser with spatial intensity profile with a full width at half maximum of $400\pm 50$ nm. In order to simplify the numerics, we consider instead injection on a single B850 complex at the center of the chain. In this way, we can simulate shorter chains (up to 1001 subunits) and be safely protected from systematic errors introduced by the finite size of the array. The results for spatially broad excitation can be recovered as shown in Appendix A in the case of the exactly solvable model, i.e. by convolving the result for local injection with the desired excitation profile.

We clarified that our injection profile is local, i.e. $I^{n}_{\alpha}\propto\delta_{n_{0}n}$ where $n_{0}=501$ is the central site of the chain, but we have not specified yet in which states $\alpha$ of the ring excitations are injected. Optical excitation at 800 nm cannot be absorbed by the B850 subunit (as seen in Fig. 2b), nevertheless it can excite the B800 ring. We are not explicitly considering the B800 subunit in our model. However, excitations enter the B850 ring upon downhill B800 $\to$ B850 energy transfer. Due to the small coherent coupling between the two concentric rings, this transfer process is largely controlled by incoherent rates of the form of Eq. (D), therefore we can think of an indirect excitation of the B850 excitons with an energy distribution given by

[TABLE]

where

$\varepsilon_{\textrm{B800}}$ and $\gamma_{\textrm{B800}}$ are optical gap and dephasing rate of the B800 pigments, whose values are given in Appendix E.

The disorder-averaged injection profile (73) is shown in Fig. 8a.

The resulting stationary probability profiles shown in Fig. 8b for different lattice steps $l$ ranging from 8.5 nm to 6.5 nm, after averaging over $10^{3}$ realizations of static disorder. The probabilities $\bar{p}^{n}_{\alpha}$, shown as a function of the position $x=(n-n_{0})l$, exhibit a clear exponential decay around the injection site in the middle of the array. A convolution of these exciton population profiles with a Gaussian injection profile with 400 nm full width at half maximum yields the distributions reported in Fig. 3b.

All simulations are performed for coplanar arrangements of B850 rings. As discussed in the main text, the nano-engineered arrays analyzed here allow for coplanar B850 rings. However, in the main text we also note that in biological light-harvesting membranes, an angle of about 5∘ relative to the aggregation plane is observed. This small tilt leads to slower energy transfer, especially for short inter-ring distances, as shown in Fig. 9.

Lastly, we note that all simulated LH2 arrays feature a fixed inter-complex distance $l$ throughout the length of the array. This clearly represents an approximation, since every realistic macromolecular assembly will show fluctuations in the lattice constant. While these large-scale geometric defects are expected to severely limit the range of ballistic energy transfer and eventually lead to diffusive transport, they only act as an additional source of static disorder in our model, which already describes energy transfer in the diffusive regime. Thus, in our case, this further noise source will lead to changes in the transfer rates but will not affect the scaling of diffusion length with respect to the lattice constant, leaving our main conclusions unaltered. Nevertheless, the effect of lattice constant fluctuations can be estimated by analyzing the dependence of the diffusion length $l_{\text{diff}}$ on the inter-complex distance $l$, shown in Fig. 3d. This seems to be described by a monotonically decreasing convex function. In other words, this means that the increase in diffusion length that we obtain by shortening the lattice constant $l$ by an amount $\delta l$ exceeds the decrease that we would obtain by expanding the lattice constant by the same amount to $l+\delta l$. Thus, when assuming a symmetric distribution of lattice constants around $l$, the averaging would lead to a slightly larger diffusion length. However, we note that this effect might be too small to be observed in our case, due to the presence of other sources of static disorder.