Molecular Origin of Correlated Bath Effects in Photoinduced Charge Transfer Dynamics in Polar Solvents

TL;DR

This paper explores how environmental correlations in polar solvents influence charge transfer dynamics after light absorption.

Contribution

A geometric metric is proposed to quantify bath correlations in multistate nonadiabatic dynamics.

Findings

Correlated baths synchronize energy fluctuations to facilitate competing reactions.

Anticorrelated baths suppress one reaction when favoring another.

Molecular conformation and dipole moments modulate energetic correlations.

Abstract

Nonadiabatic dynamics
in the condensed phase often involve correlated
environments shared by multiple electronic states, challenging the
traditional isolated bath assumption. We investigate these effects
using the multistate harmonic (MSH) model and atomistic Hamiltonian
applied to photoinduced charge transfer in a trimer consisting of
a methylperylene donor and two tetracyanoethylene acceptors dissolved
in a polar solvent. We propose a geometric metric based on the angular
relationship of reorganization energies between transitions sharing
an initial state to quantify bath correlation. Our analysis identifies
distinct regimes: a correlated bath where synchronized energy gap
fluctuations facilitate competing reactions, and an anticorrelated
bath where fluctuations favoring one reaction suppress the other.
These energetic correlations are modulated by molecular conformation
and charge distribution, specifically through changes in dipole moments
and solvent-accessible surface area. This study provides a connection
between the energetic perspective of environmental correlations and
the molecular details governing nonadiabatic dynamics in polar solvents.

Figures

Full text

Nonadiabatic dynamics plays a fundamental role in the photophysical and photochemical processes governing solar energy conversion, ranging from artificial photovoltaic cells to natural photosynthetic systems. ?−? ? ? ? ? ? ? ? The complex, intertwined nature of electronic and nuclear degrees of freedom (DOF) in these systems requires approaches that go beyond the widely adopted Born–Oppenheimer approximation underlying modern quantum chemistry. Consequently, there is a rich history of developing nonadiabatic dynamics methods for complex molecular systems in the condensed phase, aiming to elucidate the mutual influence of electronic and nuclear DOF within effective model or atomistic molecular Hamiltonians. ?−? ? In standard system-bath models, the relevant electronic states typically constitute the basis for the system Hamiltonian, while the nuclear DOF consisting of intramolecular and intermolecular vibrations form that for the bath Hamiltonian, with system-bath coupling mediating the nonadiabatic dynamics. ?−? ? ? ? ? ? A widely used example is the Frenkel exciton model for excitation energy transfer (EET), ?,? which traditionally relies on the “isolated bath assumption”: that each chromophore couples to its own local environment, and that these local baths are statistically independent. ?−? ? ? ? ? ?

However, the isolated bath assumption breaks down when chromophores are in close proximity, sharing a common solvent environment, or when charge transfer (CT) states form between moieties embedded in the same solvation environment. ?−? ?,?,?−? ? In such scenarios, one expects significant shared or correlated bath effects to influence the nonadiabatic molecular dynamics (NAMD). Theoretical treatments of solvation dynamics have long recognized that solvent motions are not merely independent local fluctuations but organize into collective modes. Seminal work by Stratt and collaborators utilizing instantaneous normal mode (INM) theory demonstrated that the solvent shell possesses specific symmetries that can selectively couple to solute electronic transitions. ?−? ? ? These studies highlighted that the solvent response is inherently collective, with solvent motions capable of promoting or suppressing specific reaction pathways depending on how the solvent’s modal symmetry aligns with the electronic transition dipoles.

The dynamical consequences of such environmental correlations have been a subject of significant theoretical inquiry. Model studies by Strümpfer and Schulten? demonstrated that correlated bath fluctuations can modify energy transfer efficiency in photosynthetic complexes. Similarly, Dutta et al.? showed that spatial and temporal correlations critically shape excitonic spectral line shapes. However, characterizing the magnitude of these effects in realistic systems remains nontrivial. For instance, in atomistic molecular dynamics simulations of the light-harvesting II complex, Olbrich and Kleinekathöfer? found that the cross-correlations between the excitation energy fluctuations of different pigments were negligible. This contrast suggests that the manifestation of bath correlation is highly sensitive to the molecular geometry, the nature of the electronic states, and the specific correlation functions employed to quantify the environment. For instance, a recent study on the EET of two nonfullerene acceptor Y6 molecules in chloroform demonstrated that the bath transitions from correlated to uncorrelated as the intermolecular distance increases.? Similarly, photoinduced CT dynamics in a carotenoid-porphyrin-fullerene (donor-bridge-acceptor) triad dissolved in tetrahydrofuran reveal that the locally excited ππ* state and the various CT states share a common solvent environment, leading to conformation-dependent correlated or even anticorrelated bath effects. ?,?,?

