Probing exciton diffusion dynamics in photosynthetic supercomplexes via exciton–exciton annihilation Open Access

Photosynthesis harnesses and converts solar energy into chemical energy through a coordinated process involving light-harvesting complexes (LHCs) and reaction centers (RCs). In higher plants, photoexcited LHCs generate excitons that transfer to Photosystem II (PSII) RCs, driving water-splitting to produce oxygen and facilitating electron transfer for CO₂ fixation. With a quantum efficiency exceeding 80% under ideal conditions,¹ understanding exciton motion is crucial for optimizing energy transfer efficiency. Although fluorescence lifetime measurements provide valuable insights into exciton dynamics, they are limited in resolving the behavior of individual excitons and the structural design buried in the pigment–protein complex.2–4 Tracking single exciton trajectories and mapping their distributions across PSII remain experimental challenges. The exciton diffusion length (L_D), a key parameter, bridges the gap between exciton lifetimes and their intricate pathways, representing the average distance excitons travel before relaxation or trapping. In photosynthetic complexes, the diffusion length is closely tied to the high quantum efficiency of PSII5–7 and influences energy transfer kinetics, especially under changing membrane morphologies and RC closure.⁸ It also affects the efficacy of nonphotochemical quenching, which protects plants from excess sunlight.⁷ 

The exciton diffusion length can be experimentally obtained using a diffusion model based on exciton–exciton annihilation (EEA), where excitons act as mobile quenchers for other excitons.9–13 EEA has been used in previous studies of photosynthetic systems to obtain the exciton hopping time.14–16 It involves a non-radiative transition of two excitons, where one exciton relaxes to the ground state (S₀), while the other is excited to a higher energy state (S_n) before returning to the first excited state (S₁), as shown in Fig. 1(a).15,17,18 Since the transition occurs rapidly when two excitons reach close separation [which is defined as annihilation radius (R)], exciton mobility is the rate-determining step of the overall EEA dynamics, and the annihilation rate (k_A) represents the reaction rate of a diffusion-limited reaction.14,15,19,20 The exciton diffusion length can be estimated from the annihilation rate using a diffusion model. By analyzing the exciton diffusion length, we learn about the exciton transport characteristics and their connection to the structural arrangement in PSII.

Various techniques, such as transient absorption (TA) spectroscopy, have been utilized to infer diffusion lengths based on the annihilation rate k_A by observing exciton population changes under varied excitation powers.21–23 However, high excitation power introduces higher-order interactions, complicating the analysis of two-exciton dynamics directly related to EEA. Intensity-cycling TA spectroscopy, developed by Malý et al.,¹¹ enables numerical cancellation of higher-order signals [Fig. 1(b)] by adjusting the excitation power to specific fluences. The capability of separating higher-order signals with this technique has been useful for several excitonic systems, including organic dyes,11,24 quantum dots,11,25 LHCII,¹¹ and thylakoid membranes,²⁶ but has not yet been applied to individual photosynthetic complexes.

Despite these advancements, using a diffusion model to analyze the diffusion length still introduces uncertainties due to various assumptions, such as those related to the annihilation radius. To address these challenges, kinetic Monte Carlo (kMC) simulations offer a detailed approach to modeling exciton trajectories as they migrate through the system,28–30 enabling the estimation of key parameters, such as the annihilation rate, diffusion coefficient, and diffusion length.31,32 By comparing simulation results with experimental data, we can gain a deeper understanding of exciton transport. This approach is particularly valuable for well-defined systems, such as the C₂S₂-type PSII-LHCII supercomplex²⁷ [Fig. 1(c), hereafter referred to as C₂S₂], where the availability of a high-resolution structure and a quantum-dynamical model of exciton dynamics, which has been benchmarked against experimental data,2,3,6 enables us to explore various complexities that arise in using EEA to characterize exciton motion.

