Efficient light harvesting and photon sensing via engineered cooperative effects

Efficient devices for light harvesting and photon sensing are fundamental building blocks of basic energy science and many essential technologies. Recent efforts have turned to biomimicry to design the next generation of light-capturing devices, partially fueled by an appreciation of the fantastic efficiency of the initial stages of natural photosynthetic systems at capturing photons. In such systems extended excitonic states are thought to play a fundamental functional role, inducing cooperative coherent effects, such as superabsorption of light and supertransfer of photoexcitations. Inspired by this observation, we design an artificial light-harvesting and photodetection device that maximally harnesses cooperative effects to enhance efficiency. The design relies on separating absorption and transfer processes (energetically and spatially) in order to overcome the fundamental obstacle to exploiting cooperative effects to enhance light capture: the enhanced emission processes that accompany superabsorption. This engineered separation of processes greatly improves the efficiency and the scalability of the system.

Since the discovery of coherent features in natural light-harvesting complexes [1][2][3][4][5][6][7][8] and subsequent studies of the functional role of these features, there has been great interest in engineering biomimetic devices for photon sensing or light harvesting, able to exploit coherent quantum features [9][10][11] even in ambient conditions. Natural lightharvesting complexes are composed of organic chromophores, each characterized by a dipole moment that determines its coupling to the electromagnetic field (EMF) and its interaction with neighboring systems. Once light is absorbed, the induced photoexcitation is transmitted to another molecular aggregate, called the reaction center (RC), where charge separation occurs, which drives subsequent steps in the photosynthesis chain. The main quantum coherent effects that are thought to contribute to the high efficiency of natural photosynthetic complexes are induced by the delocalization of the excitation over many molecules [12][13][14][15][16][17][18][19][20][21]. Such delocalized excitonic states can have an enhanced dipole strength that strongly couples them to the EMF. Thus, these states are able to absorb light at a rate much larger than the single-molecule absorption rate. Indeed, the absorption rate of a single delocalized excitonic state can increase linearly with the number of molecules over which the excitation is delocalized. On the other hand, the states that absorb light efficiently also exhibit enhanced emission rates (termed superradiance) since the process is time reversible. The time-reversal character of absorption/emission processes is one of the main reasons for the Shockley-Queisser efficiency limit for photocell devices [22].
Several recent proposals have aimed to suppress re-emission in systems composed of few molecules. Specifically, it has been discussed how coherent effects can suppress re-emission leaving absorption intact for a molecular dimer [23][24][25], and a mechanism for suppressing re-emission, exploiting bright and dark states and fast thermal relaxation, has been devised for two [9,26] and three coupled molecules [27]. Moreover, in Ref. [11] it is shown how to maximize cooperative absorption and engineer super-absorbing manyatom structures that avoid superradiance by controlling the structure of a collection of atoms/molecules and engineering their vibrational environment to achieve delicately tuned thermal transition rates. The above results suggest that one can in principle exploit cooperative effects to enhance light capture and transfer by designing engineered structures that avoid detrimental effects such as superradiance.
In this work we propose a molecular architecture that is able to suppress re-emission, while leaving absorption intact. The general idea of our device is based on engineering a super-absorbing state at high energy. Once the excitation is absorbed, it is transferred to the low-energy states by thermal relaxation, which, being much faster then re-emission, prevents radiative losses. The low-energy states transfer the excitation to the central core absorber with reduced radiative losses due to the fact that their dipole strength is smaller than the high-energy absorbing state.
Specifically, our light-harvesting device is composed of a ring of N molecules surrounding a central core absorber, similar to photosynthetic complexes found in purple bacteria. The molecular arrangement is engineered so that we have three bright excitonic states with orthogonal dipole moments, one at high energy and two at low energy.
Only the low-energy states are coupled to the central core absorber. By changing the orientation of the molecular dipoles, we can control the brightness of these states. If we make the high-energy state the brightest, absorption mainly occurs through it. That state is not directly coupled to the central core absorber, but fast thermal relaxation funnels the excitation towards the low-energy states. The brightness of these cannot be zero if we want to exploit the radiative coupling with the central core absorber to transfer and trap the excitation. Nevertheless, to minimize re-emission, it is convenient to keep the low-energy brightness rather small. Balancing these two requirements, one can find optimal parameter ranges in which transfer is maximized.
The efficiency of our proposed architecture is analyzed as a function of the ring radius, which determines the number of molecules needed to keep a given density, and the orientation of the molecules in the ring. Related work providing theoretical insight into the role of fold symmetry in promoting efficient energy transfer in LH2 can be found in [28]. Since our proposed architecture suppresses radiative recombination, the advantages of our design are present only if radiative recombination is the main cause of efficiency losses. For instance, for very large trapping rates at the RC, radiative losses become negligible. Note that the mechanism proposed here extends the design of Ref. [9] for two molecules to a device composed of an arbitrary number of molecules. Moreover, while the transfer mechanism in Ref. [9] was not radiative, here we consider that the coupling to the central core absorber is radiative in nature, as it happens in natural photosynthetic systems.
Our model device can be tuned to mimic natural light-harvesting systems, and we compare the performance of the optimized device to that of a model of a natural system, under weak laser excitation and a realistic model for natural sunlight. The efficiency of our device, in the optimal size range N ≈ 50, is found to be more than two orders of magnitude larger than that of a single absorber and enhanced by a factor larger than N . In the same regime, as we will show below, the efficiency of models mimicking natural systems is enhanced only by a factor equal to the number of absorbers. We also show that, under natural sunlight excitation, the decoupling mechanism of our model leads to an improvement in the efficiency of about five times with respect to a non-decoupled configuration.

The structural model
We consider a molecular complex composed of N molecules (or point absorbers) placed on a ring of radius R with constant density d, sketched in Fig. 1a. The excitation energy of each molecule, ω 0 , is assumed to be the same while their dipole moments µ n have the same modulus, µ, but possibly different orientations, indicated by the unit vector p n . The properties of the system will be studied for different system sizes, keeping the density d fixed, while varying the radius. The role of the ring structure is to absorb the electromagnetic radiation and transfer the excitation to a central core absorber, such as the reaction center of natural photosynthetic complexes, that we mimic with an additional central site, coupled to an external environment where the excitation can be irreversibly trapped. The whole system is also coupled to a thermal bath at fixed temperature T . The molecules lie on a ring in the xy plane with their dipole moments tangential to the circumference and tilted by an angle ±θ with respect to the xy plane, see Fig. 1a. Moreover, the vertical components are alternated upwards and downwards, so that the normalized dipole moment orientations arê p n = cos θφ + (−1) n sin θẑ , whereφ andẑ are the unit vectors corresponding, respectively, to the azimuthal and to the vertical direction of a cylindrical coordinate system. Due to the discrete rotational invariance of the system around the z axis, the eigenstates and their dipoles, in the single excitation manifold (that is appropriate to describe the weak-field limit) take the form |E α = n c n (E α ) |n , with c n (E α ) = 1 √ N exp(i2παn/N ), and the corresponding dipoles are p α = n c n (E α )p n , with a square modulus ranging from zero to N . Here, |n is the state where the n-th molecule is excited and all the other ones are in their ground state. For any θ, the dipole strength | p α | 2 is non-vanishing only for three excitonic states: |E 2 , |E 3 and |E N , where we have ordered the excitonic states by increasing energy. Actually, |E 2 and |E 3 span a degenerate subspace and, without loss of generality, we choose them to be two orthogonal states in this subspace such that their dipole moments are withx andŷ being the unit vectors of the planar axes. The third bright eigenstate is the highest-energy exciton, whose dipole moment is perpendicular to the ring plane. Since the emission rate from a state is proportional to | p α | 2 , the states |E 2,3,N are also superradiant, with an emission rate proportional to N . All the other states are subradiant with zero emission rate (dark states). Moreover, the enhanced dipole of the high-energy state |E N is orthogonal to the xy plane, and to the dipoles of the other bright states: this guarantees the separation of excitation and transfer processes. We add a RC to this model by an additional site |rc placed at the center of the ring, with excitation energy ω rc , and coupled to an external environment (sink) where the excitation can be trapped at rate κ. The RC is dipole coupled to the choromophores on the ring, and we choose the dipole moment of the RC along the y axis,p rc =ŷ. We also set the energy of the RC site to be resonant with the first excited state of the ring. As a consequence, only the state |E 2 , with dipole moment along the y direction, has a non-vanishing coupling strength to the reaction center (see Appendix C), Since the density of dipoles on the ring d = N/(2πR) is kept constant, then Ω C scales with N as This coupling determines the transfer between the ring and the RC. Since the sum of the dipole strengths of all the eigenstates must be constant, α | p α | 2 = N , increasing θ from 0 to π/2, the dipole strength of the high-energy state (|E N ) is also increased, so that the smallest dipole strength of the low-energy states is decreased [see Eqs. (2)]. This configuration limits radiation losses together with ring-RC transfer. Nevertheless, as we discuss in the manuscript, the transfer is much faster than the radiative losses below a critical ring size, thus preserving the trapping efficiency.
In order to describe this ring+RC system interacting with an electromagnetic field and with a thermal phonon reservoir, we use a master equation, see section 2. The reservoir coupling is assumed to be Markovian, and each molecule is assumed to couple to an independent Ohmic bath at the same temperature. The parameters of the bath, see Appendix E, have been chosen so that thermal relaxation among exciton states occurs in about a picosecond at room temperature for N = 32. This is comparable with estimates for natural photosynthetic systems reported in literature [29,30], and it is much faster than the emission timescales, of the order of nanoseconds.
As a measure of efficiency of our device we use the stationary current transmitted from the central site to the sink, while the system is driven by the EMF, defined as where ρ rc (t) is the population of the RC at time t. The current is further divided by the maximal stationary current, I s , of the RC alone (in absence of the ring) under the same illumination conditions. In this way the normalized current I/I s measures the increased efficiency of our network of dipoles with respect to a single site. Since the excitation can only be trapped in the central site, a normalized current larger than unity indicates increased effectiveness of the network of sites in absorbing and transferring the excitation.