It appears that correlated bath effects in nonadiabatic processes are dictated not merely by physical proximity or shared environments, but fundamentally by how this shared environment couples to the multiple electronic states. Broadly, bath correlation describes a regime where the environments of distinct electronic states cannot be treated independently. Yet, the precise definition of a “correlated bath” varies in the literature and often lacks a clear, quantitative molecular picture. In this Letter, we focus on establishing a rigorous description of correlated bath effects in nonadiabatic dynamics and clarifying the underlying molecular physical picture.

The recently developed multistate harmonic (MSH) model ?,?,? provides an ideal framework for this task, allowing for an unambiguous definition of the degree of bath correlation. In contrast to standard excitonic models that partition the environment into independent local baths, the MSH framework employs a “globally shared bath”, where a single set of global normal modes couples to all electronic states simultaneously. This globally shared bath is parametrized from realistic anharmonic atomistic simulations to systematically satisfy all pairwise reorganization energy constraints. This globally shared bath is rigorously defined within an extended (F – 1)-dimensional nuclear space for an F-state system. The MSH model’s global bath naturally encompasses the full spectrum of environmental correlations, ranging from correlated to anticorrelated bath regimes, and recovers the uncorrelated or isolated bath limit as a specific critical case.?

To quantify the bath correlation, one must consider at least three electronic states and two distinct electronic transitions, for example, |1⟩→|2⟩ and |1⟩→|3⟩. In the extended nuclear space, the square roots of the pairwise reorganization energies for these transitions form a triangle, allowing the correlation to be defined by the angle θ_23_ between the transitions:

$$[eqn]$$

Here, E _ r _ ^(12)^, E _ r _ ^(13)^, and E _ r _ ^(23)^ represent the reorganization energies for the transitions |1⟩↔|2⟩, |1⟩↔|3⟩, and |2⟩↔|3⟩, respectively. These values are determined directly from the energy gap fluctuations in all-atom simulations via:

$$[eqn]$$

In this expression, σ_ XY _ ^2^ = ⟨U _ XY _ ^2^⟩ – ⟨U _ XY _⟩^2^ is the variance of the energy gap U _ XY _ = V _ X _ – V _ Y _, and k _ B _ T is the thermal energy with Boltzmann constant k _ B _ and temperature T. For EET, |1⟩ is typically selected as the ground state while |2⟩ and |3⟩ are locally excited states; for photoinduced CT, |1⟩ may be a locally excited state while |2⟩ and |3⟩ are two CT states between different molecules or of distinguished CT character. Importantly, the MSH model preserves the exact reorganization energy values obtained from the all-atom Hamiltonian.?

To provide a physical interpretation of this geometric metric, it is instructive to consider the distinct regimes of bath correlation. Mathematically, cosθ_23_ corresponds to the normalized cross-correlation coefficient between the energy gap fluctuations of the two transitions. When θ_23_ = 90°, the cosine is zero, implying that the fluctuations of the energy gaps U 12 and U 13 are statistically uncorrelated. In the MSH framework, this signifies that the equilibrium nuclear displacements for the two transitions are orthogonal in the extended-dimensional nuclear space; physically, this corresponds to a scenario where the bath motions coupled to one transition are entirely independent of those coupled to the other, as often observed in EET between spatially well-separated chromophores.? Deviations from orthogonality indicate bath correlation effects: an acute angle (θ_23_ < 90°) implies a positive correlation, where solvent fluctuations that stabilize state |2⟩ (relative to |1⟩) simultaneously stabilize state |3⟩. Conversely, an obtuse angle (θ_23_ > 90°) implies a statistical anticorrelation, where solvent motions favoring one transition energetically penalize the other. We emphasize that this classification is transition-specific; a globally shared bath may be correlated for one pair of transitions while being anticorrelated for another within the same MSH model.

To this end, we propose using the angle θ_23_ = θ_CT1,CT2_ as a quantitative metric to categorize the degree of bath correlation for the transitions of interest (e.g., EX → CT1 and EX → CT2). Specifically: (1) θ_23_ < 90° indicates a correlated bath; (2) θ_23_ > 90° indicates an anticorrelated bath; and (3) θ_23_ = 90° indicates an uncorrelated (independent) bath.? This geometric definition highlights that the degree of bath correlation is an intrinsic property of the system’s energy gap fluctuations, dependent on the interplay of all three pairwise reorganization energies.