In this work, we investigated exciton dynamics via singlet EEA in the C₂S₂ supercomplex by combining intensity-cycling TA spectroscopy and kMC simulations. Using intensity-cycling TA, we extracted the fifth-order response signal corresponding to EEA and obtained an annihilation rate k_A,exp of 0.235 ps⁻¹. To explore the assumptions underlying diffusion analysis, we performed kMC simulations to model exciton trajectories using a benchmarked rate matrix2,3,6 and an approach using a distance-dependent annihilation probability. This method allows optimizing the simulated TA profiles to match the experimental results and provides the parameters such as the annihilation radius and system size for describing the two-exciton interaction. Using a 3D normal diffusion model, we obtained a diffusion length L_D,exp = 10.9 nm and a diffusion coefficient D_exp = 1.39 × 10⁻⁷ m²/s, consistent with prior reports (D_exp = 1–5 × 10⁻⁷ m²/s).5–7 While the 3D normal diffusion model allows a direct comparison of the diffusion coefficient with previous studies, trajectory simulations revealed sub-diffusive exciton motion due to strong RC traps and spatial constraints in the chlorophyll (Chl) energy networks. Applying an anomalous diffusion model, we estimated a diffusion length $LD,expsub$ = 9.7 nm, closer to the L_D,sim = 7.4 nm from simulation. Overall, the experimental diffusion length of ∼10 nm suggests efficient exciton mobility across the entire supercomplex for energy transfer. Further analysis of exciton residence time and EEA in subunits showed inhomogeneous distributions, which may relate to proposed quenching sites for photoprotection.33,34 This combined experimental and theoretical approach provides insight into the exciton dynamics and the structure–function relationship within a key photosynthetic complex.

We fitted the third-order response using an exponential decay with two components, $PP(3)=A1⁡exp−t/τ1+A2⁡exp−t/τ2$⁠, which yields an average lifetime τ_avg of 143 ps [Fig. 2(c)]. The fifth-order response shown in Fig. 2(d) was fitted using the fifth-order response function $PP(5)=1−exp−kAtPP(3)$ based on Ref. 11. k_A represents the annihilation rate, which describes the rate of the process where two excitons diffuse, encounter, and annihilate. In a diffusion-limited annihilation, where annihilation happens much faster than diffusion, the annihilation rate k_A can be considered as the reaction rate of a diffusion-limited reaction (further discussion in Sec. II C).19,20 We initially assumed this k_A to be time-independent based on a 3D normal diffusion model where the displacement of excitons changes linearly with time. The fitted annihilation rate k_A is ∼0.235 ps⁻¹, corresponding to an annihilation time constant of 4.3 ps. Compared with LHCII with an annihilation rate of 0.05 ps⁻¹,11,15 the annihilation in C₂S₂ happens at a much faster rate.

To investigate wavelength effects on EEA dynamics, we used a 650 nm pump to primarily excite Chl b while detecting Chl a (Fig. S4). The comparison of PP⁽³⁾ and PP⁽⁵⁾ signals for 650 and 674 nm excitations shows that both wavelengths yield similar annihilation rates (Figs. 2 and S4). This wavelength independence of the annihilation rate aligns with previous LHCII findings,¹¹ where the rate is 0.05 ps⁻¹ whether Chl a or Chl b was excited. Excitation at 674 nm minimizes the Chl b-to-Chl a energy transfer, allowing a more accurate characterization of EEA (more details in the supplementary material, Note 5). Additionally, the measured k_A of the C₂S₂ supercomplex does not depend strongly on the pump power (Figs. S2 and S3), supporting the robustness of our results.

While the measured k_A provides insights into the exciton diffusion characteristics, analyzing it using a diffusion model and extracting the diffusion length in a C₂S₂ supercomplex requires various assumptions: (1) The annihilation radius R needs to be determined to accurately describe the annihilation probability. (2) The system size V, defining the effective volume where diffusion and annihilation event takes place, is needed to convert the microscopic kinetics in two-exciton systems to the macroscopic diffusion length L_D. (3) The nature of exciton diffusion in C₂S₂—whether it is diffusive or sub-diffusive—determines the appropriate diffusion model, such as normal diffusion or anomalous diffusion models. To evaluate these three assumptions and better understand exciton diffusion dynamics, we performed kMC simulations of exciton annihilation, as detailed in Secs. II B–II D.

The kMC simulation of EEA dynamics is comprised of two steps [as illustrated in Figs. 3(a) and 3(b)]: first, simulation of exciton trajectories and second, simulation of the encounter of two excitons. As the first step, the exciton trajectories are generated by stochastically selecting and simulating discrete energy transfer events based on transfer rates,28–30 allowing us to trace the exciton path over time. To align with experimental conditions, the initial excitation in the simulation is distributed over the calculated spectral density of excitonic states within the 675–685 nm range, predominantly exciting Chl a. We adopted the transfer rates previously developed by our group to account for the energy transfer between chlorophylls (Chls).2,3,6 Based on a total of 50 000 exciton trajectories, we statistically analyzed the time that an exciton takes to arrive at the final states, defined as the first passage time (FPT), either by reaching a state that immediately transfers to the charge separation state in RCs or by relaxing non-radiatively to the ground state. The FPT distribution with logarithmic intervals in Fig. 3(c) shows the normalized counts of exciton trajectories after exciting Chl a, which is consistent with the results from analytically solving the master equation (red curve, more details in Sec. IV).