Hamiltonian and Master Equation
In our model we choose the values of parameters to be close to realistic values found in natural complexes such as LHI-LH2 in purple bacteria [29,30]: squared dipole moment µ 2 = 519310Å 3 cm −1 (corresponding to µ ≈ 10 D); excitation energy of the ring sites ω 0 = 12911 cm −1 (corresponding to a single-site transition wavelength λ 0 ≈ 775 nm); The dynamics of our model is described by the following master equation, written in a rotating frame with respect to the driving field mode frequency ω [31][32][33][34]: Here, the Hamiltonian H S = H 0 + ∆ + H EM captures the evolution of the system in the weak-field limit, where no more than one excitation is induced in the system. Specifically, H 0 = N n=1 (ω 0 − ω) |n n| + (ω rc − ω) |rc rc| represents the site energies of the ring chromophores and RC, and ∆ represents the Coulomb coupling between chromophores. Here we assume that the ring chromophores are distant from each other and from the RC, so that each molecule can be approximated as a point dipole. Explicitly, the matrix elements of the coupling ∆ are Herer nm := r nm /r nm is the unit vector joining the n-th and the m-th sites, and p n := µ n /µ is the normalized dipole moment of the n-th site. The dielectric constant is r = 1, which is a good approximation for molecules surrounded by air. In principle, the nearest-neighbor coupling in the ring should be computed without using the point dipole-point dipole approximation, because the distance between the chromophores is comparable to the molecular size. Nevertheless, this is a detail which does not qualitatively change our results. For instance, the nearest-neighbor dipolar couplings used in this manuscript range between ≈ 500 − 1200 cm −1 , which is comparable to the ≈ 400 − 800 cm −1 couplings estimated from detailed electronic calculations in Ref. [29].
On the other hand, the coupling with the central core absorber can be safely assumed to be a point dipole-point dipole coupling (as it has been done also in Refs. [29,30]) since the molecules in the ring are far apart from the central core absorber. The expressions in Eq. (7) are valid in the small volume limit, where the wavelength of the optical transition is larger than the system size (λ 0 R), which is the regime where natural light-harvesting complexes operate. The full expressions, without this approximation can be found in Appendix A.
The term H EM in the Hamiltonian describes the coupling between molecules and the continuous-wave (CW) driving laser and it is given by where Ω R = µE 0 / is the Rabi frequency, E 0 is the amplitude of the electric field,ˆ is a unit vector which specifies the laser polarization, k is the wave vector of the laser field, |q represents the system sites (either the ring sites |n or the reaction center |rc ) and r q is the position of the q-th site. |0 is the ground state of all molecules in the system. In our calculations we always choose | k| ≈ 2π/λ 0 and since we are in the small volume limit, we can approximate the matrix elements of Eq. (8) as We have confirmed the validity of this approximation for our ring system, as long as R/λ 0 < 0.1, see Appendix B. Now we return to the other terms in the master equation, Eq. (6): L f l and L rc are Lindblad dissipators derived under the Born-Markov and secular approximations [35] and they describe, respectively, fluorescence emission of the molecules and transfer to the RC, while R T is a non-secular Redfield dissipator [35] modelling thermal relaxation and decoherence in the presence of a thermal bath. The dissipators read explicitly where the sums over m, n run over all the system sites (ring sites or RC), a n = |0 n| (here, n| can be a ring site or the RC), a rc = |0 rc|, and Γ mn ≈ γp n ·p m in the small volume limit (R λ 0 ), with γ = 4 3 µ 2 k 3 0 / r . Here, k 0 := ω 0 n r /c, where c is the speed of light and n r the refractive index. For the realistic parameters chosen here, the decay width of a single molecule is γ = 3.7 × 10 −4 cm −1 . We also set n r = 1, which is a good approximation when the system is surrounded by air ‡. Again, for a discussion about the regime beyond the small volume limit see Appendix B. Finally, R T describes dissipation due to the coupling of each molecule to an Ohmic bath, where are the thermal rates, depending on the spectral density J(ω) and on the Bose distribution n BE (ω) of the phonons which form the bath and More details about R T can be found in Appendix E. Note that, for the coupling to the thermal bath, we use the non-secular Redfield dissipator R T instead of the commonly used Lindblad dissipator since we found that, in our model, the secular approximation is not valid and produces unphysical results, a well-known issue in molecular excitonic transfer [36]. Specifically, when the coupling Ω C is very small, see Eq. (4), the secular approximation incorrectly predicts that the transfer rate between the ring and the RC becomes independent of Ω C , while the Redfield dissipator R T correctly predicts that the transfer rate tends to zero with Ω C . Although the Redfield master equation is known to produce negative populations in the intermediate-to-strong system-bath coupling (see Ref. [9] and references therein), we checked that all the steady-state populations are positive within the parameter range that we analyzed. On the other hand, the secular approximation is valid for modeling fluorescence decay and decay to the reaction center, and thus we can keep the L f l and L rc dissipators in their Lindblad form, see Eqs. (10) and (11). ‡ This value of γ corresponds to a fluorescence time τ fl = 14 ns and differs from the excitation lifetime ∼ 1 ns found in literature [12], because here γ represents only the radiative decay processes and nonradiative decay is neglected. In the case of pure water, one should set n r = 1.33 and r = n 2 r = 1.77, thus obtaining γ = 4.9 × 10 −4 cm −1 and τ fl = 11 ns. For a proteic environment, instead, it is usually set r = 2.3 [17] which, keeping the refractive index of water, gives the same γ obtained in air. coupling Ω C = (µ 2 /R 3 ) N/2 cos θ with |E 2 ; the ring eigenstate |E 3 , whose dipole strength is p E3 = N/2 cos θx, and is decoupled from the RC; and the highest-energy ring eigenstate, |E N , whose dipole strength is p E N = √ N sin θẑ. The radiative decay rates of the states depend on the dipole strength of the excitonic states and are given by γ α = γ| p α | 2 . The absorption rates for the D-configuration (a) and LH-configuration (b) depend on the laser frequency, its intensity and polarization, see Eq. (19). (c): Effective three-level model. Schematic representation of the effective three-level model which is able to capture the main properties of the many-level system, see Eqs. (17a,17b,17c).
Our model has been derived under the single-excitation approximation. This is a good approximation of a realistic situation only if the excited state population is much smaller than unity, which is true for our choice of the parameters (see Appendix G for more details).

Illumination conditions
We consider three types of electromagnetic field states illuminating the device. Firstly, under what we call the D-configuration (Fig. 1b, more details in Fig. 2a), we consider a coherent, CW monochromatic polarized field as it was done in Ref. [23]. The polarization axis is chosen to be aligned withẑ, which means that it couples to the highest-energy excitonic state in the device, see Eqs. (2). Second, under what we call the LH-configuration (Fig. 1c, more details in Fig. 2b), we consider a coherent, CW monochromatic field polarized in theŷ direction and incoming perpendicular to the ring. Such a field only excites the low-energy ring eigenstate |E 2 . As one can see from Fig. 2b, absorbing and transfer states coincide in this set-up, and such model is a good representative of some natural light-harvesting complexes (see Appendix D).
Finally, under what we call the Sunlight configuration (Fig. 1d), we model illumination by natural sunlight, which is isotropic, unpolarized, incoherent and broad-band. This is modeled well as black-body radiation at 6000 K [37,38]. Specifically, in the sunlight configuration the Hamiltonian term H EM is not present, while we include two additional Lindblad dissipators for absorption and stimulated emission induced by sunlight, where n S ≈ 0.04 is the Bose occupation of the Sun photons at the excitation energy ω 0 and at the Sun temperature (6000 K) and f S = 5.4 × 10 −6 accounts for the Sunto-Earth distance [39]. Specifically, under sunlight illumination each eigenstate acquires absorption and stimulated emission rates, B α = f S n S γ| p Eα | 2 , with f S representing the solid angle of the Sun as seen on Earth, with r S being the radius of the Sun and R ES the Sun-to-Earth distance. Finally, the rates are proportional to the squared magnitude of the eigenstate dipole strength, so that only the states |E 2,3,N have a non-vanishing sunlight absorption rate. The intensity of natural sunlight is 1365 W/m 2 , and we consider the same intensity also in the Dconfiguration and the LH-configuration. In those configurations the light intensity is encoded in the Rabi frequency, Ω R = µE 0 / . By imposing the intensity of the CW laser to be E 2 0 /(4π) = 1365 W/m 2 (using Gaussian units), we determine the corresponding value of E 0 and, from that, the Rabi frequency, which is Ω R = 4.68γ (in units of the single-molecule radiative decay rate, γ ≈ 0.07 ns −1 ). We keep this value of Ω R fixed in all the manuscript.