In this work, we investigate correlated bath effects in nonadiabatic dynamics by focusing on a donor–acceptor-acceptor (DA_1_A_2_) molecular trimer. This system comprises one methylperylene (MPe) donor and two tetracyanoethylene (TCNE) acceptor molecules dissolved in acetonitrile (ACN) solvent at 300 K. ?,?−? ? ? Our approach goes beyond qualitative descriptions of collective solvent modes by providing: (a) a rigorously defined, model-independent geometric metric (θ_23_) for classifying bath correlation regimes; (b) the integration of this metric into the MSH formalism, which preserves atomistic bath correlations within a reduced-dimensionality framework; (c) a concrete molecular interpretation of the metric that extends beyond abstract collective solvent modes; and (d) a systematic model construction strategy that disentangles bath correlation effects from those driven by the system Hamiltonian. The intermolecular photoinduced CT process initiates with the formation of an excitonic (EX) state localized on the MPe molecule (D*A_1_A_2_). This is followed by electronic transitions to two distinct CT states, denoted as CT1 and CT2, corresponding to charge transfer from the MPe donor to the first (D^+^A_1_ ^–^A_2_) and second (D^+^A_1_A_2_ ^–^) TCNE acceptor molecules, respectively.

Figure(a1–c1) depicts three selected conformations of the trimer: the gas-phase optimized structure (R0); a modified structure (R1) where the first TCNE molecule (A_1_) is displaced 3 Å away from the donor relative to R0; and a sandwich-like structure (R2) where the two TCNE acceptors are positioned on opposite sides of the MPe donor. Our analysis centers on the transitions of interest, EX → CT1 and EX → CT2, with the corresponding reorganization energy relationships illustrated in Figure(a2–c2). Table lists the reorganization energies for all state pairs across the three conformations. Based on the angular parameter θ_23_ listed in Table, which is computed using these reorganization energies, it is evident that conformations R0 and R1 exhibit a correlated bath, whereas conformation R2 displays an anticorrelated bath.

To describe the nonadiabatic dynamics, we employ the general F-state Hamiltonian

$$[eqn]$$

where the diagonal terms represent the nuclear Hamiltonian of the jth electronic state (j = 1, ···, F), given by H _ j _ = P ^2^/2 + V _ j _(R). Here, R and P are the mass-weighted nuclear positions and momenta, and V _ j (R) is the potential energy surface (PES) of the jth electronic state. The off-diagonal terms Γ jk _ (j≠k) describe the electronic coupling between states j and k. For the all-atom multistate Hamiltonian, we construct the PES of the F states using a modified general Amber force field (GAFF), where the bonded interaction parameters are kept the same according to GAFF. However, the atomic charges and excitation energies of different electronic states are parametrized via the electronic structure calculations, as detailed in the Supporting Information. This potential of the jth state is expressed as?

$$[eqn]$$

Here, V _ j _ ^FF^(R) represents the force field (FF) potential energy of the jth state. The term W _ j _(r ^tri^) = E _ j _(r ^tri^) – V _ j _ ^FF^(r ^tri^) acts as an energy correction that incorporates the excitation energy E _ j _(r ^tri^) obtained from electronic structure calculations for the trimer geometry r ^tri^, and the all nuclear positions R includes the solute trimer configuration r ^tri^ and the solvent configuration. This correction explicitly removes the double-counted intramolecular trimer energy in the FF term V _ j _ ^FF^(r ^tri^), which has been already included in the excitation energy E _ j _(r ^tri^) evaluated at a higher quantum mechanical level. For more details, refer to refs ? and ?.

Alternatively, within the MSH model Hamiltonian, the PES of the jth state takes the form?

$$[eqn]$$

where the nuclear DOF span an extended (F – 1)-dimensional space. The coordinate R _ a,i _ corresponds to the ith normal mode with frequency ω_ i _ (i = 1, ···, N) within the ath nuclear subspace (a = 1, ···, F – 1). The parameter S _ i _ ^(aj)^ defines the equilibrium shift, subject to the constraint that S _ i _ ^(aj)^ ≡ 0 when a ≥ j. The set of discretized frequencies {ω_ i _} is determined by the spectral density, which is derived directly from the energy-gap time correlation function calculated in the all-atom simulations.? A fundamental requirement of the MSH model is that it must simultaneously satisfy all pairwise reorganization energies derived from all-atom simulations. This constraint is expressed as

$$[eqn]$$

The energy minima parameters, {ε _ j _}, and their differences, defined as the reaction free energies ΔE ^(XY)^, are determined from the averaged energy gap and the reorganization energy via

$$[eqn]$$

The MSH model Hamiltonian can be cast into system-bath form, H = H S + H B + H BS. The system Hamiltonian is given by

$$[eqn]$$

and the sum of bath and system-bath coupling terms is expressed as

$$[eqn]$$

where

$$[eqn]$$

represents the shifted harmonic potential surfaces in the extended nuclear space without the vertical shifts ε _ j _.

To investigate the contributions of different Hamiltonian components to the nonadiabatic dynamics, we designed three distinct parametrization schemes for both the all-atom description and the MSH model construction.