The annihilation radius R is defined as the distance at which the annihilation probability p(r) reaches 99% [1 − exp(−(ɛ/R)⁶) = 99%], corresponding to a diffusion-limited annihilation process. In most studies of diffusion-limited annihilation, the annihilation radius is assumed to be 1–5 nm, representing the mean inter-pigment distance or efficient energy transfer radius,³⁷ rather than being explicitly determined by experimental data or theoretical considerations, which may cause inaccuracies in modeling the annihilation process. In our simulation, we optimized the annihilation radius R by varying it from 0.7 to 6.1 nm to minimize the discrepancy between the simulated and experimental k_A while obeying the constraints imposed by the system size and diffusion length (see Note 2 of the supplementary material for optimization details). The annihilation radius R used in our simulation is 2.7 nm. It is worth noting that while our simulation used a fixed R for all EEA events, its actual value may vary between subunits due to structural inhomogeneity.

Based on the optimized annihilation radius and the distance-dependent annihilation probability, we simulated the annihilation events and analyzed the exciton encounter time, which is defined as the time that two excitons take to travel from the initial excitation location to encounter and annihilate. Figure 3(d) shows the distribution of exciton encounter times of more than 22 000 EEA events. Similar to single-exciton dynamics, the EEA time scale runs from 0.1 ps to a few hundred ps, showing a multitude of exciton transfer pathways. The presence of fast EEA events (around 0.1–1 ps) indicates fast energy transfer between adjacent pigments.

With the simulated single- and two-exciton dynamics, we can reproduce PP⁽⁵⁾ and extract the annihilation rate from the simulation. Since the FPT distribution represents the exciton population that has reached the RC or relaxed, the integral of FPT distribution, ∫dN(t)/dt, can be used to construct the kinetic profile of PP⁽³⁾ involving single-exciton dynamics, as shown in Fig. 3(e), where PP⁽³⁾ = −(1 − ∫dN(t)/dt). The bi-exponential decay fit of the simulated PP⁽³⁾ profile shows an average lifetime of 137 ± 0.2 ps, matching the measured τ_avg = 143 ps for the C₂S₂ supercomplex in our TA experiment and τ_avg = 100–150 ps in a prior report.³⁸ 

Then, we move on to the simulation of PP⁽⁵⁾. Recall that the fifth-order response is given by $PP(5)=1−exp−kAtPP(3)$⁠.¹¹ This equation indicates that singlet EEA contributes to the rise component in PP⁽⁵⁾, i.e., [$1−exp−kAt$], and the decay component represents the single-exciton dynamics described by PP⁽³⁾. In our kMC analysis, the exciton encounter time distribution represents the number of EEA events over a certain period. Therefore, PP⁽⁵⁾ can be constructed by integrating the exciton encounter time distribution over time and multiplying it by PP⁽³⁾. Figure 3(f) shows the simulated PP⁽⁵⁾ profile based on the simulated EEA events using kMC trajectories. By fitting the simulated PP⁽⁵⁾ using the fifth-order response function, we obtained a time-independent annihilation rate of k_A,sim = 0.115 ps⁻¹ from the simulated PP⁽⁵⁾, which is smaller than the value from the experiment (k_A,exp = 0.235 ps⁻¹). This difference may arise from differences in the underlying assumptions between simulations and experiments, such as the distance-dependent annihilation probability or the inhomogeneous diffusion, which will be discussed in more detail in Sec. II C.