Effective three-level model
Here we show that the dynamics of the complex structure described above, under all the illumination conditions considered, can be mapped to the dynamics of an effective threelevel incoherent model with the relevant quantum effects encoded in few parameters. The three-level model is described by the zero-excitation state |0 , a single-excitation state |e for the whole ring, and a single-excitation state |rc for the RC, see Fig. 2c. The excitation pumped by the EMF into the ring is quickly funneled to the low energy states by thermal relaxation. Therefore, we determine the rates between |0 , |e and |rc assuming that the ring is always at thermal equilibrium with the phononic reservoir. Under this assumption (see Appendix I), the emission rate from |e to |0 is the thermal average of the ring eigenstate emission rates, γ = α e −Eα/(k B T ) Z γ α , while the transfer rate from the ring to the RC is also a thermal average, T RC = α e −Eα/(k B T ) Z T RC α , involving the transfer rates T RC α between each α eigenstate and the RC, that are proportional to the squared coupling between |α and |rc . Note that the T RC rate is equivalent to the well-known multi-chromophoric Förster resonance energy transfer (MC-FRET) [21,40] or generalized Förster theory [41] rate, as we show in detail in Appendix I. In our specific case, due to the ring symmetry, only the |E 2 eigenstate has a nonvanishing transfer rate, T RC 2 = τ −1 RC (32/N ) 5 cos 2 θ, where τ RC = 3.9 ps is the transfer time between the ring and the RC at θ = 0 and N = 32 (more details in Appendix I), while T RC α = 0 for all α = 2. On the other hand, the absorption rate is the sum of all the absorption rates, B T OT = α B α (where the absorption rates, B α ∝ | p α | 2 , have different expressions whether the excitation is induced by a CW laser or by sunlight, see Appendix I), and the transfer rate from the RC to the ring is also the sum of all the transfer rates, T RC T OT = α T RC α . Finally, the stimulated emission rate is again a thermal average, B = α e −Eα/(k B T ) Z B α . The RC also can absorb the incoming radiation with an absorption rate B RC , that accounts also for stimulated emission, while its emission rate is γ. This approach yields the following rate equations for the populations of the three levels, Solving for the steady state of these equations (details in Appendix I), we obtain an approximation to the steady-state transmitted current, The validity and effectiveness of this three-level model is discussed in the next section, see also Ref. [39]. Note that the effective three-level model presented in Eqs. (17a,17b,17c) is able to describe the whole system, composed by the ring and the central core absorber, under both the pumping from a light source (laser or sunlight) and thermal relaxation. Our effective three-level model is based on the assumption of fast thermal relaxation, incoherent pumping and incoherent transfer between the ring and the core absorber. Specifically the coupling between the ring, assumed at thermal equilibrium, and the central core absorber is described with an approach equivalent to the the generalized Förster theory, see Appendix I. Finally note that in literature three-level models describing exciton transport have been widely used [20,38,[42][43][44]. In particular our approach is similar to the one used in Ref. [38] where the pumping of sunlight on a dimer system has been considered.

Super-absorption in the low-fluence regime
First we demonstrate that the molecular device developed above is capable of exploiting cooperative effects to enhance the absorption from a weak-intensity EMF. For any θ > 0 and N , only three ring eigenstates have a non zero dipole strength: two in the low-energy region (|E 2 and |E 3 ) and one with the highest energy (|E N ). Concerning the lowenergy states, |E 2 has a polarization along y while |E 3 along x. In contrast, the high-energy state has a polarization along the z-axis. Thus for an EMF polarized in the z direction (D-configuration) the absorbing and the transferring states are separated: only the highest-energy state |E N is coupled to the EMF, while only the low-energy state |E 2 can transfer the excitation to the RC.
At high (e.g., room) temperature the pumping rate for this system under CW laser excitation can be described semi-classically by the Förster rates [43] where (ω − ω abs ) is the detuning frequency of the laser with respect to the absorption frequency and p abs is the dipole strength of the absorbing state. Γ T is the dephasing rate of the coherences between the absorbing state and the ground state. We compute Γ T analytically in Appendix F, and show that it depends only on the density of states of the system and on the parameters of the bath. Critically, Γ T is independent of N and very weakly dependent on θ. Thus T L ∝ | p abs | 2 ∝ N . This demonstrates what we call superabsorption: the absorption is concentrated in a very specific system eigenstate characterized by a giant dipole, and the absorption rate grows proportionally to the system size. Note that our definition of superabsorption refers to the low-fluence regime, which is the focus of this manuscript. Under high fluence, cooperative absorption is instead characterized by a super-linear absorption rate, as it has been shown in Ref. [11].

Scalability and efficiency
Now we demonstrate that superabsorption can work in concert with the engineered supertransfer from ring to RC, to result in a photocurrent that scales with the system size.
Here we analyze the efficiency of our device under a laser field polarized along the z direction and under the action of a thermal bath at room temperature. In particular we will analyze the dependence of the current on the laser frequency and of the peak currentĪ (maximal current obtained at the optimal laser frequency) as a function of the system size. Next, we compare the efficiency of the D-configuration with an alternative illumination condition, the LH-configuration, where the absorbing and the transferring states coincide, thus mimicking natural light-harvesting complexes more closely. Finally, we evaluate our device in the Sunlight configuration, which is a realistic model of illumination by natural sunlight. At the same time, we show that the results of the full N -level system can be captured by the simpler three-level system introduced in the previous section. Fig. 3 shows the dependence of the photocurrent on the CW laser frequency both for the D-configuration (Fig. 3a, where the field is assumed polarized in the z direction) and for the LH-configuration (Fig. 3b, with the field along y), for different system sizes at room temperature. For the D-configuration we choose an angle θ = 0.475π that gives the optimal current for N = 64 and a close-to-optimal current for N = 16, 32 (see dashed line in Fig. 3c). On the other hand, for the LH-configuration we choose the optimal angle θ = 0, where the current is maximal (see Fig. 3d).
The combined effect of superabsorption (at high-energy), thermal relaxation and transfer (at low energy) results in a peak in the transmission spectrum at the high energy of the absorbing state (that is higher than ω 0 , see vertical dashed lines in Fig. 3a) and not to that of the transferring states (that would be lower than ω 0 ). Since the high-energy state is totally decoupled from the RC, the peak at its frequency can only be explained by thermal relaxation after absorption. Note also that the height of the peak increases with the system size due to cooperative absorption. In Fig. 3a we also show as continuous curves the results of our analytical three-level model Eq. (18), which reproduces the current across the entire frequency range.
In Fig. 3b we show the ratio between the current I and the single-site current I s [see Eq. (G.5)] for the LH-configuration as a function of the laser frequency at room temperature, for different system sizes. The intensity has a peak (vertical dashed lines) when the laser frequency is resonant with the low-energy bright eigenstates of the ring, that are resonant with the RC. It is interesting to note that even for the LH-configuration the height of the peak increases with the system size due to cooperative absorption. However, the peak current obtained with the LH-configuration is about three times smaller than with the D-configuration (compare Fig. 3a and Fig. 3b). Also for the LHconfiguration, the three-level model [Eq. (18)] reproduces very well the results of the full system (see continuous lines in Fig. 3b).
Therefore, in order to understand whether the separation of the absorbing and transmitting states that we engineered can improve the scalability of the system, we compute the peak transmitted current as a function of θ and of the system size.
In Fig. 3(c,d,e) we plot the peak currentĪ vs. θ for different values of the system size. In each panel, the results of the master equation are compared to the effective three-level model, Eq. (18). Various illumination conditions are considered. For the D-configuration (Fig. 3c), the peak current increases with θ up to a maximum, close to π/2. This is a consequence both of the absorption rate T L increasing with θ for the D-configuration, and of the emission rate being suppressed on increasing θ. On the other hand, for the LH-configuration (Fig. 3d) the peak current is maximal for θ = 0 and it decreases with θ. This is a consequence of the absorption rate decreasing with θ in the LH-configuration and of the fact that absorption and emission are not decoupled. Then, for the Sunlight configuration (Fig. 3e) the current is enhanced by increasing θ, as a result of the suppression of emission. In all the three cases shown, the current is enhanced on increasing N from 16 to 64. Finally, as one can see from the figure, the three-level model (lines) gives a good approximation of the master equation results (symbols) in all the ranges considered. Small deviations for large N and for θ close to π/2 can be explained by a variation of the value of τ RC used in the three-level model, which in the figure has been kept fixed as N and θ vary. Indeed we set τ RC = 3.9 ps which has been obtained by fitting the master equation for N = 32 and θ = 0. Nevertheless τ RC can vary by up to 20% as it is shown in Appendix I. If we account for those variations in τ RC in our calculations of the current, we obtain a perfect agreement also for large N , see shaded areas in panels (c,d,e). Deviations for small N , see panel (d), are due to the fact that the couplings between the ring and the RC are large and the energy transfer is not fully incoherent as it is discussed in Appendix I.
Since we are interested in the scalability of the device at large N , in the following we use the three-level system, that is much less computationally expensive than the master equation at large N . In Fig. 4 the normalized maximal currentĪ 3 /I s is shown vs. θ and N at room temperature, as obtained from the effective three level model, see Eq. (18) for the trapping rate κ = 10γ and for different pumping mechanisms (i.e. D-configuration, LH-configuration and Sunlight configuration). In all cases, we can see that the current increases at first with N . Moreover, for the D-configuration the efficiency improves with increasing θ, it reaches an optimal value for 40 < N < 80 and it ultimately decreases with N for very large ring sizes. Such improvement with θ can also be seen for natural sunlight pumping and has been observed and commented above in Fig. 3(c,e).
The scaling of the current with the system size can be understood as follows. For small size, N , the excitation is cooperatively absorbed by the ring and efficiently transferred to the RC where it is trapped. Indeed, a small ring radius implies a strong dipole coupling to the RC and, therefore, a fast transfer. So, the trapped current for small N ultimately scales as the absorption rate, increasing with N . On the other hand, for large sizes N , a large ring radius implies a weak coupling to the RC, that decreases as ∼ |Ω C | 2 ∼ 1/N 5 , see Eq. (4). Such suppression of the transfer to the RC acts as a bottleneck for large N , so that the trapped current decreases with N for large N in all cases. Moreover for large N the thermal population in the superradiant state coupled to the central absorber decreases as 1/N , thus quenching the current for large system sizes.
Moreover, in Fig. 4 the normalized currentĪ 3 /I s of the LH-configuration (ypolarized) is shown as a function of N and θ and it is compared with that of the D-configuration for the trapping rate κ = 10γ, that is of the same order of the emission rates ( γ ≈ γ).
Finally, we also analyze the model in the Sunlight configuration. In this case, the pumping is incoherent, broad-band and isotropic. As a consequence, the total absorption rate of natural sunlight, B T OT = N γf S n S , is proportional to N and independent of θ. Nevertheless, even in this case for κ = 10γ we see an increment of the current on increasing θ, because the system benefits from the suppression of emission.
About the choice of the trapping rate we note that this is critical to the efficiency of our set-up. Indeed a very large trapping rate prevents re-emission, since the excitation is quickly trapped, and thus makes less useful the suppression of re-emission which we consider here. The trapping rate can vary a lot depending on the specific system. In photosynthetic antenna complexes a charge separated state is created very quickly once the excitation is absorbed in the reaction center (few picoseconds, corresponding to κ ≈ 10 4 γ). On the other hand, charge transfer in the RC is much slower and the RC is not active until the charge-separated state is neutralized again, and this occurs on the order of 100 µs, corresponding to κ ≈ 10 −4 γ. In the figures presented here an intermediate trapping rate has been considered, κ = 10γ, corresponding to a trapping time of ≈ 1 ns. Even if this trapping rate is only slightly faster than the emission rate (≈ 10 ns), a considerable advantage is obtained in the D-configuration. As one can see from Fig. 4, for the trapping rate κ = 10γ the maximal current obtained from the D-configuration is about three times higher than the LH-configuration. Note that our proposed design is useful only if the emission plays an important role. If, instead, the trapping rate κ is so fast to overcome any emission process, there is no advantage in the decoupling mechanism. Specifically, the efficiency of the D-configuration, as opposed to the LH-configuration and to a single site, improves on decreasing the trapping rate κ, as we show in the following.
In Fig. 5 we study how the efficiency of our proposed device depends on the trapping rate κ. We analyze a broad range of reasonable values of the trapping rate: the lower bound, κ = 10 −4 γ ≈ (100 µs) −1 , corresponds to the reset time in purple bacteria reaction centers [12], while the upper bound, κ = 10 4 γ ≈ (1 ps) −1 , is the charge separation rate in purple bacteria [12].
Specifically, in Fig. 5 we plot the ratioĪ 3 /(N I s ) between the maximal steady-state current obtained from the three-level modelĪ 3 , Eq. (18), divided by N times the singlesite current, I s , for N = 32 as a function of the trapping rate κ/γ and of the angle θ. Different panels represent different illumination conditions. For the D-configuration, see Fig. 5a, for any κ, the current is enhanced by increasing θ, because the absorption is increased and at the same time the emission is suppressed. It is interesting to note that, for any fixed θ, the current increases with κ. Indeed, for these parameters, the emission rate in our device is of order ≈ N γ cos 2 θ ≈ γ (for N = 32, as in Fig. 5), and therefore a fast trapping rate κ γ allows to overcome the emission rate. Similar comments can be applied to the sunlight configuration, see Fig. 5c. Also in such case, for fixed θ, the current increases with κ, because the fast trapping rate overcomes the emission. In the sunlight configuration, for κ = 10 4 γ the current becomes basically independent of θ and equal to N times the single-site current, see also Fig. J1c and Appendix J for more details.
For the LH-configuration, see Fig. 5b, a similar pattern can be seen: the current increases with κ for θ fixed. This configuration does not have a decoupling mechanism between absorption and transfer, and therefore the LH-configuration shows the best performance under the trivial conditions: (a) optimal absorption rate at θ = 0 and (b) maximal trapping rate (in the figure, κ = 10 4 γ).
We also point out that, both for the D-and the LH-configuration, there is a broad range in the (κ, θ) parameter space where the current is more than N times larger than the single-site current I s . Interestingly, the normalized currentĪ 3 /(N I s ) in Fig. 5a-b can be larger than unity for the LH-configuration only for large trapping rates, while for the D-configuration the normalized current can be larger than unity even for small trapping rates.

Robustness to disorder
The efficiency of our proposed device will be affected by disorder, that can be due for instance to fluctuations in the positions of the sites, in the orientations of their transition dipoles, or in the site energies. Here we study how the efficiency of our device, measured by the peak steady-state currentĪ/I s , is affected by disorder in the dipolar orientation. We introduce the angular disorder as follows: each dipole is given a random orientation inside a cone centered on the precise orientation of Eq. (1). All the cones have the same solid angle, that can take values from 0 to 4π. The magnitude of the solid angle represents the disorder strength: a vanishing solid angle represents no disorder, while the maximal solid angle 4π represents completely disordered dipoles. The energy of the RC is always equal to the energy of the first excited state of the ring without disorder and, similarly, the frequency of the CW laser is determined by the energies of the ring at zero disorder: the CW laser is resonant to the highest-energy ring eigenstate for the D-configuration, and it is resonant to the first-excited state of the ring for the LH-configuration.
In Fig. 6a we plot the average normalized peak steady-state current, Ī /I s , against the strength of angular disorder for N = 32, different illumination conditions and the corresponding optimal angles θ: for the D-and sunlight configurations we consider the optimal angle θ = 0.475π discussed above (see Fig. 3c), while for the LH-configuration we show the case θ = 0 (which is optimal, see Fig. 3d). As one can see from Fig. 6a, in all cases the current is suppressed by disorder. Such suppression is very sharp for the D-and LH-configurations, and milder for the sunlight configuration. Specifically, the efficiency of our system is robust to very high disorder for the sunlight configuration, but it is also robust for the D-configuration. Indeed, as one can see from Fig. 6a, the current for the D-configuration remains very high for angular disorder as large as 10% of the full solid angle. On the other hand, the LH-configuration appears less robust to disorder, and the current quickly decreases for angular disorder larger than 1% of the full solid angle (see Fig. 6a).
Moreover, in Fig. 6b we consider the D-configuration for θ = π/2 compared to the optimal case (θ = 0.475π). In the case θ = π/2, at zero disorder, the only ring eigenstate with non-vanishing dipole strength is the highest-energy ring eigenstate, with polarization perpendicular to the ring plane, see Eq. (2). Therefore, at zero disorder the ring is decoupled from the RC and there is no current. However, one can see that a small amount of angular disorder is able to increase the current up to very high values, which are comparable to the optimal configuration θ = 0.475π in absence of disorder, see Fig. 6b. The reason is that a small angular disorder slightly increases the dipole strength of the lowest-energy states of the ring: in this way, the coupling between the ring and the RC is activated giving rise to a current. Actually, the emission from the ring remains low, because most of the dipole strength remains concentrated in the highestenergy state (for small disorder). These results, therefore, suggest an alternative way to engineer our device: instead of a finely-tuned configuration with θ very close to the optimal value 0.475π, one can use an existing ring where all the dipoles are perpendicular to the plane (θ = π/2) and a bit of angular disorder is present. In such configuration, the emission is suppressed while it is possible to transfer the excitation to the RC via the low-energy ring states. So far, we kept the energy of the RC fixed and at resonance with the first and second excited states of the ring, in the absence of disorder. However, it is known that the spectral width is affected by disorder, therefore one may ask whether putting the RC at resonance with the zero-disorder energy levels is the best choice in the presence of disorder. Therefore, in Fig. 7 we show how the peak current changes as a function of the angular disorder and of the energy of the RC. To guide the eye, in each panel we plot the average energies of the three lowest eigenstates of the ring as a function of disorder: the lowest excitonic state |E 1 (black continuous line), and the first and second excited states |E 2 and |E 3 (black dashed lines). Analyzing Fig. 7, one can see that, for the LH-and sunlight configurations, the optimal efficiency is always obtained with the RC energy red-shifted with respect to the ring states by some thousands of cm −1 . Indeed if the RC is at resonance there is a large backward transfer to the ring states and excitation can be lost by re-emission. If the detuning is too large, then transfer to the RC is suppressed and again re-emission lowers the efficiency. The optimal position of the energy of the RC follows the ring ground state energy (black continuous lines in each panel) as it is modified by disorder. As a final note we stress that in absence of disorder the LH-configuration with optimal detuning can even be more efficient than our proposed device configuration, compare panels a) and b) in Fig. 7. This shows that several paths are available to suppress re-emission, and red-shifting the central core absorber energy can also be very effective in suppressing re-emission. This mechanism has also been discussed and exploited in Ref. [39] by some of the authors of this manuscript. As a final remark, let us note that exploiting thermal relaxation to suppress re-emission requires detailed knowledge of the system-bath coupling and bath structure, while the mechanism proposed by us in this manuscript is more direct and does not depend on the details of the system-bath coupling to be effective.

