Temporal correlations of sunlight may assist photoprotection in bacterial Photosynthesis

Photosynthetic systems utilize adaptability to respond efficiently to fluctuations in their light environment. As a result, large photosynthetic yields can be achieved in conditions of low light intensity, while photoprotection mechanisms are activated in conditions of elevated light intensity. In sharp contrast with these observations, current theoretical models predict bacterial cell death for physiologically high light intensities. To resolve this discrepancy, we consider a unified framework to describe three stages of photosynthesis in natural conditions, namely light absorption, exciton transfer and charge separation dynamics, to investigate the relationship between the statistical features of thermal light and the Quinol production in bacterial photosynthesis. This approach allows us to identify a mechanism of photoprotection that relies on charge recombination facilitated by the photon bunching statistics characteristic of thermal sunlight. Our results suggest that the flexible design underpinning natural photosynthesis may therefore rely on exploiting the temporal correlations of thermal light, manifested in photo-bunching patterns, which are preserved for excitations reaching the reaction center.

Natural photosynthesis thrives under constantly changing environmental conditions thanks to active regulatory mechanisms. These mechanisms must balance the interplay between thermal light absorption, transport of the resulting electronic excitations (excitons) to a reaction center (RC), charge separation and metabolic output. However, the sheer complexity of these four individual processes makes their combined description challenging and fundamental questions remain open. In particular, photosynthetic bacteria can survive in light intensities that exceed predictions from current theoretical models [1][2][3][4][5][6][7]. To resolve this discrepancy, we propose a novel approach to describe the coupled dynamics resulting from those fundamental processes in real photosynthetic vesicles, as depicted in Fig.1.
Illumination. The coherence length of sunlight at the earth's surface is approximately 50 µm [8,9]. This value exceeds the size of any unicellular photosynthetic organism, e.g. a purple bacterium, as illustrated in Fig.1(a). Hence, the optically active molecules of photosynthesis, which arrange in light harvesting complexes (LH1 and LH2 in purple bacteria) within photosynthetic vesicles, will absorb thermal photons with significant temporal correlations leading to interspersed bursts of photons, as noted decades ago [10][11][12].
Charge separation dynamics. Once an exciton arrives at a RC it triggers an electron transport chain, whose outcome is the production of quinol molecules (Q B H 2 ). This transport goes through intermediate metastable states of charge separation with relevant recombination times (t crit in Fig.1(c)) [15][16][17], which in combination with the bursted structure of the incoming excitons may effectively reduce the conversion efficiency of the photosynthetic machinery. These processes are described in detail in the next section.
Metabolism. Quinol reduction facilitates the formation of an electrochemical gradient that drives the rotation of a transmembrane macromolecule ATP-synthase (ATPase) to synthesize Adenosine Triphosphate (ATP) (Fig.1(d)). The single AT-Pase present per vesicle may synthesise at most 100 molecules (about 200 excitations/s) of ATP for bacterial photosynthesis [2]. Based on this analysis, typical membranes with 400-600 LHs [1,25] and 90% transport efficiency [26,27], will accumulate hydrogen, increasing dangerously the cytoplasmic acidity for intensities as low as I =0.1-0.2 Watt/m 2 . This estimate is in stark contrast with the proven survival of purple bacteria for light intensities as high as 100 Watt/m 2 [1,27]. In principle, charge recombination could provide the means to avoid the formation of excess quinol, but it has not been considered in previous attempts to provide a global description of the primary steps in bacterial photosynthesis [26][27][28][29][30].
In this article, we show that a detailed but comprehensive model of the elementary processes in the photosynthesis of purple bacteria unveils a dynamical interplay between the bunched statistics of thermal light, the process of charge recombination and quinol production. In particular, the long time-intervals between photon bursts resulting from thermal light absorption provide an essential element to determine the photoprotective role of metastable states in RCs.