• Scheme 1 [S1]: This scheme represents a fully consistent, conformation-specific parametrization.? For each conformation (R0, R1, and R2), the all-atom Hamiltonian is constructed using electronic couplings ({Γ_ jk _}), excitation energies ({E _ j _}), and multistate force field (FF) parametersincluding atomic charges and energy corrections ({W _ j _}) that are calculated individually for that specific geometry. Subsequently, the MSH model parameters for each conformation, including electronic couplings, energy minima ({ε _ j }), mode frequencies ({ω i _}), and equilibrium shifts ({S _ i _ ^(aj)^}), are determined directly from the corresponding all-atom simulations.
• Scheme 2 [S2]: This scheme is designed to isolate the effect of solute geometry on the bath while keeping electronic parameters fixed to the reference state. For constructing the all-atom Hamiltonians of conformations R1 and R2, we utilize the electronic couplings, excitation energies, FF parameters, and atomic charges derived from the R0 conformation. However, the trimer geometry itself is changed to R1 or R2, which requires recalculating the energy corrections to the force fields. Correspondingly, the MSH models for R1 and R2 under Scheme 2 retain the electronic couplings and excitation energies from the R0 electronic structure calculations, while reparameterizing the energy minima and the combined bath/system-bath coupling terms (H BBS) based on all-atom simulations of the R1 and R2 conformations within this mixed parametrization.
• Scheme 3 [S3]: This scheme effectively combines the electronic system of the reference state with the full environmental bath of the target state. To construct the MSH models for conformations R1 and R2, we employ the system Hamiltonian H S (including electronic couplings and energy minima) taken directly from the R0 MSH model. In contrast, the bath and system-bath coupling terms (H BBS) are adopted from the fully consistent MSH models [S1] of the respective R1 and R2 conformations.

Consequently, conformation R0 is described solely by Scheme 1, whereas conformations R1 and R2 utilize Schemes 1–3. Scheme 1 is the standard consistent approach for parametrizing both all-atom and MSH Hamiltonians, where electronic structure and molecular dynamics calculations are derived directly from the specific geometry of each conformation, particularly regarding the charge distributions of different electronic states. Conversely, Scheme 2 for conformations R1 and R2 employs the electronic couplings, excitation energies, and force field parameters from conformation R0. This design allows us to isolate the influence of the trimer’s geometric changes on the nonadiabatic dynamics. Scheme 3 acts as an intermediate bridge between Schemes 1 and 2. Comparing the R0 reference with Scheme 3 for a given conformation (R1 or R2) highlights the differences arising from the bath and system-bath coupling parameters, as their system Hamiltonians are identical. Furthermore, comparing Schemes 1 and 3 for the same conformation (R1 or R2) isolates the impact of the system Hamiltonian, specifically the electronic couplings and excitation energies, as the bath parameters remain constant. We included Table S2 in the Supporting Information to summarize the main features of these schemes.

Table details the all-atom parameters for the three trimer conformations under Schemes 1 and 2, including the reorganization energies, and the magnitudes of dipole moment changes accompanying electronic transitions. Table summarizes the MSH model parameters for these conformations across Schemes 1 through 3, alongside the Marcus rate constants for all potential excited-state electronic transitions. Table lists the solvent-accessible surface area (SASA) of the trimer segments calculated via the Shrake-Rupley approach,? and the minimal intermolecular distances between the donor and the acceptors. For a generic transition X → Y, the Marcus rate constant is defined as ?−? ?

$$[eqn]$$

and the corresponding activation energy is given by

$$[eqn]$$

It is important to note that these Marcus rate constants are utilized here primarily as a diagnostic tool to interpret the dominant kinetic pathways and activation barriers, rather than as quantitative predictions of the full nonadiabatic dynamics, which are rigorously simulated using NAMD approaches. Figure visualizes the Marcus parabolas and the associated rate constants for the three conformations. ?,? These rate constants reveal two primary kinetic pathways for the photoinduced CT process:

$$[eqn]$$

This identification stems from the observation that only the rate constants k 1→2, k 1→3, and k 2→3 exhibit significant magnitudes.

The first question we address is identifying which molecular properties significantly influence the reorganization energies and, consequently, the degree of bath correlation, θ_23_. Tables and ? suggest that the change in the solute’s permanent dipole moment upon transition plays a primary role, while the solvent-accessible surface area (SASA)? plays a secondary role. For example, comparing Scheme 2 against Scheme 1 for the same conformation (R1 or R2) reveals overestimated reorganization energies for the EX → CT2 (E _ r _ ^(13)^) and CT1 ↔ CT2 (E _ r _ ^(23)^) transitions, alongside an underestimated reorganization energy for the EX → CT1 (E _ r _ ^(12)^) transition. As defined in eq, the reorganization energy is determined by the thermal fluctuations of the energy gap between the two PESs, which microscopically corresponds to solvent reorganization upon the electronic transition. Therefore, a more significant rearrangement of solvent molecules results in a larger reorganization energy.