Conclusions and perspectives
In many natural light-harvesting complexes, most of the dipole strength is concentrated in few states that absorb light and, at the same time, transmit the excitation to an external trapping environment. The large dipole strength of such states favour the absorption of light but it also induces losses by re-emission of the excitation, thus limiting the efficiency of the energy transfer.
Here we propose a light-harvesting device in which the absorbing and the transferring states are engineered to be different by structural arrangement of chromophores. We proved that our engineered device is able to improve the efficiency of light-harvesting complexes by several orders of magnitude both when the interaction with a polarized monochromatic field is considered and under natural sunlight. Since the solar spectrum is broad, to use the proposed device as solar light-harvesting complex, an ensemble of devices absorbing at different frequencies should be considered. The proposed device can also be extended to unpolarized light by arranging the rings on a spherical shell, in a similar way to the arrangement found in the chromophores of purple bacteria [12].
Our approach allowed us also to study the scalability of our device as the number of light-harvesting chromophores is increased by increasing the ring radius and keeping a fixed chromophore density. We have shown the existence of an optimal size. The reason for that is that even if increasing the ring radius improves absorption, it also suppresses transfer towards the central core absorber. As a future development we plan to consider different architectures to make the system scalable to larger sizes without losing efficiency. In particular, instead of increasing the ring radius, a network of smaller rings where the excitation is efficiently transferred between them and finally concentrated in the central core absorber should be able to improve the efficiency of our device at larger sizes. Indeed this is precisely the architecture of several photosynthetic natural systems [12,29,30]. A similar idea has been successfully employed by some of the authors of this paper in Ref. [39] where a bio-inspired sunlight-pumped laser has been proposed.
In order to realize the proposed complex molecular structure a precise control over molecular orientation and structure is required. In this context, modern molecular synthesis and modification techniques can easily meet these needs. Several nanostructures, including nanotubes [45], DNA proteins [46], and viruses [47,48] can be precisely functionalized with organic molecules using an impressive variety of bioconjugation tools. Engineering of synthetic molecular aggregates in linear, circular, and other geometric configurations, with nanometer separation between molecules, is commonplace. The major challenge facing these techniques is controlling energetic and structural disorder after functionalization. However, recent experiments demonstrating the use of functionalized DNA proteins for light harvesting [46] constitute a proof-ofprinciple confirmation of the promise of such synthetic molecular engineering techniques for light capture technologies. We also note that the effect of structural disorder could also be exploited at our advantage. As shown above, in a perfect H-aggregate, i.e. a molecular aggregate characterized by a single bright state above the energy of the monomer absorption peak (Fig. 6b, θ = π/2), all the dipole strength is concentrated in the highest excitonic state, and structural disorder is able to add some dipole strength to other states allowing the lowest excitonic state to also transfer energy to the central core absorber. In this way we achieve the same separation of absorption from transfer shown to be so effective in improving light-harvesting efficiency. Finally, a lesson can be drawn about natural photosynthetic antenna complexes from our analysis: natural LH2 systems can be very efficient at absorbing and transporting light excitations if emission is efficiently suppressed either by a large trapping rate, see Fig. 5c, or by properly detuning the RC energy, see Fig. 7d. In this sense it is likely that natural systems also exploit suppression of superradiant emission while using super-absorption to enhance their efficiency.