The kMC simulation allows us to use the simulated diffusion length L_D,sim and annihilation rate k_A,sim to calculate the system size, which will be used later for analyzing the experimental results to obtain the experimental diffusion parameters (L_D,exp, D_exp). Specifically, L_D,sim is calculated using trajectory displacement. For each exciton trajectory, we define the displacement (σ) as the distance between the density-weighted centers of the initial and current states at a specific time. Then, the diffusion length L_D,sim is represented by the root mean squared displacement of surviving excitons $⟨σ2⟩$ when 63.2% of trajectories have reached the final state, which yields L_D,sim = 7.4 nm and D_sim = 0.67 × 10⁻⁷ m²/s based on the 3D normal diffusion model. We then obtained a system size V of 78.6 nm³/system by inputting D_sim from trajectory displacement and k_A,sim = 0.115 ps⁻¹ from the fitting of the simulated PP⁽⁵⁾ to Eq. (3). Our method of estimating the system size yields a consistent value across different annihilation radii we tested (supplementary material, Note 4). This V value provides a basis for calculating the experimental diffusion coefficients using the same 3D normal diffusion model with k_A,exp = 0.235 ps⁻¹ and R = 2.7 nm. We obtained the experimental diffusion coefficient D_exp = 1.39 × 10⁻⁷ m²/s and diffusion length L_D, exp = 10.9 nm (Table I).

Previous studies, such as those by Amarnath et al.,⁵ have modeled the diffusion coefficient for PSII membranes with different morphologies, including one case where LHCII pools are mixed with PSII supercomplexes (mixed) and another case where LHCII pools are segregated and PSII supercomplexes form arrays (segregated). The diffusion coefficients for the two cases with and without rapid charge separation in the RCs range from 1 to 5 × 10⁻⁷ m²/s.⁵ Notably, the segregated PSII arrays with charge separation exhibit the lowest diffusion coefficient, closely aligning with a previous singlet–singlet annihilation measurement with D = 1 × 10⁻⁷ m²/s.³⁷ This is in good agreement with our experimental diffusion coefficient of 1.39 × 10⁻⁷ m²/s analyzed using a 3D normal diffusion model. Furthermore, the segregated PSII arrays with charge separation also display a diffusion length of around 12 nm similar to the L_D,exp of 10.9 nm obtained from our TA experiment.

Compared with the experimental diffusion length L_D,exp = 10.9 nm analyzed by the normal diffusion model, the diffusion length from kMC simulations L_D,sim = 7.4 nm is slightly shorter. This difference may stem from approximations in the simulation, such as the use of single-exciton energy levels, instead of double-exciton energy levels for modeling two-exciton dynamics. The Pauli exclusion principle and Coulomb repulsion are likely to shift the double-exciton energy away from twice the single-exciton energy.39,40 We accounted for Coulombic interactions when estimating the EEA probability, but the actual exciton annihilation may involve additional mechanisms. If other types of interactions, such as exchange coupling due to orbital overlap, are present in the EEA process, they will alter the EEA rate and the obtained diffusion length. The inhomogeneity of the C₂S₂ supercomplex may also contribute to this difference as variations in EEA probability across different subunits could lead to deviations in the simulated diffusion length. In addition, assuming a 3D normal diffusion when analyzing the experimental dynamics can overestimate the diffusion length, especially in the C₂S₂ supercomplex with spatial constraints and RC traps that result in sub-diffusion. Therefore, in Sec. II D, we will analyze the simulation and experiments using an anomalous diffusion model, where the annihilation rate is time-dependent, k_A(t).

While a 3D normal diffusion model allows us to quantify the diffusion coefficient and compare it with previous reports, such an analysis is based on the assumption of diffusive exciton motion, i.e., the mean squared displacement depending linearly on time. In reality, the exciton motion in a C₂S₂ supercomplex may be sub-diffusive due to the presence of RCs, acting as energy traps that restrict excitation transport. Additionally, structural boundaries and percolated Chl network may also confine exciton motion, making the diffusion process sub-diffusive by limiting the spatial freedom of excitons. This may lead to further discrepancies between experimental and simulation results due to the assumption of a 3D normal diffusion model in fitting the data. To quantify the sub-diffusivity, we analyzed the time dependence of trajectory displacement ⟨σ²⟩(t) ∝ t^α, where α is the anomalous diffusion exponent. A linear dependence indicates a diffusive motion (α = 1), while a sublinear dependence suggests a sub-diffusive motion (α < 1). Instead of a linear dependence, the time dependence of the mean squared displacement in the C₂S₂ supercomplex becomes sublinear [Fig. S6(a)], with α = 0.53.