In the context of photoinduced CT within a polar solvent, a substantial change in the electrostatic field induced by the electronic transition of the solute drives a significant rearrangement of the surrounding acetonitrile (ACN) molecules. The magnitude of the solute’s permanent dipole change upon transition, denoted as |Δμ _ XY _|, serves as the primary metric for characterizing these electrostatic interactions. This metric aligns with the trends in reorganization energy observed between Scheme 1 and Scheme 2 in Table. Specifically, |Δμ 12| decreases, whereas |Δμ 13| and |Δμ 23| increase for both the R1 and R2 conformations. We also observe that for a fixed trimer conformation, such as the comparison between R1 [S1] and R1 [S2], a larger mean absolute difference in atomic charges generally corresponds to a larger |Δμ _ XY _| (see Table S1).

Furthermore, when the set of atomic charges is held constant and only the solute conformation varies, the dipole change upon transition remains strongly correlated with the reorganization energy trends. For instance, the transition from R0 to R1 [S2] exhibits a significantly increased |Δμ 12|, which implies an increase in E _ r _ ^(12)^. This shift is directly attributable to the displacement of the TCNE A_1_ molecule 3 Å further away from the MPe donor. However, distinct behaviors emerge when the dipole change is identical. Both R0 and R1 [S2] possess the same |Δμ 13| values, yet the reorganization energy E _ r _ ^(13)^ is larger in R1 [S2]. In such cases, the SASA serves as a secondary descriptor where a larger SASA generally leads to a larger reorganization energy. Since the transition EX →CT2 involves the MPe donor and the TCNE A_2_ acceptor, the observed increase in SASA­(D, A_2_) from R0 to R1 [S2] in Table aligns with the corresponding trend in reorganization energy.

Figure shows the photoinduced CT dynamics over the first 2 ps for the three trimer conformations, employing both all-atom and MSH model Hamiltonians parametrized via Schemes 1 and 2. The NAMD simulations were conducted at 300 K using the classical mapping model (CMM) with a zero-point-energy parameter of γ = 0.309. ?−? ? The CMM approach, alongside linearized semiclassical (LSC) dynamics ?,? and symmetrical quasiclassical (SQC) dynamics, ?,? has been established as a robust method for NAMD with mapping Hamiltonian ?,? in condensed-phase systems.? The dynamics obtained from the MSH models align well with the all-atom NAMD results, reaffirming the MSH model’s capacity to accurately capture the essential nonadiabatic dynamics of complex molecular systems in condensed phases, consistent with findings in other systems. ?,? The Marcus rate constants listed in Table and visualized in Figure provide a kinetic basis for interpreting the NAMD results in Figure. For example, comparing R1 [S1] to R0 reveals that CT1 population growth accelerates while CT2 population growth decelerates. This shift can be attributed to the rate constant k 2→3 for the CT1 →CT2 transition, which is 2 orders of magnitude lower in R1 [S1] than in R0 due to a reduced driving force and increased reorganization energy. Thus, the sequential pathway (i) in eq is hindered. Simultaneously, the rate k 1→3 for the EX →CT2 transition increases in R1 [S1] relative to R0, enhancing the direct pathway (ii). In contrast, the population transfer from the initial EX state to the CT states in conformation R2 [S1] is markedly slower than in the other conformations under Scheme 1. This sluggish dynamics arises because all three primary rate constants are low, a consequence of the weak electronic couplings inherent to the bilateral A-D-A arrangement of conformation R2.

The second question we address is determining the dynamical signatures associated with different bath correlation scenarios. We compare the photoinduced CT dynamics over the first 10 ps for MSH models constructed using Schemes 1–3 for conformations R1 and R2 against the reference R0 conformation, as plotted in Figure.

We begin by analyzing the correlated bath case in conformation R1. To isolate the pure bath correlation effect, we compare the dynamics of the reference conformation R0 with that of R1 [S3]. By construction, these two models share an identical system Hamiltonian (electronic couplings and site energies ε _ j ), differing only in their bath parameters. The geometric metric θ_23 decreases from 77.7° in R0 to 71.0° in R1 [S3], indicating a more strongly correlated bath in the latter. This enhanced bath correlation is intrinsically linked to the increase in the reorganization energy E _ r _ ^(12)^, which rises from 0.823 eV in R0 to 1.432 eV in R1 [S3]. Physically, this increase stems from the greater donor–acceptor A_1_ separation in R1 geometry, which requires more extensive solvent rearrangement upon charge transfer. The dynamical consequences of this shift are evident in Figure(c1): the population decay of the EX state and the corresponding rise of the CT1 population are significantly faster in R1 [S3] than in R0, while the population gain in CT2 is slower. This acceleration of the primary EX →CT1 step is driven by a two-order-of-magnitude increase in the rate constant k 1→2. The enhancement in rate is a direct consequence of the activation energy E act ^(12)^ dropping from 0.17 eV in R0 to a nearly barrierless 0.0034 eV in R1 [S3] (using eq), a change ultimately caused by the increased reorganization energy pushing the reaction closer to the activationless regime as shown by cyan square and orange circle in Figure(c1).