Acknowledgments
This research was supported, in part, by the Center for Research Computing of the University of Notre Dame through access to key computational resources. We acknowledge financial support from Fondazione E.U.L.O. in the frame of the project "Trasporto quantistico in sistemi nanostrutturati con applicazioni ai biosistemi". Sandia National Laboratories is a multimission laboratory managed and operated by National Technology & Engineering Solutions of Sandia, LLC, a wholly owned subsidiary of Honeywell International Inc., for the U.S. Department of Energy's National Nuclear Security Administration under contract DE-NA0003525. This paper describes objective technical results and analysis. Any subjective views or opinions that might be expressed in the paper do not necessarily represent the views of the U.S. Department of Energy or the United States Government.

Appendix A. Effective Hamiltonian
The effective Hamiltonian of an aggregate interacting with an electromagnetic field can be written as [31][32][33]49] where the latter expression refers explicitly to the site basis |n . In Eq. (A.1) ω 0 is the energy of each site. In the limit where only one excitation is present in the system [31][32][33]49] the diagonal and off-diagonal matrix elements are given respectively by, where γ = 4 3 µ 2 k 3 0 / r , µ = | µ| is the transition dipole matrix element, r the relative dielectric permittivity,p n := µ n /µ the normalized dipole moment of the n-th site, x nm = k 0 r nm ,r nm := r nm /r nm the unit vector joining the n-th and the m-th sites, and k 0 = ω 0 /c. A derivation of the expressions in Eqs. (A.2c) and (A.2d) is presented in Appendix H.
In the small volume limit x nm 1, that is when the system size is much smaller than the wavelength λ 0 = 2π/k 0 connected with the optical transition, the matrix elements can be approximated as Appendix B. Emission rates and coupling to the laser beyond the small volume limit Eigenvalues and eigenstates can be obtained by diagonalizing the effective Hamiltonian, Eq. (A.1). From them it is possible to define the emission rate γ α associated with each eigenstate with eigenvalue E α − i γ α /2. In the small volume limit one can approximate the emission rate as For any value of the angle θ between the dipoles and the ring plane, the dipole strength of the ring eigenstates is non-vanishing for just three bright eigenstates. Two of them correspond to the degenerate subspace of the first and second excited states. Without loss of generality, we choose two combinations of those two states so that their dipole moments are withx andŷ being the unit vectors of the planar axes (note that E 1 is the excitonic ground state energy). The third bright eigenstate is the highest-energy one, whose dipole moment is perpendicular to the ring plane. For the above three states, the dipole strength increases with the system size due to cooperative effects induced by the symmetric arrangement of the dipoles in the ring. In order to study the range of validity of the small volume approximation, in Fig. B1 we compare the normalized maximal decay width Γ max /γ obtained from the complex eigenvalues of H eff (open circles) with that obtained in the small volume limit using Eq. (B.1) (continuous line) for different values of R/λ 0 where R is the radius of the system ring. Note that we increase R keeping a fixed density, meaning that the number of sites N increases proportionally to R. As one can see significant deviations appear already for R/λ 0 > 0.1.
In the previous paragraph we analyzed the emission rate of the system eigenstates. Now we focus on the absorption rate. In particular, let us assume that the ring system is coupled to a laser, described as a monochromatic electromagnetic wave. The interaction Here Ω R = µE 0 / is the Rabi frequency (E 0 being the amplitude of the electric field),ˆ the laser polarization and k the wave vector of the incident field. H EM describes the photon absorption and stimulated emission, and it induces coherent oscillations between the ground state |0 (i.e. the state without any excitation) and the single excitation states |n . The absorption rate of each eigenstate |E α is proportional to the square of the following matrix element: where we used the decomposition of the ring eigenstates in the site basis. When the size of the system is smaller than the wavelength of the laser ( k · r n 1), the coupling can be expressed using the dipole strength of an eigenstate Since | k 0 | ≈ | k|, for small ring sizes (| k 0 |R 1) we are in small volume limit approximation, therefore both Eq. (B.7) and Eq. (B.1) are satisfied. This means that for each eigenstate, the absorption and the emission rate are both proportional to the squared dipole strength. In Fig. B1 we plot the squared dipole strength of the eigenstate with highest energy (as a continuous line) vs. the ratio R/λ 0 . In the same figure we also plot the maximal normalized emission Γ max /γ (open black circles) and the maximal normalized absorption | 0| H EM |E α | 2 max / Ω R 2 2 (full red dots). We can observe that, for R/λ 0 < 0.1 both the maximal emission Γ max /γ and the maximal absorption are well approximated by the squared dipole strength. Since the maximal dipole strength is proportional to N (see Eqs. (B.2)) we have a cooperative coupling to the laser. On the other hand, when k · r n ≈ k 0 · r n 1, the approximations (B.1) and (B.7) are not valid and both the maximal absorption and the maximal emission rate grow slower than linearly with N .