In the case where exciton motion becomes sub-diffusive, the anomalous diffusion model can be used to analyze the diffusion length. For this anomalous diffusion model, we used a time-dependent annihilation rate k_A(t) = k₀t^α−1, instead of a time-independent annihilation rate, for the PP⁽⁵⁾ fitting equation in Eq. (1).41–43  k₀ is a proportionality constant representing the initial annihilation rate. This time-dependent annihilation rate reflects the time-varying nature of exciton diffusion in a sub-diffusive system. Figure S6 demonstrates that the anomalous diffusion model with α = 0.51 yields better fits of both the experimental and simulated $PP5$ compared to the normal diffusion model. The parameters for the sub-diffusive process were obtained using the same methodologies, substituting the standard diffusion coefficient with anomalous diffusivity (A), calculated as k₀ = 4π(A_m + A_n)R × 2/V. The experimental diffusion length was determined using $LDsub=6Aτavgα$⁠, resulting in an $LD,expsub$ of 9.7 nm, closer to the simulation result of L_D,sim = 7.4 nm. While a normal diffusion model provides analysis that can be compared with prior work, our results show that it will slightly overestimate the diffusion length for the sub-diffusive exciton motion in the C₂S₂ supercomplex (L_D,exp = 10.9 nm).

The sub-diffusive behavior in the C₂S₂ supercomplex can be attributed not only to the energy difference greater than k_BT between charge separation states in the RC and other sites in the C₂S₂ supercomplex but also to the boundary effects due to the relatively small size of C₂S₂ compared to the thylakoid membrane. To demonstrate that the boundary effect limits exciton diffusion, we performed trajectory displacement analysis on the C₂S₂M₂ supercomplex,⁴⁴ which has two additional moderately bound LHCIIs (Fig. S7) and 108 more Chl molecules, in comparison with the C₂S₂ supercomplex. The simulated diffusion length L_D,sim of the C₂S₂M₂ supercomplex increases to 9.7 nm (Fig. S8) with a slightly larger α = 0.56 (Fig. S9) in contrast to the C₂S₂ supercomplex (7.4 nm), suggesting that the size and construction of photosynthetic systems can constrain the exciton diffusion. For even larger systems, such as thylakoid membrane, TA experiments and theoretical modeling have obtained a diffusion length of 30–50 nm as a result of the diffusive exciton transport in the well-connected system.5,26 

To understand how exciton transport and annihilation vary between subunits in the C₂S₂ supercomplex, we summarize the exciton residence time and exciton encounter time across the subunits, including s-LHCII, RC, CP 47, CP43, CP29, and CP26. The residence distribution in Fig. 4(a) shows the percentage of excitons residing in each subunit in different time ranges upon excitation based on trajectory analysis. We observed that a large portion of the living excitons tend to dwell in s-LHCII and CP47, indicating that these subunits serve as primary sites for exciton accumulation. The percentage of living excitons residing in s-LHCII further increases to ∼40% at longer times around 100 ps.

We also analyzed the probability distribution of annihilation events for each subunit at different encounter times based on the EEA trajectories [Fig. 4(b)]. A strong similarity was observed between the two distributions, as expected: the longer an exciton resides in a subunit, the higher the probability of encountering another exciton. Approximately 30% of EEA events occur in s-LHCII over the initial 0–200 ps, and CP47 consistently accounts for 20% of EEA due to the longer residence times of the excitons. Although EEA is typically negligible in nature, we showed that the frequency of EEA is correlated with the exciton residence time in different subunits, demonstrating the potential to use the EEA location as a probe to track the exciton residence time.

The inhomogeneous distribution of the exciton residence time and EEA can reveal the functional roles of the subunits. As a core antenna, CP47 facilitates efficient energy transfer to the PSII RC to ensure optimal photochemical activity. In contrast, s-LHCII primarily functions as a peripheral antenna, capturing and funneling light energy to the core while also serving as an important site for non-photochemical quenching (NPQ).33,34,45 The prolonged exciton residence times and high frequency of EEA events in s-LHCII suggest that quenchers, such as carotenoids, may be more effective when bound to s-LHCII, increasing the likelihood of exciton dissipation through chlorophyll-to-carotenoid energy transfer.26,46,47 These results highlight the synergy between light harvesting and photoprotection, where the exciton motion is regulated by the inhomogeneous protein structure to balance energy transfer efficiency with photodamage mitigation.