Moreover, comparing Schemes 2 and 3 allows us to investigate the effects of bath parameters and thermodynamic driving forces. For the transitions of interest, EX →CT1 and EX →CT2, the angle θ_23_ decreases from 74.3° in R1 [S2] to 71.0° in R1 [S3], indicating a stronger bath correlation in the latter. Figure(b1,c1) shows that the population decay of the EX state and the growth of the CT1 state proceed more rapidly in R1 [S3] than in R1 [S2]. This acceleration correlates with an increased rate constant k 1→2 in R1 [S3], driven by the fact that the thermodynamic driving force (ΔE 12) almost cancels with the larger reorganization energy (E _ r _ ^(12)^), which leads to an activation energy E act ^(12)^ of 0.0034 eV in R1 [S3], significantly smaller than 0.213 eV in R1 [S2]. As illustrated by the orange circles in Figure(b1,c1), these factors push the reaction to nearly barrierless top of the Marcus turnover in R1 [S3] compared to R1 [S2].

Alternatively, comparing R0 with R1 [S2] isolates the dynamical differences arising solely from the change in solute geometry, as R1 [S2] retains the Hamiltonian parameters of R0. The bath correlation angle θ_23_ shifts from 77.7° in R0 to 74.3° in R1 [S2], suggesting a more highly correlated bath in the modified R1 geometry. Both the short-time dynamics in Figure(a2,b2) and the long-time behavior in Figure(b1) reveal faster EX decay and CT1 growth in R1 [S2] relative to R0. This behavior is supported by the kinetic data in Table and Figure(a1–3,b1–3), which show elevated rate constants k 1→2 and k 1→3, alongside a significantly reduced k 2→3 in R1 [S2]. The enhancement of k 1→2 and k 1→3 results from solvent rearrangement following photoexcitation, which induces larger reorganization energies across the three state pairs, coupled with a larger thermodynamic driving force in the inverted regime in R1 [S2] than R0. Conversely, the suppression of k 2→3 stems primarily from a reduced thermodynamic driving force ΔE 23 of 0.047 eV and an amplified reorganization energy E _ r _ ^(23)^ of 2.222 eV, leading to a high activation energy of 0.579 eV in R1 [S2] in the normal regime, effectively disrupting the CT1 →CT2 step of pathway (i) defined in eq.

Beyond bath parameters and solute geometry, we assess the impact of the system Hamiltonian, specifically electronic couplings and excitation energies, by comparing R1 [S1] and R1 [S3]. Because the TCNE A_1_ molecule is displaced further from the MPe donor in R1, the electronic coupling |Γ_12_| is lower in R1 [S1] compared to the optimized R0 structure used in R1 [S3]. However, |Γ_13_| and |Γ_23_| are higher in R1 [S1] according to Table. Figure(a1,c1) demonstrates that in R1 [S1], the EX decay and CT1 growth are slower, whereas CT2 growth is faster compared to R1 [S3]. This trend aligns with the kinetic pathway analysis; although the backward reaction CT2 → CT1 is faster than the forward reaction CT1 → CT2, both are orders of magnitude smaller than the suppressed k 12 and enhanced k 13 rates, as shown in Table and Figure(a1–3,c1–3).

Next, we examine the anticorrelated bath case represented by conformation R2. To isolate the pure bath correlation effect in the anticorrelated regime, we compare the reference R0 with R2 [S3]. As with the previous comparison, the electronic Hamiltonian remains fixed, so all dynamical differences arise solely from the bath. The correlation metric θ_23_ shifts from 77.7° in R0 to 93.9° in R2 [S3], marking a clear transition to an anticorrelated bath. This regime change is energetically manifested by a significant increase in the reorganization energy E _ r _ ^(23)^, which rises from 1.934 eV in R0 to 2.444 eV in R2 [S3]. The molecular origin of this increase lies in the geometry of R2: with acceptor A_2_ positioned on the opposite side of the donor relative to A_1_, transitioning between the two CT states requires a massive solvent rearrangement to invert the polarization field across the donor plane. Dynamically, Figure(c2) indicates that while the EX decay and CT1 population growth are faster in R2 [S3] than in R0, the accumulation of the CT2 population is similar. This behavior is primarily governed by the suppression of the sequential transfer channel; the rate constant k 2→3 drops by 2 orders of magnitude in R2 [S3] compared to R0. This kinetic blockade is directly caused by the elevated activation energy (rising from 0.228 to 0.346 eV), which is itself a consequence of the large reorganization energy E _ r _ ^(23)^ characteristic of the anticorrelated bath. The corresponding Marcus parabolas illustrating this barrier increase are shown in Figure(c3), comparing cyan squares for R0 and purple triangles for R2 [S3].