Appendix C. Analytical coupling to the RC
Here we discuss the coupling between the ring eigenstates and the reaction center (RC) in the small volume limit. The reaction center is modeled by a site with excitation energy ω rc and the same transition dipole moment µ of the ring sites. The dipoledipole coupling between two sites is given by Eq. (A.3a) and the ring eigenfunctions by Eq. (B.5). The coupling between a ring eigenstate and the RC is In our model, the radial component of the dipoles is vanishing, so that p n ·r n = 0 ∀n = 1, . . . , N .
This leads to a simplification of the expression (C.1), which becomes We have the reaction center dipole oriented along the y axis, i.e.
so that only one of the N eigenstates (which belongs to the doubly degenerate subspace of the first and second excited states and which we call |E 2 ) is coupled to the reaction center, with a coupling  Figure D1. Top panels: schematic representation of the directions of the transition dipoles for the purple bacteria LHI complex taken from Ref. [50] (left) and for the LH-configuration (right). Lower panels: normalized decay rates Γ α of the eigenstates vs. their energy E α and vs. the index α = 1, . . . , N . The Hamiltonian for the LHI complex is taken from Ref. [29]. Here the system size is N = 32 for both the LHI complex and the LH-configuration.
In this manuscript, the density d = N/(2πR) is kept constant, so that Ω C scales with N as Appendix D. The LH-configuration as a representative of natural light-harvesting complexes The molecular structure introduced in the main text in the LH-configuration is a good representative of some natural light-harvesting complexes. In Fig. D1 on the left panels we show the LHI complex of purple bacteria [29,50], while on right panel we show the LH-configuration with θ = 0. On the top panels we show a schematic representation of the directions of the transition dipoles in the models, while in the lower graphs we show the normalized decay rates of the eigenstates of the systems vs. their energy. The positions and dipole orientations of the LHI complex have been taken from Ref. [50] and the Hamiltonian parameters from Ref. [29]. In both configurations the decay rates are concentrated in the first and second excited states, which are degenerate and have the same decay rate ≈ N γ/2.
The interaction with the laser field is described by Eq. (B.3). Note that we neglected the counter-rotating terms, according to the rotating wave approximation (RWA, [31,32]). Moreover, by a unitary transformation, the time dependence of the laser term has been removed, leading to a shift of the diagonal terms in the Hamiltonian (i.e. the site energies), which become In order to describe the dynamics of our model coupled to a thermal bath, we consider the following Hamiltonian: where ∆ is defined in Eq. (A.3a) and we use the following master equation [31][32][33][34], where the last three terms describe, respectively, the fluorescence, the trapping in the RC and the thermal dissipation in presence of a thermal bath. They are given by where Γ mn are given by Eq. (A.2d), a n = |0 n|, a rc = |0 rc| and R T [ρ] is the thermal dissipator. In this scheme, each site is assumed to be coupled to an independent Ohmic bath with linear coupling. Nevertheless they all follow the same dynamics, having the same temperature, spectral density and coupling strength [51,52]. Specifically, the Hamiltonian of the system coupled to the independent thermal baths reads where g k is the linear coupling strength of a site with the harmonic oscillator with frequency ω k . Here, X k,n = /(2M ω k ) b † k + b k is the position of the harmonic oscillator (M here is the mass of the oscillator). The Redfield dissipator, Eq. (E.6), is derived under the Born-Markov approximations [35], assuming weak system-bath coupling and fast bath relaxation time, but without applying the secular approximation, as discussed in the main text. Taking the continuum limit for the sum over the bath frequencies ω k and defining the spectral density J(ω) through we obtain Eq. (E.6) [34,53], with In Eq. (E.9) n BE is the standard Bose-Einstein distribution of the phonons n BE (ω) = 1 e ω/k B T − 1 (E. 10) and the spectral density is chosen as [34] The linear dependence of the spectral density at small frequencies, J(ω) ∼ ω for ω ω c , results from the ∼ 1/ω scaling of the squared system-bath coupling (see Eq. (E.8), lefthand side), multiplied by the density of modes of the oscillators in 3D, ∼ ω 2 dω. We use an exponential cut-off at large frequencies, J(ω) ∼ e −ω/ωc for ω ω c , because it has been used to reproduce spectroscopic results in similar molecular aggregates [34,53]. As regards the bath parameters, we set the reorganization energy to E R = 200 cm −1 and the cut-off frequency to ω c = 333 cm −1 . With this choice of the parameters the thermal relaxation among exciton states occurs in about 1 ps at room temperature for N = 32, which is comparable with the estimates for natural photosynthetic systems reported in literature [29,30], and which is much faster than the times obtained by the radiative emission rates γ α ∼ 1 ns −1 . The operators in Eq. (E.6) can be expressed as Here |Ẽ α are the eigenstates of (H 0 + ∆), according to andc n (Ẽ α ) = n|Ẽ α . Note that here the system includes both the ring and the RC, thus the expression of the coefficientsc n (Ẽ α ) is different from the one given in Eq. (B.5), where just the ring is considered.

Appendix F. Thermal dephasing rate
Using the expressions in Appendix E, we have that the dynamics of the coherences between the ground state |0 and the ring eigenstates in the small volume limit follow with the overlap coefficients Eq. (F.1) can be simplified in the present case, under polarized CW laser excitation. Coherences between the ground state and the ring eigenstates are created by the laser hamiltonian term H EM . In our analysis, such coherent coupling involves just one ring eigenstate: either |E N for the D-configuration, or |E 2 for the LH-configuration. No coherent coupling with the ground state is present in the Sunlight configuration. Therefore, in Eq. (F.1) the only non-zero term is where α = β, coinciding with the absorbing state, that we label "abs" in the following. Thus, after re-labelling the summed index δ → β, Eq. (F.1) simplifies to with the simplified overlap coefficients In a more compact expression, we can write where we define the dephasing rate Using the spectral density, Eq. (E.11), we obtain: where E R is the reorganization energy, ω c is the cut-off frequency, T is the bath temperature, and the coefficients Λ abs,β are given in Eq. (F.4).
In the case of the device described in the main text (D-configuration), when the absorbing state is the highest-energy state, we can make some approximations. The coefficients (F.4) have the value Λ abs,β = 1/N for any β so, using Eq. (B.2) for p abs and integrating over the spectrum we have where (E) is the density of states. Since in our case the spectral extension E N − E 1 is independent of N , then the average density of states is Moreover, in realistic situations (T = 300 K and the parameters E R and ω c from the literature) the first term in Eq. (F.8) can be neglected for N 10 6 . Thus, under these approximations, we can claim that Γ T is independent of N .
For the LH-configuration we need to use the full expression given in Eq. (F.6) which can be obtained by the parameters of the system and by the diagonalization of the system Hamiltonian. Also in this case we have found (numerically) that Γ T is independent of N , as shown in the following. In Fig. F1 we plot Γ T vs. N for the LH-configuration (θ = 0, open circles) and the D-configuration (squares for θ = π/3 and crosses for θ = 5π/12). As one can see, for N 20 the thermal dephasing rate is weakly dependent of N in all cases. Moreover, the two values of θ considered for the D-configuration exhibit the same dephasing rate for any N 20.
We also study the dependence of Γ T on the angle θ for the D-configuration, see Fig. F2, for N = 100. We can see that Γ T has a very weak dependence on θ (less than 3% variation), too.
So, our analytical approximation (supported by numerical simulations) show that the dephasing rate depends very weakly on both N and θ. For instance at T = 300 K and for the parameters considered in the main text, we have that Γ T ≈ 6 ps −1 for LH-configuration (F.10a) Appendix G. Analytical results for a single site Here we analyze the explicit expression of the Master Equation (6) in the case of a single site which plays the role of an absorber, emitter and trapping state. We call respectively |0 and |1 the ground state and the excited state of the site, so that Eq. (6) reads where ∆ 0 = ω 0 − ω is the detuning between the laser frequency and the transition frequency of the site. As one can see from Eq. (G.1c), the coherence term has a dephasing rate which can be determined by neglecting the terms proportional to iΩ R and i∆ 0 , which are related to oscillations. Note that this dephasing rate is exactly the same given by Eq. (F.6), in the trivial case β = 1. The dephasing rate of the coherence term ρ 01 can also be interpreted as the dephasing rate between the ground state and the (unique) absorbing state, which we call Γ T . In the following we will show that it is possible to obtain a stationary solution for Eq. (G.1). Using (E.9), Eq. (G.2) can be rewritten as 3) The stationary current I s of a single molecule at fixed temperature T is defined as: where ρ ∞ 11 is the steady-state value of ρ 11 (t). The explicit expression for I s can be derived analytically by setting the derivatives in Eq. (G.1) to zero, .