However, it is important to note that the measurement and kMC simulations for C₂S₂ were performed at 77 K. Although previous theoretical studies have shown that increasing to room temperature does not significantly alter the energy transfer direction or overall exciton dynamics,³ it can affect the residence time on specific pigments. Additionally, the initial excitation condition for the simulation, designed to match the experiment, represents a small portion of the solar spectrum. Including the full spectrum may change the ratio of excited Chl a and Chl b molecules, thereby affecting the annihilation dynamics. While our analysis provides insights into exciton transport and annihilation dynamics, the connection between annihilation and photoprotection requires further investigation.

Future investigations can further explore the distribution of EEA events across different chlorophylls to elucidate their coordinated functional roles in photosynthesis. For example, our simulations reveal that Chl a612 in s-LHCII consistently exhibits a notable level of EEA, as indicated by the blue dots in Figs. 4(c) and 4(d). In the thylakoid membrane, this Chl a612 is located near lutein 620 (Lut620), a carotenoid previously proposed as a possible quenching site.48,49 This spatial proximity combined with the higher likelihood of exciton residence at this site suggests an increased probability of NPQ. Similarly, Chl a602 in CP26 shows the highest EEA percentage at early times. The low-lying exciton state associated with Chl a602 in CP26 has been shown to facilitate triplet–triplet energy transfer when xanthophyll is bound to CP26, implicating this chlorophyll as a key player in one of the quenching pathways.50,51 In contrast, Chl a602 and Chl a603 in CP29 exhibit the highest number of EEA events at later times, which may be related to the role of CP29 in regulating excitation energy to CP47 and the RC for efficient photosynthesis. Future studies, including simulations at ambient temperature with broad-spectrum excitation, should focus on these site distributions to deepen our understanding of how exciton dynamics reflect the balance between light harvesting and photoprotection in PSII.

The combined experimental and modeling study of exciton motion in the C₂S₂ supercomplex enables us to explore a number of assumptions employed in the application of EEA measurements to synthetic light-harvesting systems. We first explored the dependence of the annihilation rate on the separation of two excitons, which introduces a parameter termed annihilation radius for optimizing EEA simulation (Sec. II B). Next, we examined the approximation of the two-exciton system size for evaluating diffusion-related parameters (Sec. II C). Finally, we evaluated the diffusion length based on the anomalous diffusion model, which accounts for the sub-diffusive exciton motion (Sec. II D).

If a 3D normal diffusion model is used to analyze the experimental fifth-order signal, a value of the L_D,exp is obtained as 10.9 ± 1.6 nm. Our simulation, however, suggests that a sub-diffusive description of the motion is more appropriate with a value of α = 0.51. If this description is used on the experimental data, a $LD,expsub$ of 9.7 ± 0.1 nm is obtained. Our simulation yields a value of 7.4 ± 0.2 nm for L_D,sim using the trajectory displacement method. Since this method uses the actual motion of the excitons in the simulation, which is well fitted by a mean squared displacement proportional to t^0.53, the model value should be compared with $LD,expsub$ = 9.7 nm obtained from the experimental analysis. Given the uncertainties in the analysis of both the experiment and simulation identified in Secs. II B–II D, the agreement between the two diffusion lengths seems reasonable. A diffusion length of 7–10 nm is consistent with the notion that an excitation can explore the entire supercomplex no matter where its initial location is. If one of the two RCs is closed by an earlier excitation, an exciton initially created close to that RC still has a good probability of reaching the other, open, RC. Moreover, we showed that boundary effects from protein structural constraints should also be considered a factor in sub-diffusive behavior. When considering larger systems, such as the thylakoid membrane, the work of Bennett et al.6,7 and Amarnath et al.⁵ suggests that in such systems where excitons have access to much larger numbers of chlorophylls, the motion is more likely to be characterized by diffusive rather than sub-diffusive motion.

Analysis of the simulation trajectory data indicates that the length of time an exciton spends in the various subunits that make up the C₂S₂ supercomplex is highly inhomogeneous, showing that EEA probabilities will also be distributed. This analysis also provides insight into the putative sites for non-photochemical quenching by suggesting locations where activatable quenchers are likely to be most effective.

Future studies should develop more quantitative, structure-based models of exciton–exciton annihilation so that the annihilation probability in the simulation can be input rather than being a parameter. Overall, this comparatively simple system demonstrates the value of higher-order spectroscopy for characterizing exciton motion with its combination of both high temporal and spatial resolution and exposes the current gaps in our knowledge of exciton–exciton annihilation.