A. Thermal light illumination and excitons statistics
The photon-statistics of thermal light absorption shows typical intensity fluctuations that lead to patterns of bunched photo-detections separated by waiting times (where no detection happens), when the detection takes place within the On the left side, we illustrate the parameters describing the incoherent light source provided by sunlight. Here a denotes the source diameter, D the distance to the light reception system and λ the mean wavelength of the incident light. The large orange circle with diameter l c ≈ 50µm, [9] represents the coherence area of sunlight, encircling a complete photosynthetic bacterium of length ≈ 40µm [31]. It is within this coherence region defined by l c where thermal temporal correlations will become manifest. On the right side of panel (a) we depict the photosynthetic apparatus of the bacterium, embedded in characteristic vesicles. They contain the light harvesting complexes LH2 and LH1 and the encircled reaction centres (RCs), with typical inter-complex distance r k,l l c . (b)-(d) Flow diagram illustrating the elementary processes in bacterial photosynthesis. Photons (depicted as orange bolts) are detected by the light harvesting complexes (LH1 and LH2) and converted into photo-excitations totally or partially delocalized along individual harvesting complexes (excitons), which travel through the network of LHs until reaching an open RC as depicted in (b). Once in the RC, the exciton induces charge separation and a chain of electron transport steps that end up in quinol (Q B H 2 ) production. Charge can be recombined, as represented by the red arrow, with recombination lifetime t crit in (c). In (d) quinol oxidizes and the liberated hydrogen (together with two additional hydrogens pumped from the cytoplasm) produce the electrochemical gradient required to make the ATPase macromolecule rotate and produce ATP. The ATPase molecule catalyzes the reaction between adenosine diphosphate (ADP) and inorganic phosphate (P i ). coherence length and coherence time of the light [32][33][34]. This pattern -also called burstiness of the detected photontraces -arises as a consequence of the light-matter interactions and hence depends on the joint properties of the light field (frequency spectrum, intensity and correlations) and the detection system provided by the LH1-2 complexes (absorption efficiency, absorption bandwidth, spatial configuration and detection time window). We have implemented a statistical model of thermal light detection where the probability distribution of absorption events is characterized by its factorial moments generating function [35]. This framework, whose main features are revised in the Methods, allows for the in-clusion of light correlations in a natural manner, showing that the degree of burstiness can be tuned by the ratio between the detection time window and the coherence time of the detected light T/τ c . This parameter allows the continuous statistical tuning of the absorbed excitations from a thermal (T/τ c 1) field with maximal burstiness, to a Poissonian field (T/τ c 1) without burstiness, as shown in Fig.2(a).
The statistics of bursts in a sequence of absorption events can be characterized by the so-called burstiness parameter, B = (σ t − t )/(σ t + t ). This quantity vanishes in a fully Poissonian process where the mean inter-event time t equals the standard deviation of waiting times σ t , but has a positive  In the inset, we show the characterisation of low-B signals using as figure of merit the quantities t , σ t , t intra , and t inter (descriptions in the text). All quantities are displayed in milliseconds on a semi-logarithmic scale. We considered10 6 light absorption events.
value whenever the times series exhibit bursts of events [36]. Unsurprisingly, in Fig.2(b) we show that the ratio T/τ c and B are related monotonically. This implies that either of them can be used to quantify the bunching of the light absorption, and supports our use of B in what follows as a measure of the burstiness in traces of events, in particular, on the arrival of excitations to the RCs. To characterise the statistical properties of the waiting times t between consecutive absorption events beyond the average t and variance σ t , we consider the average waiting times between and within bursts, namely t inter and t intra , respectively. For the simulated events, all four statistical quantities increase as B increases and approaches 1 for a thermal distribution (see inset in Fig.2(b)), while preserving the average intensity of the source (which is obtained by energy conservation and not by the reciprocal of t for thermal fields [14]). It is worth stressing from the inset in Figure 2(b) that σ t grows faster than t . This reflects the faster growth of t inter compared with t intra for larger burstiness B. It can also be observed in the same figure the quick convergence among t and σ t for T/τ c 1, as expected from the trend towards independent random events following a Poisson distribution.
Chan and collaborators recently showed that for the conditions of sunlight absorption, namely ultraweak chromophore-light coupling and intensity, the key parameter determining the quantum dynamics of absorption in photosynthetic systems is the ratio between the chromophores absorption spectral bandwidth ∆ω, and their concomitant emission rate Γ [38], which characterizes the recovery time of the pigments.
Since the coherence time of absorption events is approximated as the inverse of the spectral bandwidth of the absorption (τ c ≈ 2π/∆ω), and the recovery time of the chromophores sets the maximum of the detection time (T max ≈ 1/Γ); the ratio is analogous to the maximum of T/τ c , i.e., ∆ω/Γ ≈ 2πT max /τ c (see Appendic A1 for further discussion). Regarding the recovery time, experiments have shown that excited state absorption occurs in LH complexes with about 90% of the intensity of ground state absorption after ∼0.1 ps [39]. Therefore, at the chromophores absorption spectra (λ LH1 = 875 ± 20nm and λ LH1 = 850 ± 20nm) the coherence time (τ c ∼ 0.1ps) is comparable to their recovery time (1/Γ), reflecting the similar energetic scales of exciton-exciton, exciton-phonon, phononphonon couplings and energetic disorder of this system [38]. In practice, due to the rather low physiological light intensity, doubly excited states will seldom arise in LHs and hence, the photosynthetic absorption set -composed of hundreds of LHs-should be capable to discriminate individual photons within the coherence time. These results indicate that thermal light detection in photosynthetic systems might operate in a regime where the detection time is smaller or comparable to the coherence time, which makes it at least plausible that the temporal correlations present in thermal light may affect absorption statistics.