(G.5)
It is well known that, for sufficiently large dephasing (Γ T Ω R /2), the full quantum master equation is well approximated by a set of rate equations with suitably defined rates [43]. Specifically, following Ref. [43], we derive the pumping rate T L between the states 0, 1 as With such rate we proceed to write a probability-preserving rate equation for the probabilities P 0 and P 1 to be in ground and in excited state respectively, dP 0 dt = T L P 1 − T L P 0 + (γ + κ)P 1 (G.7a) where we have taken into account the pumping rate T L (G.6). The incoherent transmitted current is the stationary value of the probability P 1 (t), can be easily obtained from (G.7) and it is given by Interestingly, the result in Eq. (G.9) coincides with the exact quantum result given in Eq. (G.5) for any value of the parameters. Note that in general (when the number of sites is larger than one) a set of effective rate equations like Eqs. (G.7) gives different results from a quantum master equation like Eq. (G.1). Now, having computed the stationary current, we analyze the stationary population of the excited state |1 . Our interest is to determine when the single-excitation assumption (under which the master equation (6) has been derived) is valid, as a function of the parameters. Let us start by considering the resonance case, ∆ 0 = 0 so that the stationary solution of Eq. (G.1) is . (G.10) Let us recall that the effective Hamiltonian (A.1) and the master equation (6) have meaning only for low excitation (ρ 11 1). In our simulations we have Ω R = 4.68γ, (γ + κ) > γ and 2Γ T > γ (p) (0) ≈ 2 × 10 6 γ. This implies that ρ ∞ 11 10 −5 , so that we are always in the single-excitation regime.

Appendix H. Natural sunlight
In this section we derive the master equation for a generic molecular aggregate coupled to natural sunlight.
Let us consider an aggregate of N two-level systems all having the same excitation energy ω 0 . In these calculations we use = 1. The aggregate interacts with the radiation emitted by the Sun, that we model as a black body at temperature T S , with a correction accounting for the Sun-to-Earth distance. In this approach, the spontaneous emission process will come out naturally from the interaction with the vacuum mode of the electromagnetic field. The full Hamiltonian iŝ Here the site Hamiltonian iŝ withσ z j being the z Pauli matrix for the j-th site. The black body Hamiltonian iŝ where the summation runs over the modes k and the polarizations λ = 1, 2 of the field, the dispersion relation is ω k = ck (c is the speed of light) and the creation/annihilation operators follow to the commutation rules [b k,λ ,b † k ,λ ] = δ k, k δ λ,λ . Finally, the lightmatter interaction Hamiltonian iŝ is the dipole operator on the j-th site, d j is the transition dipole moment of the same site,σ ± j = (σ x j ± iσ y j )/2 and is the electric field in the position r j , with e k,λ being a unit vector which specifies the polarization.
Since the coupling to the EMF degrees of freedom is weak, we perform the Born-Markov and secular approximations [35]. The density matrix is therefore factorized aŝ ρ(t) ≈ρ S (t) ⊗ρ B , and we get the Lindblad master equation whereÂ j (ω) are operators acting on the sites, The complex rates G ij (ω) are , (H.9) where we use the notation . . . B = tr B {. . .ρ B }. Now we assume that the black body, that in our case is the Sun, is at thermal equilibrium, i.e.
where T S = 6000 K is the temperature of the Sun. In this case the expectation values of the operators in Eq. (H.9) are b k,λb k ,λ B where we have defined the Bose-Einstein occupation of the Sun photons and we have introduced the factor that represents the ratio of the Sun solid angle as seen from the Earth over the full solid angle (r S is the Sun radius while R ES is the Sun-to-Earth distance). Such factor is needed to correct the model to describe natural sunlight: in the present calculations, in fact, the molecular aggregates exchange energy with all the modes of the EMF in all directions, both absorbing and emitting photons. This is correct for the spontaneous emission process, but it needs to be corrected by the factor f S for the absorption and stimulated emission processes, which depend on the number of thermal photons n S (ω k ). Indeed, the Sun photons hit the system from a fraction f S of the whole solid angle, so that n S is multiplied by f S . Thus, defining r ij = r i − r j , Eq. (H.9) can be written as . (H.14) As regards the sum over k, we take the continuum limit Now, if we assume that the dipoles have all the same magnitude µ but different orientation, namely d j = µp j , and defining the function choosing a frame for the integration over k where the z axis has the same direction as r ij , we have Now we perform the integral over τ using the relation where P is the Cauchy principal value. So, we can split the rates into their real and an imaginary parts, which are, respectively, Appendix H.1. Real Part of the rates: absorption and decay Let us start from the real part (H.20). The two integrals are easily performed, taking into account that the only possible values of ω are ±ω 0 . By defining k 0 = ω 0 /c we get To have the explicit dependence of Γ ij (ω) on the parameters, we evaluate F ij (x), that results Note that F ij (x) is an even function of x which, in our case, gives the useful equality F ij (−k 0 r ij ) = F ij (k 0 r ij ). Moreover, one can see that F ji (x) = F ij (x), which implies that both the matrices Γ ij (ω) and S ij (ω) are symmetric for i ↔ j. As regards the diagonal terms (i = j) we can analytically extend the function to x = 0 thanks to the limit The real parts of the rates are then where the single-molecule spontaneous decay rates are Using the symmetry properties of Γ ij , the contribution from the real part of the rates to the master equation reads which is in the Lindblad form and where we have defined the coefficients (H.28) The first double sum of Eq. (H.27) describes the spontaneous and stimulated emission processes, while the second double sum describes the absorption process of excitation from the Sun.
Appendix H.2. Imaginary Part of the rates: radiative coupling Let us now focus on the imaginary part of the master equation. Thanks to the symmetry S ji (ω) = S ij (ω) we have where we have defined the real part of the radiative Hamiltonian using (H.8) Thanks to the commutation rules σ where we have defined the matrix elements The diagonal terms give a divergent renormalization energy that is constant for all the molecules, so we disregard it here. Then, the off-diagonal matrix elements are independent of the Sun temperature, namely where x 0 = ω 0 r ij /c. The integral can be computed using contour methods, resulting in

Appendix H.3. Final Expression and Single-excitation approximation
The final expression of the master equation is and, defining the parameters x ij = ω 0 r ij /c and γ = 4 As regards the Hamiltonian term, neglecting the renormalization of the site energies, we haveĤ where the coupling terms are given by (H.34). The master equation (H.35) acts on the full Hilbert space spanned by the N sites, having dimension 2 N and including all the possible numbers of excitations (from none to N excitations). However, in this manuscript we focus on the weak fluence regime, where the photon absorption rate is much smaller than the excitation decay rate, so that no more than one excitation at a time is present in the system. In this regime we perform the single-excitation approximation, i.e. we neglect all the states with more than one excitation and we consider only: the state |0 , where all the sites are in their ground state, and the N single-excitation states of the form |j =σ + j |0 , where only the jth site is excited while all the other ones are in their ground state. In this (N +1)-dimensional subspace, each σ + j operator acts only on |0 resulting in σ + j |0 = |j , while each σ − j operator acts only on |j giving σ − j |j = |0 . Therefore, we write the master equation in the single-excitation approximation by replacingσ + j → |j 0| and σ − j → |0 j| into (H. 35). For readability, we also drop the subscript "S" fromρ S , and we have the single-excitation master equation with the single-excitation Hamiltonian Finally, let us consider the particular case where there exist a common eigenbasis |α for both Ĥ 0 +∆ and i,j Γ ij |i j| such that We can then write (H.38) in that basis: If we consider the diagonal elements, which describe the dynamics of the populations of |0 and of the eigenstates |α , that part of the master equation can be mapped into a Pauli master equation, which reads where we have defined the absorption and stimulated emission rates B α = f S n S (ω 0 )γ α .

Appendix I. Effective three-level model
Here we provide a detailed derivation of the three-level model introduced in the main text. Let us consider a system made of N sites (molecules) and let us perform our analysis in the single-excitation approximation, meaning that we consider just the ground state of the whole system |0 and the states |j where only the j-th site is excited. Let us now consider a basis where the single-excitation subspace is diagonal, and let us call |α the eigenstates of the single-excitation subspace.
The system is assumed to be coupled to an incoming radiation field, which induces absorption and stimulated emission to each eigenstate |α with a pumping rate B α . Moreover, we consider the spontaneous emission of excitation by radiation from each eigenstate |α , with a fluorescence rate γ α . Finally, we add one level (labelled "RC") which is coupled to each |α state with a transfer rate T RC α and where excitation can be collected to a trapping environment (sink), modeled by a trapping rate κ. The RC can also absorb excitation from the radiation field with a rate B RC , and it also has stimulated and spontaneous emission rates, B RC and γ.
Neglect coherences As a first approximation, we assume that the coherences between the eigenstates do not play a role in the transport process. We then write a rate equation for the population of the ground state and for the RC state: where P 0 is the population of the ground state, P rc is the population of the RC and P α is the population of the α-th excitonic eigenstate.
Thermal equilibrium Secondly, let us assume that the excitonic subspace is at thermal equilibrium with a temperature T . Formally, we define the population in the whole excited subspace as so that the trace preservation condition is P 0 (t) + P e (t) + P rc (t) = 1 .