To prepare the samples, all procedures were carried out in the dark to minimize light exposure. PSII-enriched membranes were extracted from spinach leaves, following established protocols52,53 with slight modifications.54,55 For the isolation of the C₂S₂-type PSII-LHCII supercomplex, PSII-enriched membranes (0.5 mg of Chl/ml) were solubilized in a 1.0% (w/v) α-DDM (n-dodecyl-α-D-maltopyranoside, Anatrace) solution, prepared in a buffer containing 25 mM MES-NaOH (pH 6.0), and incubated on ice for 30 minutes. The solubilized membranes were then centrifuged at 21 000 × g for 5 min at 4 °C. The resulting supernatants were loaded onto sucrose gradients, consisting of 2.1 ml layers of 0.1, 0.4, 0.7, 1.0, and 1.3M sucrose, with 0.5 ml of 1.8M sucrose at the bottom, in ultracentrifuge tubes (14 × 89 mm², Beckman Coulter). The gradients were prepared with a buffer containing 0.03% (w/v) α-DDM, as described above. Ultracentrifugation was carried out at 154 300 × g for 24 h at 4 °C using a SW 41 Ti swinging-bucket rotor (Beckman Coulter). After centrifugation, the separated bands were collected dropwise from the bottom of the tube. The fraction containing the C₂S₂-type PSII-LHCII supercomplex was concentrated using a 100 K MWCO centrifugal filter. The concentrated sample was then diluted with a buffer containing 25 mM MES-NaOH (pH 6.0), 10 mM NaCl, 3 mM CaCl₂, 400 mM sucrose, and 0.03% α-DDM prepared with D₂O. Finally, the concentrated C₂S₂-type PSII- LHCII supercomplex was flash-frozen and stored at −80 °C for later TA measurements.

The transient absorption spectroscopy used a common pump–probe scheme. A Ti/sapphire regenerative amplifier (RegA 9050, Coherent) was pumped by a diode laser (Verdi V-18, Coherent) and seeded by a Ti/sapphire oscillator (Vitara-T, Coherent), generating 800 nm pulses with a repetition rate of 250 kHz. The pulse was modulated by an external stretcher/compressor and then split into pump and probe beams by a beam splitter. For the pump beam, an optical parametric amplifier (OPA 9450, Coherent) generated pulses centered at 674 nm with a 50 nm full-width-at-half-maximum (FWHM) and then passing through a 675 nm longpass filter for exciting the Q_y bands of Chl a. The pump beam was compressed by prism pairs to an FWHM of 82 fs. The pump fluence was adjusted between 2.5 and 10 nJ by a neutral-density filter wheel. For the probe beam, a visible continuum was generated by a 1 mm sapphire crystal and filtered by a 700 nm shortpass filter. The pump and probe pulses passed through the cryostat (OptistatDN, Oxford instrument) where the C₂S₂ supercomplex sample was cooled down to 77 K (path length = 400 μm, OD = 0.8 at 674 nm). The pump and probe beam sizes were 160 and 80 μm in diameter, respectively, and the beams were overlapped at the sample at the magic-angle (54.7°) polarization. The systematic error of pump fluence due to the Gaussian beam profile, comparing the average fluence to the edge, is ∼20%. After passing through the sample, the probe was filtered by a polarizer to remove the scattered pump beam. A monochromator (SpectraPro 300i, Acton Research Corp.) was tuned at 680 nm to select the probe wavelength. The exit pulses were collected by a diode detector (DET10A, Thorlabs), generating analog signal input to a lock-in amplifier (SR830, Stanford Research), which synchronized the pump–probe signals with a chopper positioned in the pump beam path. For a set of intensity-cycling TA measurements in Fig. 2, we set the pump fluences to be 2.5, 7.5, and 10 nJ. The fifth-order signal captures EEA, rather than saturation and cascade effects, considering the number of Chls in the C₂S₂ supercomplex and the fact that the relative intensities of PP⁽³⁾ and PP⁽⁵⁾ do not show significant differences for pump powers of 2.5/7.5/10 and 6/18/24 nJ.