B. Coupling thermal light absorption with exciton dynamics and charge separation on bacterial vesicles
After photon absorption, the exciton dynamics across the membrane vesicle can be modeled based on the transfer rates between nearest neighbour complexes. Figure 3 summarizes the excitation kinetics and spatial arrangements of the biomolecular complexes LH1, LH2 and RC in purple bacteria. In typical vesicles, the LH2 antenna outnumbers the LH1 complexes, which form an ellipse that encircle the RC protein. Excitation transfer between RCs, LH1s and LH2s occurs from induced dipole transfer on a picosecond time-scale [30], while vibrational dephasing destroys coherences in hundreds of femtoseconds at most [18][19][20][21], resulting in incoherent energy transfer between the RC LH1 LH2 sustained, however, by delocalized excitons over single complexes [13-21, 35, 40, 41]. Transfer rate measures from pump-probe experiments are available for LH2↔ LH1 and LH1↔ RC transfer steps [14][15][16][17]30], whereas transfers LH2↔LH2 and LH1↔LH1, i.e. between the same type of complexes, have been estimated via generalized Förster rates [30,42]. Fluorescence, inter-system crossing, internal conversion and further dissipation mechanisms, have been included with a single conservative lifetime 1/γ D of 1 ns [42]. All transfer rates are depicted in Fig.3(a) and their values are given in Table I. Here t 1,1 is the transfer rate between LH1s, t 2,2 labels the rate between LH2s, t 1,2 (t 2,1 ) denotes the transfer rate from an LH1 to an LH2 complex (viceversa). The transfer rate between an LH1 and its encircled RC is t 1,RC and the respective reverse rate is t RC,1 . Excitonic transfer between LHs has an associated dissipation rate γ D . (b) Spatial arrangement of LH complexes in purple bacteria. Configuration of 44 core complexes (LH1 + RC) (blue circles), and 356 LH2 (green cricles), for the intensity I = 10 −3 photons/(LH·ms) used in simulations.
The formation of a "special pair" (P) in the RC initiates the first electron transfer step P * → P + , followed by several reactions ending up with the reduction of the quinone Q − A or auxiliary quinone Q − B . After P becomes neutral following a rather fast cycle of 1µs, (See Appendix B for details on other photo-chemical transfer times scales involved) [15,18], it can initiate a second electron transfer that results in the reduction of both quinones Q − B Q − B , which triggers the uptake of two intracytoplasmic protons to produce quinol Q B H 2 . The cycle starts over when P is neutral again, and quinol is exchanged by Q B via affinity reactions [15] (see Appendix B for a more detailed explanation). Previous work showed that the interplay between exciton dynamics and the cycling of quinol and the auxiliary quinone affects the quinol output in vesicles [26,28], which were observed to adapt to low and high light intensities [1]. However, the actual output of quinol production in such a model is larger than the turnover of the ATPase molecular rotor discussed in the introduction, and would lead to bacterial death. The subsequent consideration of burst statistics in thermal light did not result in a notable reduction of the quinol output [29], which is the achievement of the present work. Using the photon absorption time-traces exemplified in Fig.2, we initiate the excitonic dynamics with the rates (γ D , t i, j for i or j =LH1,LH2,RC) depicted in Fig.3(a), over a typical photosynthetic vesicle composed of a few hundreds of LHs, illustrated in Fig.3(b). In our simulation we include the charge dynamics occurring within the RCs, therefore, excitations reaching RCs initiate charge separation P * → P + and, crucially, going beyond previous considerations made in references [26,28,29], we also include the recombination of the intermediate metastable state P + QB − . This recombination takes place in about 100 ms to 1 s, a time-scale commensurable with photon absorption waiting times in physiological conditions.
A crucial observation of this work is that the bursted structure of thermal field's statistics is preserved for excitations reaching the RCs. This happens because the transfer of excitons to the reaction center happens in a time scale (∼ps) that is much shorter than the time between bursts (∼s). In detail, the calculation of the burstiness B for the simulated waiting times between absorbed excitations, and for the simulated waiting times between excitations reaching each of the RCs, shows that the bursts characteristic of thermal light drive a bursted arrival of excitations to individual RCs, as shown in Fig.4(a). Since the charge recombination depends on a pairwise exciton arrival to the RC, it might be sensitive to the burstiness of the initial absorption and therefore quinol production may also be influenced by the arrival statistics.
In Fig.4(b) we show how temporal correlations of thermal light can tune the quinol production, as quantified by the quinol yield η = 2Ṅ QH 2 /I tot , whereṄ QH 2 represents the rate of quinol production in the stationary state and I tot is the photon absorption rate of the simulated photosynthetic vesicle. Notice that the quinol yield decreases for higher temporal correlations of the absorbed photons. Thermal light exhibits long waiting time intervals t inter which allow the relaxation of the P + Q − B . In this way, higher correlations in the absorption events result in lower efficiency of quinol production. Temporal correlations arise whenever spatial correlations of thermal light are appreciable. Although spatial correlations are high on the length-scale of bacterial vesicles (cf. Fig 1), their almost constant spatial profile across a full vesicle explains the robust performance of vesicles to changes in the specific geometrical arrangements of LH1 and LH2s (see Appendix C for details).
Under average physiological light intensity I = 10 −3 photons/(LH·ms), which is near to the maximum ATPase capacity, and recombination lifetime t crit between 100 and 1000 ms, the quinol efficiency of the membrane varies across the full range 0 − 90%, depending on the ratio T/τ c (see Fig. 4(b)). For this light intensity, we obtain that the efficiency is strongly dependent on the degree of temporal correlation from absorbed photons. The quinol yield is kept on very small levels for increasing values of t crit the lower the burstiness, as shown in Fig.4(b). Moreover, Fig. 4(c) depicts the quinol production rate relative to the ATPase maximum capacity. For physiological light intensity, this figure shows an interesting crossover for T 0.7τ c where the maximum ATP turnover is achieved (the exact crossover depends on the particular value of t crit ). This last observation identifies a possible photoprotection mechanism. For the same membrane vesicle under much lower light intensity I = 10 −4 photons/LH·ms, this maximum turnover is never reached. In this situation a very low 0% quinol yield is due to the recombination mechanism, independent on the degree of temporal correlations in the absorbed light. The lifetime of charge recombination within the RCs seems to place an important constraint on the physiological light intensity for the survival of photosynthetic bacteria. The sensitivity of a metastable pairwise charge separation to thermal light correlations can therefore help tune quinol production to the ATPase turnover capacity in order to avoid excess acidity and cell damage. The strong photo-protection mechanism is a consequence of the metastable state lifetime and delayed arrival of bursts, and does not result from additional dissipative processes like annihilation. When a burst of photons arrives to the membrane, two photo-excitations have a higher probability to coincide in the same LH structure and annihilate [44]. The dissipation by annihilation, shown in the inset of Fig. 4(c) increases with temporal correlations for the physiological intensity, but its effect remains minimal (a decrease in quinol yield of 0.1−1%) in contrast to the decrease in performance due to the temporal correlations of photo-excitation events shown in Fig.4(b).