To investigate the effects of bath parameters and thermodynamic driving forces, we compare R2 [S3] with R2 [S2], as shown in Figure(b2,c2). The degree of bath correlation θ_23_ increases from 93.9° in R2 [S3] to 96.9° in R2 [S2], indicating stronger anticorrelated bath effects in the latter. Figure(b2,c2) reveals that while the EX decay rate remains comparable in both schemes, CT1 growth is faster and CT2 growth is slower in R2 [S3] than in R2 [S2]. This behavior can be rationalized by the kinetic data: the rate constant k 1→2 is larger, while k 1→3 is smaller in R2 [S3] compared to R2 [S2], as depicted by the purple triangles in Figure(b1,b2,c1,c2). It is noted that k 2→3 is higher in R2 [S3] than in R2 [S2], a consequence of the more evenly distributed charge distribution in R2 [S3] leading to a smaller reorganization energy in the normal regime (Figure(b3,c3)).

Alternatively, comparing R0 with R2 [S2] isolates the dynamical impact of the solute geometry change, as both use the same R0 electronic parameters. The bath correlation angle shifts from θ_23_ = 77.7° in R0 (correlated bath) to θ_23_ = 96.9° in R2 [S2] (anticorrelated bath). Both the short-time dynamics in Figure(a2,c2) and the long-time behavior in Figure(b2) show that EX decay and the growth of both CT1 and CT2 are faster in R2 [S2] than in R0. This acceleration is supported by the fact that the sum of k 1→2 and k 1→3 is higher in R2 [S2], despite k 1→2 and k 2→3 being lower and k 1→3 being higher individually according to Table and Figure(a1–3,b1–3). The dramatic four-order-of-magnitude reduction in k 2→3 for R2 [S2] relative to R0 is primarily attributed to a massive increase in the reorganization energy E _ r _ ^(23)^. This increase reflects the significant solvent rearrangement required when the solute geometry changes from having two acceptors on the same side (R0) to having them on opposite sides of the donor (R2), calling for the accommodation of two widely separated acceptor molecules with flipped charge distributions.

Finally, comparing R2 [S1] with R2 [S3] elucidates the effects of the system Hamiltonian, specifically the electronic couplings and excitation energies. While the thermodynamic driving forces remain similar, the electronic couplings are significantly smaller in R2 [S1], which features the bilateral A-D-A arrangement, compared to R2 [S3], which retains the couplings computed for the R0 geometry (Table). Consequently, Figure(a2,c2) shows that EX decay and the growth of both CT1 and CT2 are slower in R2 [S1] than in R2 [S3]. This slowdown is expected, as all dominant reaction rates in the two kinetic pathways defined in eq are reduced, as detailed in Table and Figure(a1–3,c1–3).

Beyond analyzing population dynamics, we can extract atomistic details from the nonadiabatic simulations of the three trimer conformations. Figure plots the changes in the radial distribution function (RDF) for various solute–solvent components at finite time points (0.5, 1, and 4 ps) relative to time zero (photoexcitation). ?,? These time slices were derived directly from the all-atom nonadiabatic simulation of conformation R1, conducted using SQC dynamics with a triangle window function. ?,? In the first row of Figure, we observe a general tendency for the entire solute trimer to slightly repel ACN solvent molecules within a 3–5 Å range from the molecular surface. Specifically, the MPe donor repels ACN molecules, whereas the two TCNE acceptors attract the polar solvent, inducing a shell-like rearrangement at distances of approximately 2–7 Å. The second and third rows of Figure decompose these interactions by displaying the solute moieties relative to the negative N atom and positive H atoms of ACN, respectively. It is evident that as the MPe donor becomes positively charged during the NAMD, it attracts the ACN nitrogen atoms and repels the hydrogen atoms (Figure(b2,c2)). Conversely, the two TCNE acceptors (A_1_, A_2_), which gradually acquire negative charge, exhibit the opposite behavior (Figure(b3,b4,c3,c4)). We propose a semiquantitative relationship for the RDF changes, positing that the intensity of the TCNE acceptors’ interaction with the solvent atoms correlates directly with the product of the CT1 or CT2 population at a given time and the SASA of the corresponding acceptor. For instance, in conformation R1 at 4 ps, the populations of both CT1 and CT2 are approximately 0.1, and the SASA values for both TCNE A_1_ and A_2_ are roughly 0.7. Consequently, we anticipate the dominant peak height ratio A_1_:A_2_ to be near 1:1, a prediction confirmed by the data in Figure(b3,b4) and (c3,c4). Additional RDF data for the equilibrated ground state (time 0) of all conformations, as well as the RDF evolution for other conformations, are provided in the Supporting Information.