The kinetic rate matrix used for the kinetic Monte Carlo simulations has been reported by Leonardo et al.² and Yang et al.³ Briefly, the approach used by Bennett et al.⁶ to develop the structure-based model was applied similarly here to describe the excitation energy transfer network in the PSII supercomplex. Based on semi-empirical Hamiltonians of the subunits in the PSII supercomplex, excitonic states were separated into domains according to electronic couplings and degrees of delocalization, allowing the classified calculation of the energy transfer rate using generalized Förster theory for inter-domain and modified Redfield theory for intra-domain energy transfer. The TrEsp method⁵⁶ was employed to evaluate the interprotein coupling based on the cryo-EM protein structure of the C₂S₂ supercomplex (PDB: 3JCU)²⁷ and the C₂S₂M₂ supercomplex (PDB: 5XNL).⁴⁴ The exact parameters used can be found in Refs. 2 and 3. To match the experimental conditions, the rates for the C₂S₂ supercomplex were calculated with a temperature set as 77 K. The initial excitation location is calculated based on the simulated absorption probability (absorption lineshapes) of all excitonic states within the 675–685 nm range, which predominantly leads to Chl a excitation. The rates for the C₂S₂M₂ supercomplex were calculated at room temperature. To improve simulation accuracy, exciton in the lowest excitonic state in the RCs was assumed to undergo charge separation instantaneously rather than to leave the state.⁶ A non-radiative relaxation with 1 ns lifetime was set for all excitonic states.

The kMC simulation was performed using the energy transfer rate matrices under the exciton basis. The master equation for the time dependence of the exciton population is determined by the rate matrix (K): dP/dt = KP, where P(t) is the exciton population at time t, and K is a rate matrix describing the energy transfer between states. The overall P(t) of non-final states provides the time-dependent population decay of excitons. Initially, an excitonic state was randomly populated according to the absorption probability of the excitonic states in the range of 675–685 nm. The kMC trajectory can stop at any of the final states, including the lowest-energy state in the RCs and the ground state via nonradiative relaxation. We generated a total of 50 000 trajectories for each kinetic model to process the trajectory analysis. To ensure the validity of the 50 000 trajectories, the FPT distribution obtained from the trajectories is compared with that obtained from eigendecomposition. The results generated by the two approaches are consistent.

To simulate the two-exciton dynamics in a C₂S₂ supercomplex, we randomly picked two independent trajectories from the kMC trajectory pool, with the constraint that the initial state of the trajectories are different excitonic states. During the simulation, the excitons may encounter each other and undergo EEA before reaching their final states. The probability of annihilation is defined based on the Coulomb interaction between states as a function of distance (see Note 1 of the supplementary material for detailed definition). We introduced a parameter, annihilation radius, which defines the range of 99% annihilation probability for optimizing EEA simulation. If EEA occurs between two trajectories, we randomly remove one exciton from the simulation, while a surviving exciton continues along its trajectory. The surviving exciton undergoes rapid S_n → S₁ relaxation within a few femtoseconds after receiving the energy from the other exciton,17,18 allowing us to assume that EEA causes no time delay and the surviving exciton can maintain its trajectory to reaching its final state. The exciton encounter event (i.e., occurring time and states) of each EEA simulation was recorded for processing further statistical distribution analysis.

The supplementary material includes the assumptions used in the kMC simulation, transient absorption results of the C₂S₂ supercomplex under different excitation conditions, and simulation results of the C₂S₂M₂ supercomplex.

This work was supported by the U.S. Department of Energy, Office of Science, Basic Energy Sciences, Chemical Sciences, Geosciences and Biosciences Division, under field work proposal 449A. We thank Professor Doran I. G. Bennett Raccah for sharing the code for the rate matrix calculation.

The authors have no conflicts to disclose.

K.Z. and T.-Y.L. contributed equally to the work.

Kunyan Zhang: Conceptualization (equal); Formal analysis (equal); Investigation (equal); Writing – original draft (equal); Writing – review & editing (equal). Tsung-Yen Lee: Conceptualization (equal); Formal analysis (equal); Investigation (equal); Software (equal); Writing – original draft (equal); Writing – review & editing (equal). Shiun-Jr Yang: Software (equal); Writing – review & editing (equal). Trisha Bhagde: Investigation (equal); Writing – review & editing (equal). Masakazu Iwai: Resources (equal); Writing – review & editing (equal). Graham R. Fleming: Conceptualization (equal); Funding acquisition (equal); Writing – original draft (equal); Writing – review & editing (equal).

The data supporting the findings of this study are available within the article and supplementary material or can be obtained from the corresponding author upon reasonable request. The code used for generating the rate matrices is available at https://doi.org/10.5281/zenodo.13346121. The code used for trajectory analysis is available at https://doi.org/10.5281/zenodo.14790382.