II. DISCUSSION
This work studies the role of the statistics of absorption events in the performance of photosynthesis. In our model we have shown that the performance of photosynthesis is affected by the spatio-temporal correlations present in thermal sunlight. The interplay between these correlations and the dynamics within RCs, sets constraints to the light intensities appropriate for bacterial survival, while it provides insight into a pathway for photo-protection, which balances the long intervals between thermal light bursts, and charge recombination taking place in the RCs. Our work underline that not only the average light intensity and the exciton dynamics, but also the statistics of the absorption events and the RC charge dynamics, are relevant aspects for the versatile performance of photosynthetic membranes.
We note that this photo-protective mechanism might operate in different ecological niches, since it encompasses the properties of the exciting field common to all forms of photosynthesis (thermal light), and the dynamics on the RCs which is, except for minor changes, conserved across species. This may suggest a concept of evolutionary fitness of pairwise charge separation for quinol production and synthesis of ATP, which has endured due to the spatio-temporal correlation present in thermal light. Arguably, the actual quantification of the statistics of absorption by a finite absorption linewidth requires further development and provides an interesting perspective for future work, but we should stress that beyond this specific quantification, the presented unified analysis of absorption, excitonic transfer and charge dynamics already provides a glimpse on a regulation mechanism of biological photosynthesis that is driven by the coherence of thermal sunlight.
Besides opening new avenues for further theoretical studies, we stress that our findings lead to predictions that can be experimentally tested with present resources to measure cellular fitness and regulation, e.g., bacterial growth rate, metabolic activity, or gene expression among others. Comparison of fitness indicators for bacteria grown under illumination with attenuated laser and pseudo-thermal sources [13,45] will shed light on the response of these organisms to the spatio-temporal correlations of the incident light. Better fitness indicators for pseudo-thermal light will confirm the significance of the interplay of light absorption statistics and charge recombination, in agreement with our predictions.
Moreover, the statistics of illumination have been shown to improve the rectification stage, i.e. Alternating-to-Direct current conversion in solar antennas [47][48][49]. Rectification nowadays drops the solar antennae collection efficiency (> 80%) to about 0.01% [50,51], but might be improved by the coherent combination of the light detected by different antennas, i.e., the spatio-temporal correlations [47][48][49]. In view of our results, it is an interesting perspective to understand to what extent the spatio-temporal correlations of thermal light could also assist the rectification process in artificial light-harvesting technologies. The absorption of light by a set of receptors is a stochastic process, in which the absorption events are not independent of each other. The quantum photon-interference leads to spatial and temporal correlation between absorption events, affecting the respective photon-counting statistics [2]. Our aim is to introduce the spatio-temporal correlations into the calculation of the probability f (t) for the time elapsed between consecutive detections t. This is used to generate the photon-absorption events required at the beginning of our photosynthesis simulations. That is, our algorithm generates a random number between 0 and 1 for F(t) = t 0 f (t )dt and resolves the corresponding t value, which is the elapsed time since the last detection event until the current one.
In summary the calculation workflow reads as follow: The spatial γ k,l and temporal γ(t 1 , t 2 ) normalized correlation functions are plugged into the generating functional G(s, T ), which in turn allows the calculation of the cumulative probability distribution for inter-event time F(t). Finally F(t) is used to generate the photon-absorption events that will be subsequently converted into photo-excitations in the Montecarlo simulations. Schematically, this is: For the convenience of the reader, we describe below all the quantities involved, while we refer to the Appendix A1 for a review of their derivation for the relevant field; namely a gaussian light field in the far-field approximation, assuming a low absorption bandwidth.
The calculation starts with the spatio-temporal normalized correlation functions (also called coherence functions) γ k,l (t 1 , t 2 ) = Ê + ( r k , t 1 )Ê − ( r l , t 2 ) / √ n k √ n l . HereÊ ± is the electric field operator, k,l refer to a couple of absorbers, r k(l) is the position of the detector labelled k(l), which receive photons at time t 1 (2) ; and n i is the average number of photon-absorptions at each detector. Assuming that the superposition of the light field coming from different points at the source, does not affect its spectral properties at the reception positions -a feature called cross-spectral purity -the spatial and temporal contributions can be separated Γ k,l (t 1 , t 2 ) = √ n k √ n l γ k,l γ(t 1 , t 2 ) [7]. These coherence functions are obtained for quasi-monochromatic light in the far field approximation [1,55] (see Appendix A1 for derivation): The spatial correlations are described by eq. 1, where J(ν k,l ) is the first order Bessel function, and the argument ν k,l = 2π| r k,l |/l c compares the inter-detectors distance with the transverse coherence lenght l c . The spatial correlation function γ k,l quantifies the decay of the spatial correlations as a function of the distance between the detectors k and l. Spatial coherence is relevant for distances inside the transverse coherence length | r k,l | l c . In turn, the transverse coherence length l c = λ D a , relates the average wavelength λ of the field, the distance to the light source D, and its diameter a, as depicted in Fig.1(a) [35]. Within the coherence area, temporal correlations of thermal light are accounted for by a temporal correlation function (eq. 2), whose argument for a stationary field is τ 1,2 = 2π(t 2 −t 1 ) τ c . Similarly to spatial case, γ(t 1 , t 2 ) describes the decay of the correlations as a function of the inter-detections time (t 1 − t 2 ), normalized by the coherence time of the light τ c . The maximum possible value that (t 1 − t 2 ) can get is the detection time T . These temporal and spatial normalized correlation functions have been experimentally well tested, see references [3,4].
Correlation functions (eqs. 1 and 2) are introduced into the statistical description of the light-absorption process by means of the factorial moments generating function G({s i }, T ), which encodes all the statistical information of the process [8][9][10]12]. This formalism accounts for the statistics of the photodetection events recorded by a set of N detectors, each of which we identify as individual LH units in our case. Here we use the generating functional formalism to calculate the photocounting joint probability P(n 1 , n 2 , ..., n N , T ), that n i absorption events occur at the i-th LH (i = 1, ..., N) during the time window T [8, 10, 12, 35, 61]: The probability to absorb n photons in the total set of detectors, regardless of the specific counting record of any individual detector, is then The calculation of the generating functional for Gaussian light leads to [10,12]: such that the spatial and temporal coherence functions enter into the model by means of the eigenvalues j and b k of the spatial and temporal Fredholm equations (see Appendix A2 for the derivation), Finally, it can be demonstrated [14] that the cumulative distribution function F(t) = 1 − P(n = 0, t) (see Appendix A3 for calculation). Therefore, F(t) just depends on the generating functional G(s, T )| s=1 , which in turn depends on the aperture time T (time window for photo-detection), and on the temporal and spatial normalised correlation functions [35].
This photo-detection theory was developed for avalanche photo-diodes where the detection time is a controllable experimental parameter. Therefore the main limitation of this formalism to describe an actual LH complex is the definition of this time window T for photon-arrivals counting. Despite the large amount of existent works regarding the theory of photon-statistics of finite band-width light absorption [11,12,[32][33][34], a fully comprehensive theory applicable to chromophores is yet to be developed. Notwithstanding, the ratio T/τ c 1 or T/τ c 1, sets the limits of Bose-Einstein (maximum correlations) or Poissonian (no correlations) statistics for light absorption, respectively. For T/τ c 1, this framework results in waiting time distributions f (t) that exhibit longer tails than the exponential distribution expected from independent events [14]. The consequence of these slow-decaying tails is a bursted structure [36], with photo-excitation traces that have very long waiting times scattered between bursts of clustered events (typical traces shown in Fig.2(a).
We simulate a full photosynthetic vesicle depicted in Fig.3(b), consisting of 400 LHs (44 core complexes and 356 LH2 antenna complexes), in agreement with realistic stoichiometry for the simulated light intensity [1]. Simulations were performed for two fixed light intensities ( I = 10 −3 photons/(LH·ms) or I = 10 −4 photons/(LH·ms)), scanning different values of T/τ c always verifying that the average rate of photon arrivals is preserved and equals the light intensity.
Although equation 3 allows to obtain the probability of photon arrivals at each LH, the calculation of the photodetection traces for the total set of 400 LH units is a very expensive computational task. Instead of that, we simulate the photo-detection traces for the full membrane (eq. 4) and in a next step we spatially distribute the detected photons in the membrane according to the second-order spatial coherence function of the intensity (see Appendix C for details).

B. Photosynthesis simulations
The generated photo-excitation traces are coupled to a dynamical model of the photosynthetic vesicle, in which the excitons travel through the membrane by individual hopping processes between LH rings until reaching an available RC. This transport process is based upon experimental estimations of excitonic transfer rates and charge transfer dynamics [28]. Here we also include the recombination of the intermediate metastable state P + Q − B , which is described by a stochastic process with associated lifetime t crit after charge separation P * → P + occurred.
After a detection event is simulated, the photo-excitation is located on the corresponding ring (LH1 or LH2), and travels to a reaction center with some probability of being dissipated during the transference. Once the excitation has arrived at the RC, the latter performs the reduction-oxidation cycle described before, and two of these photo-excitations are necessary to form a quinol molecule (Q B H 2 ). In the simulations, when a RC receives the first photo-excitation there is a check for the metastable condition; if the next one does not arrive during the life-time window t crit , the charge is recombined and the exciton is wasted (in this case the next photo-excitation will occupy its place). Otherwise, the reaction center is closed and not available during some time τ RC = 10ms while finishing the process and produces a quinol molecule [26]. Although the turnover of the RCs depends on the acidity, we did not include that dependence and only chose an intermediate rate of quinol production (100 per second at pH 5 [18]), in such a way that the closing time of the reaction centers obeys an exponential distribution with that average throughput. The transfer rates used in our simulations are summarized in Table I, with values taken from pump-probe experiments or Föster calculations assuming intra-complex delocalization, which is in good agreement with experiments [30,63,64].
In the photosynthetic membrane, after a quinol is produced, its two hydrogens are dissociated and used by the ATP-ase enzime to produce ATP. It is worth to mention that as little as 5 RCs already saturate the capacity of the ATP-ase, in the case of slowest RCs turnover equal to 50 quinols per second for a pH = 2 [18]. The mismatch between the overall quinol production and the ATPase turnover determines the bottleneck of the system, restricting the maximum quinol throughput that can be converted into ATP and the incomingphotons rate to 200 photons/s. In our simulations, we use the average intensity 0.4 photons/ms for the whole vesicle ( I = 10 −3 photons/ms·LH), which is processed without limitations by the set of 44 RCs depicted in fig.3(b), whose minimum total capacity is 4.4 photons/ms at pH=2 [18].
Appendix A: Calculations to obtain the coherence functions, the generating functional and the probability distribution of the inter-absorptions time

Spatial and temporal coherence functions
In this section we review the derivation of the second order correlation function for a thermal field in the far-field approximation for small absorption bandwidth. The correlation function is defined as the mean value of the product of the fields at different positions and times [1][2][3][4][5]: where . . . means the average over the system states. E (±) ( r i , t i ) are respectively the positive and negative electric field operators, defined bŷ Here the index µ accounts for polarization, ε k = ck 2 0 (2π) 3 , with c being the speed of the light, k = | k| the magnitude of the wave vector, and 0 the electric permittivity. The creation and annihilation operators of photons in the respective mode areâ † µ, k andâ µ, k . The Green function g k,i that propagates the longitudinal and transverse light modes from the source (in the origin) to the detection point at ( r i , t i ) is g k,i = e i( k· r i −ωt i ) , with ω = ck. Introducing eqs. A2 into eq. A1, the second order correlation of the field, we get: In the last equation polarization has been ignored because it only affects the normalization factor (cf. eq.10.8-2 of reference [1]). Notice that we have introduced a density function η d (k, k ) = (H(k 0 +∆k)−H(k 0 −∆k))(H(k 0 +∆k )−H(k 0 −∆k )) for finite bandwidth absorption in the range ∆k around the central wave number k 0 , where H(k) is the heaviside step function. Bosonic thermal correlation is included through â † kâ k = δ k, k n(k, T), where n(k, T) = 1 e ck/k B T −1 is the number of photons with frequency ω = ck, which come from a thermal source at thermal equilibrium with a reservoir at temperature T. Replacing into the last equation, the correlation function is where we have replaced r i, j = r i − r j and t 1,2 = t 1 −t 2 . The integral in eq. A4 can not be evaluated in a closed form, but reducing the number of detected modes, by assuming small absorption bandwidth and far field detection, a closed expression can be obtained. This approximations are particularly useful in the case of Sunlight absorption on Earth, where light-harvesting complexes absorb on a sharp interval of the light spectra. Lets calculate the spatial case first. In the far field approximation the product k · r i, j is written in spherical coordinates: where φ and θ are the polar and azimuthal angles for the vector indicated in the sub-indexes. Here the azimuth angle has been set to θ r i, j = π/2 for a flat detection surface, and θ s is the azimuthal angle covered by the source. In the far field approximation the angle θ k is taken to be small but not null, allowing to replace sin(θ k ) → θ k , then the correlation function can be written: where we made use of the dispersion relation ω = ck. The angular integral in eq.A7 is well known in the theory of Fraunhofer diffraction for a circular aperture [6], with J 1 (. . .) being the first order Bessel function. Now we normalize the correlation function to write Assuming cross-spectral purity [7], the spatial correlations are not affected by the spectral properties of the field, and can be taken out of the integral (evaluated at a mean frequency ω 0 ) e ω/k B T −1 . Now, applying the limit of small detection bandwidth (∆ω → 0), the factor ω 3 e ω/k B T −1 can be taken out of the integral. In the range of finite absorption bandwidth by LH complexes (λ ≈ 800 ± 50nm), this factor changes ∼ ±1.4%. After this step the last equation reads which (up to a global phase) is Here the arguments are ν i, j = θ s r i, j ω/c = 2π| r i − r j |a/λD and τ 1,2 = ∆ω(t 2 − t 1 ), where we have replaced θ s ≈ a/D (valid for small angles); that is the ratio between the radius of the source a and the distance between the source and the reception surface D.

Factorial moments generating functional of thermal light detection
In general, light has correlations at different orders. Therefore, to describe light detection it is appropriate to use the factorial cumulants generating functional C({s i } , T ), which is defined as the infinite series [8,9]  T 0 dtΓ (m) l,l (t, t). The mth order correlation function Γ (m) can be expressed in terms of lower order correlations, see [10,11]: The second order spatio-temporal correlations of the field (defined in the previous section) are also called the first order intensity correlations, since in eq. A12 Γ (1) k,l (t 1 , t 2 ) = 1 The definition of C(s) contains an infinite series which allows to describe the all-order correlations in the stochastic process, but N detectors can maximally detect N-order correlations Ê + ( r k , t 1 )Ê − ( r l , t 2 ) . For simplicity, we will write from now on Γ k,l (t 1 − t 2 ) for the stationary field. Including all of these in Eq.A11, the factorial moments generating functional gets: ; (A13) where we have defined the matrix D = 1 + m √ (s 1 s 2 . . . s m )(α 1 α 2 . . . α m )Γ (m) and Γ (m) includes spatial and temporal correlations [12]. Because of the Gaussian moments theorem, the D matrix for a Gaussian process will be a four-dimensional operator to include second order spatial and temporal field correlations: Here {κ i } is the set of eigenvalues of the integral equation . Cross-spectral purity permits to assume factorizable spatio-temporal correlations Γ k,l (t 1 − t 2 ) = √ n k √ n l γ(r k,l )γ(t 1 − t 2 ) into spatial γ(r k,l ) and temporal γ(t 1 − t 2 ) correlation functions (defined in the methods section of the main text). The corresponding eigenvalues get κ i = 1 + s i j b k , which depends on the indexes j, k when separating into temporal and spatial equations, and N l=1 α k n k γ(r k,l ) α l n l φ l (t) = b k φ k (t) .
Notice that γ(r k,l ) is the spatial correlation function normalized by the average number of absorbed photons per detector, i.e. each light harvesting ring (LH), such that Γ k,l = √ n k √ n l γ(r k,l ). The average number of photon-absorptions in the k-th detector is obtained by n k = −dG({s i }) ds k | s k =0 = α k N I k T . When α k , N and T are fixed, we will refer to n or I as "intensity" indistinctly.
In order to have correlated absorption events, the set of LHs should be able to discriminate photons inside the coherence volume, i.e. the spatial extension of the system l to be similar or smaller than the coherence length l c , and the detection time T also similar or smaller than the coherence time τ c . The spatial extension of the system is set by the maximum inter-detectors distance (l = max(| r i − r j |)). In a similar way, the temporal extension of the absorption is defined in terms of the inter-detections times (T ∼ max|t i − t j |), where t k is the time at which the detector k is available to absorb a photon. In the special case of a single absorber, one easily sees that the elapsed time for the detector to be available for two consecutive detections (i.e. the time to distinguish photons) is limited by the recovery time of the absorber. In the original experiments with avalanche photo-diodes detectors [5,8,10,13], the detection time was defined as a binning time, which hence gives the timescale needed to resolve the occurrence time of the absorption events. Therefore, T is limited by the time required by the LH complexes to absorb two consecutive photons. Since the capacity of the detection system to distinguish photons inside τ c determines its ability to reproduce the correlated structure of light, the ratio T/τ c allows the continuous tuning of the absorbed excitations from Bose-Einstein (maximum correlations at T/τ c 1) to Poissonian statistics (no correlations at T/τ c 1).

Probability distribution for the time elapsed between consecutive photo-detections
Here we present the procedure described by Rockower in reference [14], to obtain the probability distribution for a random variable t, defined as the time between consecutive photo-detections. We need first to define the quantities used below: tñ is the detection time of theñ-th photon and n(T ) is the number of photo-detection events occurring in a time window T . If the number of events occurring in T is less thanñ, theñ-th detection event occurs after T : On the other hand, there is the cumulative probability function F(T ) = T 0 dt f (t) = P(tñ < T ) and therefore, Reorganizing, we are left with F(T ) = 1 − ñ−1 i=0 P(n = i). We are interested in the time between two successive counts, i.e., the time between the arrival of one photon at the detector and the next one, defined as the "inter-detection time" t. We therefore use the last result withñ = 1 for the time passed before the first detection (cf. [9]eq.6.1): recalling that the Photocounting probability depends on the detection time itself P(n) → P(n, t). Here f (t) is the probability density function for the occurrence of a time-interval t between consecutive detection events.

Appendix B: Thermal light excitation and Charge separation
Photosynthesis captures and processes thermal light in order to produce ATP. Initial light harvesting is followed by subsequent excitonic transfer to RCs, where the electronic excitation of the special pair (P) chromophores within the RC, P → P * , induces a charge separation of the special pair P * → P + + e. This electron e is transferred along an active branch of the RC as shown in Fig.5(a), reducing several cofactors until it reduces a quinone molecule takes place. The oxidised P + , after becoming neutral, can reduce the quinone Q A again, in order to proceed with the uptake of two intracytoplasmic protons , i.e., two Q A reductions take place before a quinol molecule is synthesized. Each reduction has two possible pathways I and II depending on the availability of quinones and cytochromes, as displayed in 5(b)-(c). All the steps and time scales involved are presented in more detail in Fig.6. In a later stage, with the help of the cytochrome complex, quinol Q B H 2 releases these protons and other two protons are pumped from the outside of the vesicle. The electrochemical gradient produced by these four protons is subsequently used by a transmembrane macromolecule called ATPsynthase (ATPase) to produce an ATP molecule. Intermediate metastable states occur between the first Q − A Q B → Q A Q − B and second Q − A Q − B reduction of the quinone pair. In particular, the metastable states P + Q − A or P + Q − B will degrade to PQ A or PQ B due to charge recombination in about 100 ms or 1000 ms, respectively [15][16][17]. The P + Q − B recombination is of special relevance, since the dynamics within the RCs will more often populate this state between consecutive P photo-excitations [15]. While the absorption of photon pairs may help to avoid the decay of the metastable state, quinol production may be inhibited whenever excitations that reach a given RC are delayed longer than charge recombination.  5. (a) Scheme for photo-excitation of the special pair P and subsequent charge transfer to the accessory pigment BChl A , the bacteriopheophytin (Bphe A ) and the quinone Q A ,which take place in the active branch BChl A → Bphe A . If available and neutral, Q B is reduced, but if it was already reduced, the anionic species Q − A Q − B become a two-proton acceptor that results in a Quinol Q B H 2 . The Quinol dissociates from the RC in order to deliver protons to the intracytoplasm and increase the voltage that drives the ATP synthesis [15,18], driven concomitantly by the oxidation of the cytochrome cyt 2 → cyt 3 that restablishes the neutrality of P after photo-excitation. The participation of the main reactants of quinone and cytochrome cycles in the charge transfer dynamics along the active branch is shown in the lower part of (b).
with g (2) (ν) = 1 + |2 J 1 (ν) ν | 2 [1], and D is the diameter of the photosynthetic membrane. Since each type of light harvesting complex absorbs at different frequency, the light absorption is simulated independently for the LH1 and LH2 networks. Then, the two time series are combined to get the absorption signal that enters into the next stages of the photosynthesis simulations. In our numerical implementation, after a photo-excitation is located at the position r i of an LH1 or LH2 complex, the next photo-excitation that goes to the same LH type is placed based upon P(ν). The position for the next photo-excitation r l is selected to lie in a radius | r k − r l |, closest to ν.
The spatial correlations degree is tested by scaling the radius within which the next photon is located (ν → Kν). The PDF plotted in fig.7 shows a profile that flattens out as K → 0 to a uniform distribution. This happens because the detection area gets much smaller than the coherence area (ν 1 ∀i, j) and therefore the spatial coherence becomes very close to its maximum value where the slope of the function goes to zero (see inset in fig.7). Particularly, for the actual size of the membrane (Kν ≤ 0.01), P(ν) is practically constant throughout the full extent of the membrane because of the high degree of spatial correlations. This means, that inside of such a small detection area the arrivals are randomly distributed, but the photocounting statistics are not Poissonian. Note that only when the spatial and temporal correlations are present (l lc and T τ c ), the bunched structure of thermal light manifests (see

Path II
FIG. 6. Steps and time scales followed in paths II and II shown in Figure   5(b). fig.2(a) in the main text).
To work out the effects of spatial and temporal correlations separately, we first calculated the proportion of exciton-pairs converted into quinol, for the 4 cases depicted in Fig.8: 1) No correlations at all (Poissonian incoming signal). 2) Only temporal correlations with random spatial delivering. 3) Only spatially correlated placement of excitations, without temporal burstiness. 4) Spatio-temporal correlations occurring together as it actually is. In all the plots the non-correlation and only-spatial correlation curves overlap, in the same way  as temporal and spatio-temporal correlation curves do. These results point out how the temporal correlations are responsible for the changes in quinol production, since the spatial correlations are very high and almost constant through the membrane. These results corroborate that the spatial assignment of the arrivals does not affect reception statistics (and therefore the quinol efficiency), because the entire vesicle is very small compared with the light's coherence area (the vesicle's area is 10000 smaller than the coherence area), thus behaving like a punctual receptor. Therefore, the internal configuration of detectors is irrelevant for the absorption statistics and for the quinol production.