(I.3)
We impose thermal equilibrium in the aggregate as where E α is the energy of the α-th excitonic eigenstate and is the partition function. By substituting Eq. (I.4) into Eq. (I.1) we have Steady-state solution and current We now express Eq. (I.6) in terms of thermal averages of the rates ( X = α X α p α ) and by defining the total absorption rate, B T OT = α B α , and the total transfer rate to the RC, where the last equation is Eq. (I.3). The steady-state solution is obtained by setting the time derivatives to zero, and it is reached at long time ("t = ∞"), so that we have Finally, we define the stationary current trapped into the sink as Here we show that the transfer rates T RC and T RC T OT between ring and RC are exactly the multi-chromophoric Förster resonance energy transfer (MC-FRET) rates [21,40], also known as generalized Förster theory [41]. In the MC-FRET framework, the transfer rate from a donor aggregate "D" to an acceptor "A" is expressed in terms of their emission and absorption spectra, E(ω) and I(ω) respectively. An aggregate absorption spectrum is [21] I(ω) ∝ α | p α | 2 I α (ω), where | p α | 2 and I α (ω) are, respectively, the dipole strength and the normalized lineshape for each α eigenstate. The emission spectrum on the other hand is E(ω) ∝ α | p α | 2 E α (ω), where the emission lineshapes are multiplied by the thermal populations p α , see Eqs. (I.4) and (I.5), namely E α (ω) = p α I α (ω). The MC-FRET rate is usually expressed as [21,40,41] where α|H S |β is the Hamiltonian matrix element between the α donor eigenstate and the β acceptor eigenstate and the normalization condition is ∞ −∞ I α (ω)dω = 2π. For high temperature and short bath correlation time [21], we can neglect the phonon-induced Stokes and anti-Stokes shifts and approximate all the absorption lines as Lorentzians peaked on the eigenstate frequency ω α and with a dephasing-induced linewidth Γ φ . Under this assumption, the overlap integral in Eq. (I.12) is analytically computed as (I.14) Therefore, we can express the MC-FRET rate in Eq. (I.12) as where the transfer rates between an α donor eigenstate and a β acceptor eigenstate arē (I. 16) Note that these rates are symmetric,K β,α =K α,β , while the K D,A rate usually are non-symmetric, due to the different thermal populations p α between the donor and acceptor. The transfer rates between ring and RC from our rate equations, Eqs. (I.7), are exactly equivalent to the MC-FRET rates in Eq. (I. 15), in fact we have two cases: (i) for the transfer from the ring to the RC, the ring is the donor and the RC the acceptor, so the sum over β ∈ A runs over the single RC state and we have K ring,RC = α∈ring p αKα,rc , (I. 17) which is exactly T RC ; (ii) for the transfer from the RC to the ring we have the opposite situation, so the sum over the the donor states α ∈ D runs over the single RC state, whose normalized population is trivially p α = 1, and so we have which is exactly T RC T OT .

Appendix I.2. Transfer between ring and RC
Since the coupling between the ring and the RC is weak, we compute the incoherent transfer rates between the ring eigenstates and the RC using the MC-FRET rates in Eq. (I. 16), T RC α ∝ | α|∆|rc | 2 , proportional to the squared coupling between the two states. As we show in Appendix C, the coupling is non-vanishing only for the eigenstate |E 2 , that is resonant with the RC, where the coupling scales as Ω C ∝ (cos θ)/N 5/2 , see Eq. (4). Therefore, the transfer rate in Eq. (I. 16) can be written in the form where τ RC represents the transfer time in the reference case N = 32 and θ = 0 (representing the natural purple bacteria LHI complex, see Appendix D In order to determine the value of τ RC , now we consider the ring+RC system in absence of absorption, emission and trapping, where the excitation can only be exchanged between the ring and the RC. In such situation, Eq. (I.7) simplifies to dP e (t) dt = − T RC P e (t) + T RC T OT P rc (t) (I.23a) where the two rates are related by T RC = T RC T OT p * (here p * is the Boltzmann population of the eigenstate |E 2 ). If one excitation is initially present on the ring, the time evolution of the probability to be on the RC is (I. 24) and it allows to determine the value of τ RC , as explained below.   Table I1. The initial quadratic growths P rc (t) ≈ p * Ω 2 C t 2 / 2 [with Ω C given by Eq. (4) In Fig. I1 we compute the time evolution of the ring+RC system using the master equation (see Appendix E) in absence of absorption, emission and trapping, and in presence of a thermal bath with the standard parameters used in this work (see caption). We initialize the system with one excitation at thermal equilibrium on the ring, i.e.
and we compute the time evolution of the population of the RC, P rc (t) = rc|ρ(t)|rc (symbols in Fig. I1). We fit the results obtained this way with the three-level solution Eq. (I.24), leaving τ RC and p * as fitting parameters. As one can see from the figure, P rc (t) obtained from the master equation (symbols) grows initially quadratically in time (dashed lines) as This quadratic growth is the pure quantum-mechanical time evolution given by the ring |E 2 eigenstate (with initial occupation probability p * ) resonant with |rc . For times larger than ≈ 0.01 − 0.1 ps, the time evolution is instead well captured by Eq. (I.24) (continuous lines), with the fitting parameters reported in Table I1. We observe that, in the reference case N = 32, θ = 0 (representing the natural LHI system), there is a perfect fit with τ RC = 3.9 ps, and therefore we choose this value for all our simulations with the three-level model across the manuscript. The corresponding dephasing value for the MC-FRET rates is obtained from Eq. (I. 22), and it has the value Γ φ = 34 ps −1 .
We also note that τ RC varies by about 20% when θ is changed (see Table I1). Such ±20% variation in τ RC produces variations in the steady-state current, and we show those variations as shaded areas in Figs. 3, J1 and J2. Finally, we determine the critical time τ φ when the dynamics changes from quadratic quantum-mechanical growth to the linear growth predicted by the incoherent three-level model. This transition happens at the crossing point between Eq. (I.26) and the initial, linear growth of Eq. (I.24), that for short times is approximated as P rc (t) ≈ p * T T OT RC t. Matching the linear and quadratic expressions for P rc (t) and using Eq. (I.20) for T T OT RC , we have that Ω 2 C cancels, so that the crossing point happens at τ φ = Γ −1 φ . With the choice of τ RC and the corresponding Γ φ used in this manuscript, we have τ φ = 0.03 ps, independent of N or θ, as one can see in Fig. I1.

Appendix I.3. Parametrization for D-configuration, LH-configuration and Sunlight configuration
Here we express the relevant rates present in Eqs. (I.7) as a function of the system parameters, in absence of disorder. Due to the symmetry, simple analytical expressions can be derived, as it is shown below. Note however that, in presence of disorder, the following expressions are not valid and one must use instead the general definitions of the rates [see Eqs. (I.6), the paragraph below those equations and Eqs. (I.7)].
Let us start with the D-and LH-configurations, where absorption takes place from a polarized CW laser source. where p N and p 2 = p 3 =: p * are, respectively, the Boltzmann populations for the eigenstate with the highest energy, and for the first and second excited eigenstates, which are degenerate and, thus, they have the same Boltzmann population p * . Here, p rc =ŷ is the unit vector indicating the dipole direction of the RC, while Γ T,RC is the dephasing rate between |0 and |rc , and it is given by Eq. (G.3). Sincep rc is orthogonal to the polarizationẑ of the laser, we have B RC = 0. The phenomenological parameter τ RC represents the transfer time between the ring and the RC for N = 32 and θ = 0, as explained above. In the range of parameters that we have studied in our simulations, at T ≈ 300 K we are justified to neglect the factor p N . In the main text we have shown that the three-level model reproduces the peak current very well for the trapping rate κ = 10γ. Here we show that the three-level model works very well also for weaker and stronger trapping rates. Specifically, in Fig. J1 we plot the peak current as a function of θ for the D-, LHand Sunlight-configurations for large trapping rate, κ = 10 4 γ. As one can see, the results obtained with the master equation (symbols) are very close to the three-level model results (lines). Some discrepancies are present in the LH-configuration (Fig. J1b) for N = 16 and θ ≈ 0 because the three-level model in that regime does not capture the ring-RC transfer process. In fact for those parameters, as one can see from Fig. I1, the ring-RC population transfer follows purely quantum-mechanical quadratic dynamics, and therefore the three-level model is not expected to work. Deviations are also seen for large N in the D-and Sunlight configurations. As discussed in the main text and in Appendix I, those deviations are due to our choice of a constant, θ-independent τ RC parameter in our simulations, while Table I1 shows that τ RC can vary by up to 20% with θ. Once we account for ±20% variations in τ RC , we obtain a confidence interval for the three-level current (see shaded areas in figure) that include the master equations results. For the Sunlight configuration (Fig. J1c), for strong trapping, κ = 10 4 γ, the current is independent of θ, because in this limit the current it is ultimately determined just by the absorption rate, which is N f S n S γ, independent of θ.
On the other hand, in Fig. J2 we show the peak current as a function of θ for the D-, LH-and Sunlight-configurations for small trapping rate, κ = 10 −4 γ. Also here, the results obtained with the master equation (symbols) are nearly identical to the threelevel model results (lines), thus justifying our use of the three-level model in the main text.
Moreover, in Fig. J3 we show the peak current obtained from the three-level model as a function of θ and N for large trapping, κ = 10 4 γ. A similar figure is present in the main text, see Fig. 4, showing the effectiveness of the D-configuration for intermediate trapping values, κ = 10γ. Here, see Fig. J3, we show that for large trapping, κ = 10 4 γ, the D-configuration is actually performing worse than the LH-configuration.