As discussed above, nonadiabatic dynamics of photoinduced CT processes in polar solvents are influenced by multiple factors, including electronic couplings, excitation energies, bath parameters, and system-bath couplings. However, the primary objective of this work is to establish a molecular perspective of bath correlation, regardless of electronic couplings and vertical excitation energies. In particular, we extract the fundamental distinctions between correlated and anticorrelated baths from our findings. First, the degree of bath correlation is quantified by the angle θ_23_ between the two transitions of interest, |1⟩→|2⟩ and |1⟩→|3⟩, which is directly determined by the three pairwise reorganization energies according to eq. Physically, the reorganization energy reflects the fluctuations of the energy gap between corresponding states, which are governed by both molecular conformation and charge distribution. In the present case, polar solvent molecules respond to changes in the solute’s conformation and charge distribution upon electronic transition, making the solvent’s rearrangement a major contributor to the reorganization energy. On the molecular level, molecular conformation and charge distribution determine the dipole moment change, and molecular conformation determines the SASA of the solute, so the dipole change and SASA serve as a molecular indicator for the reorganization energy.

The physical picture of a correlated bath for the |1⟩→|2⟩ and |1⟩→|3⟩ reactions is that in-phase bath fluctuations drive the concerted progression of both reactions. This regime is characterized by a relatively small E _ r _ ^(23)^, indicating low energy-gap fluctuations between states |2⟩ and |3⟩. Within the MSH model or the equivalent multistate reaction coordinate (MRC) model, this relationship corresponds to a short equilibrium distance between these states along the collective reaction coordinate. Thus, in a correlated bath, a change in molecular conformation that increases E _ r _ ^(12)^ is expected to simultaneously increase E _ r _ ^(13)^. Generalizing this principle, we anticipate that correlated baths will be expected in cascade-type architectures (e.g., D–A_1_–A_2_ arrays) where sequential charge transfer steps involve charge displacements in similar spatial directions; here, the solvent polarization induced by the first step naturally preorganizes the environment for the second.

In contrast, the picture of an anticorrelated bath describes a scenario where a bath fluctuation favoring one reaction, such as |1⟩→|2⟩, disfavors the alternative reaction, |1⟩→|3⟩. This behavior resembles out-of-phase fluctuations and is characterized by a relatively large E _ r _ ^(23)^. It can be understood as a situation where the two reactions require solvent motions in opposing directions, reflected by a large separation between the equilibrium positions of states |2⟩ and |3⟩ along the collective reaction coordinate. Accordingly, such anticorrelated baths are likely to be common in branched or highly symmetric geometries (such as the sandwich conformation R2 studied here) where alternative pathways involve spatially opposing charge transfer patterns. In these scenarios, the solvent shell cannot simultaneously stabilize both distinct charge-transfer states, leading to the competitive frustration characteristic of the anticorrelated regime.

These theoretical distinctions of bath regimes may translate into specific, measurable experimental signatures, such as kinetic branching ratios and spectroscopic lineshapes. Kinetically, our simulations predict that an anticorrelated bath creates a dynamical bottleneck by blocking the sequential pathway; experimentally, this would be observable as a distinct alteration in the kinetic branching ratio, manifesting as the transient accumulation of the intermediate species and a suppressed quantum yield of the final product state, in contrast to the efficient sequential transfer typical of correlated-bath cascade geometries. Spectroscopically, two-dimensional electronic spectroscopy (2DES) offers a direct probe of these correlations,? as the geometric angle θ_23_ maps to the cross-correlation of energy gap fluctuations. Consequently, we predict that the correlated bath regime will be signaled by cross-peaks with a substantial positive center line slope, indicating synchronous stabilization of the EX and CT manifolds, whereas the anticorrelated regime will exhibit cross-peaks with a negative center line slope, reflecting the opposing energetic demands of the competing pathways.

In conclusion, we establish a uniquely defined metric for bath correlation between two nonadiabatic transitions sharing an initial state, |1⟩→|2⟩ and |1⟩→|3⟩. This metric is the angle θ_23_, which is fully determined by the three reorganization energies E _ r _ ^(12)^, E _ r _ ^(13)^, and E _ r _ ^(23)^ via eq. An angle θ_23_ less than 90° indicates a correlated bath, whereas an angle θ_23_ greater than 90° indicates an anticorrelated bath. These reorganization energies are intrinsically linked to the molecular properties of the electronic states and are strongly affected by changes in dipole moment and SASA upon electronic transition. Thus, the specific molecular conformation and charge distribution determine whether a molecule or aggregate in a polar solvent falls within the correlated or anticorrelated bath regime. Within a correlated bath, energy-gap fluctuations and solvent motions facilitate both reactions, associated with a small E _ r _ ^(23)^. Conversely, within an anticorrelated bath, solvent motions favor one reaction at the expense of the other, associated with a large E _ r _ ^(23)^. This study offers a clear physical picture of the correlated bath effect for nonadiabatic dynamics in polar solvents, while the investigation of bath effects on more complex dynamics in disordered heterogeneous environments ?−? ? ? will be reported in future work.

Supplementary Material