Electron transfer and spin–orbit coupling: Can nuclear motion lead to spin selective rates?

Quantum-mechanical models of electron transfer have proven extremely successful in describing key features of electronic transitions in systems such as biomolecules and solar cells. One of the simplest models that has been effective in such applications is the spin-boson model, where a two-state system (representing two electronic states) is coupled to a thermal bath (represented by a collection of harmonic oscillator modes).¹ In the most basic formulation, one assumes that the two states are coupled together via a constant coupling (V), an assumption known as the Condon approximation. If the set of modes in the thermal bath is indexed by α, the Hamiltonian of such a system is of the form

Effectively, the spin-boson model corresponds to two parabolic potential surfaces shifted spatially and in energy relative to one another, resulting in a single minimal-energy crossing point. Within the context of chemical physics, the spin-boson problem has proved itself useful as a model because (i) the Hamiltonian is simple enough such that one can solve for the dynamics at many different levels of theory and (ii) there are so many time scales of interest (temperature, driving force, diabatic coupling, reorganization energy, and the nuclear relaxation time) that one can explore many different dynamical regimes.² For instance, one dynamical regime is the case where the coupling V is large, termed the adiabatic limit. In this limit, the system evolves as if on a single electronic potential surface; one can invoke standard classical transition state theory (TST) to approximate the rate of reaction and use the Born–Oppenheimer approximation to obtain a reasonable picture of nuclear dynamics.³ A second regime is the case where V is small (the nonadiabatic limit) and the electronic dynamics become slow enough that the system dynamics cannot be described using the Born–Oppenheimer approximation. The usual approach is to apply Fermi’s golden rule and make the high temperature approximation, recovering what is often called Marcus theory.⁴ Yet, another regime, closely related to the nonadiabatic regime, is the solvent-controlled regime where nuclear motion can slow down or damp electronic transitions by trapping transitions with strong friction. Semiclassically, one often differentiates the above regimes through the framework of Landau–Zener transitions and many condensed phase problems occur in the damped and overdamped regimes; to that end, several authors have explored the interplay between electronic transitions and diffusive nuclear motion.^5–7 

Extensions of the spin-boson model Hamiltonian have been proposed and analyzed as well. For instance, several researchers have explored how Marcus theory is altered in the presence of position-dependent diabatic couplings (breaking the so-called “Condon” approximation). Perhaps the most famous example of such a model is due to Medvedev and Stuchebrukhov,⁸ who studied systems where all nuclear vibrational modes could be disentangled into two groups: modes that displace the diabat and modes that modulate the diabatic coupling. In such a regime, a simple rate expression can be derived. These expressions were effectively generalized by Jang and Newton,⁹ who calculated the consequences of non-Condon effects for a variety of different circumstances.

Below, our focus will be slightly different from the papers mentioned above insofar as our aim will be to model the effect of spin–orbit coupling on nonadiabatic electron transfer, a subject that has recently been reviewed by Lykhin et al.¹⁰ and Marian et al.¹¹ In order to treat systems with spins as rigorously as possible, we will not force the diabatic coupling to be real-valued. Instead, we will explore the implications of the fact that the $L̂⋅Ŝ$ spin–orbit coupling operator^12,13 is complex-valued. The result is that a spin-conserving process between two electronic states (say, with spin up) evolves according to a spin-dependent Hamiltonian $Ĥ$⁠, while the dynamics of the corresponding process (say, with spin down) evolve according to a spin-dependent Hamiltonian $Ĥ*$⁠. Thus, if we find differences between $Ĥ$ and $Ĥ*$ dynamics, we will effectively find differences between spin up and spin down coupled nuclear-electron transfer dynamics. A brief justification of this model is provided in Subsection 1 of the  Appendix.

Although all results below will focus exclusively on quantum dynamics (and not address semiclassical dynamics), it is helpful to analyze the differences in H vs H* dynamics in a semiclassical context. When nuclei are propagated classically along an eigensurface corresponding to a complex-valued Hamiltonian, a novel so-called Berry force^14–16 emerges from the fact that no consistent phase can be chosen for the electronic adiabats. This force is of the form

and is equivalent to the force on a classical charged particle moving through a magnetic field. Here, $Ṙ$ is the nuclear velocity, ∇_R represents a gradient taken with respect to nuclear coordinates, and D_nn is the derivative coupling, or Berry connection, from surface n to itself. From Eq. (1.5), it follows that, even though the eigenvalues of H and H* are identical, in a semiclassical context, the coupled nuclear-electronic dynamics are different because the Berry forces on each adiabatic state are equal and opposite (and could potentially lead to different, spin-dependent electron transfer dynamics in the absence of a magnetic field),

For a more nuanced interpretation of Berry force, see Subsection 5 of the  Appendix.

Now, formally, the Berry force in Eq. (1.5) is meaningful only in the limit of adiabatic nuclear motion; in such a case, one can make a simple analogy to dynamics in a magnetic field. Interestingly, recent dynamical scattering simulations have demonstrated that strong Berry-like forces can also appear in the nonadiabatic limit for systems with a few degrees of freedom and accessible conical intersections. In particular, for gas phase dynamics around a conical intersection, a Berry force can dictate which outgoing channel is populated.¹⁷ To our knowledge, however, the implications of Berry force have not been analyzed in the context of Marcus theory and nonadiabatic electron transfer; there is also a fairly small literature on how Landau–Zener and WKB theory¹⁸ changes in the presence of complex-valued diabatic couplings in multiple dimensions.^19,20 More generally, there have been few if any calculations questioning: can the differences between $Ĥ$ and $Ĥ*$ (and the corresponding Berry force effects) lead to different, spin-dependent electron transfer dynamics in the presence of nuclear friction? The present article will take a first step in this direction by using the spin-boson model.

Turning back to the quantum-mechanical description of the nonadiabatic problem at hand, we will make the standard approximation that the diabatic coupling is small so that one can use first-order perturbation theory and employ the Fermi Golden Rule (FGR) ansatz to calculate equilibrium rates,

In Eq. (1.7), $V̂$ can be a function of the nuclear bath coordinates $xα$⁠. This expression gives the rate of population decay from an initial vibronic state $1,v$ to a set of final states ${2,v′}$⁠. Here, 1 and 2 index electronic states and v and v′ index vibrational states. Now, if we assume that (i) the unperturbed Hamiltonian [$Ĥ0$ in Eq. (1.3)] is real and (ii) the initial state $1,v$ is a stationary state of the unperturbed Hamiltonian $Ĥ0$⁠, it follows then that $⟨2,v′|V̂|1,v⟩*=⟨2,v′|V̂*|1,v⟩$⁠, i.e., $⟨2,v′|V̂|1,v⟩2=⟨2,v′|V̂*|1,v⟩2$⁠. In other words, to first order in $V̂$ or $V̂*$⁠, the rate of electronic relaxation for $Ĥ$ must equal the rate of relaxation for $Ĥ*$⁠. Put bluntly, if the dynamics are initiated from quasi-equilibrated starting conditions, no spin-dependent (Berry force) effect can arise to first order in perturbation theory.

With this background in mind, it becomes clear that if such spin effects are to arise, they must enter either from from higher orders in perturbation theory or from nonequilibrium starting conditions. With regard to the former, we note that higher order effects can and do have important dynamical consequences. For instance, the famous Zusman result^5–7 is a higher order effect beyond FGR, demonstrating that, with enough friction, electron transfer rates can be dominated by solvent relaxation even for small diabatic couplings. Thus, in the future, we expect that an interesting research direction will be probing if/how Berry force effects emerge at second or higher order in the perturbation.

For the present article, our focus will be on the latter possibility, i.e., the possibility that $Ĥ$ vs $Ĥ*$ (Berry force) effects may arise from nonequilibrium initial starting conditions. Such conditions are relevant to many electron transfer processes, which occur before equilibration or a steady-state can be reached. This scenario is true particularly in the study of photo-induced electron transfer, where the timescale of light-induced electronic excitation is much faster than that of nuclear vibrational relaxation.²¹ Intuitively, we might expect that the phases of wavepackets take on greater significance for a complex-valued Hamiltonian, and thus, coherences may play a role in amplifying Berry force effects. In order to understand the short-time behavior of relaxation in such systems, equilibrium methods are not sufficient. Fortunately, nonequilibrium systems can be physically modeled through a nonequilibrium formulation of Fermi’s Golden Rule,²² a brief review of which will be provided in Sec. II.

At this point, all that remains is to introduce the spin-boson-like complex-valued Hamiltonian that we will study to explore electron transfer with spin. After investigating several different possibilities, we have settled on the following exponential spin–orbit coupling (ESOC) [with $Ĥ0$ as in Eq. (1.3)]:

A rough physical picture for how such a complex-valued, position-dependent coupling arises can be found in Sec. V. The interstate coupling in Eq. (1.8) remains bounded for all displacement of the system (i.e., $V̂$ does not diverge even when $xα→∞$⁠), making this ESOC model particularly amenable to theoretical investigation. More precisely, we will be able to apply FGR and NEFGR (non-equilibrium FGR) techniques to this Hamiltonian and gain a basic insight into how Berry force does or does not manifest in molecular electron transfer rates. Note that a change in the sign of the W_α parameters corresponds to a change in the spin and this aspect will be discussed often. Again, see Subsection 1 of the  Appendix.

Finally, as a motivation for the present research, it is worthwhile to put the present results in the context of the larger phenomenon of chiral induced spin selectivity (CISS).^23,24 Over the last 15 years, a host of experiments have shown that the electronic current running through a chiral molecule is often quite spin polarized.²⁵ At the moment, there is no uniformly agreed upon explanation for this phenomenon, given the very small magnitude of the spin–orbit coupling and Zeeman effects within simple organic molecules.^26–28 Over the last two years, several research groups^29–35 have explored if and how CISS might arise from nuclear-electronic interactions, i.e., the entanglement between electronic transitions, spin-dependent Berry forces, and nuclear motion. Indeed, the goal of Ref. 17 was to show that Berry forces can lead to strongly spin-dependent nuclear motion even for systems with small spin–orbit coupling. In this work, we will further show that such H vs H* dynamical (Berry force) effects can also survive friction for some period of time (which is consistent with recent calculations of spin polarization emerging at a molecular junction under bias³⁶ in the presence of electronic friction as caused by the production of electron–hole pairs). Although not indicative of any causation, this finding of robust spin effects is consistent with observations of the CISS effect at ambient temperatures for a wide range of molecular and solid state environments.

An outline of this paper is as follows. In Sec. II, we briefly review standard nonequilibrium Fermi’s golden rule theory; in Sec. III, we formally solve for the spin-boson dynamics [with Eq. (1.8)]. In Sec. IV, we present our numerical results (both transient and long time). In Sec. V, we interpret and discuss our findings as far as spin-dependent electron transfer dynamics. For all calculations that follow, we use units such that ℏ = 1.

Often, one is interested in the scenario whereby a vibronic system is initialized on electronic state $1$ and one interrogates the probability of reaching electronic state $2$ at time t. Let $ρ̂(0)$ be the initial density matrix, and we imagine partitioning the Hamiltonian as $Ĥ=Ĥ0+V̂$⁠, where $Ĥ0$ is the unperturbed Hamiltonian and $V̂$ is the perturbation; see Eq. (1.1).

According to standard practice, we transform to the interaction picture, perturbatively time-evolve the density matrix to first-order, and finally trace over the vibrational bath states {v′} in the $2$ electronic state. The resulting expression, which provides the two-state population as a function of time, can be conveniently expressed using a time autocorrelation function C(t′, t″) of the interaction potential operator $V̂I(t)=eiH0tVe−iH0t$⁠,

where Tr_B denotes a trace over the bath states.

At this juncture, according to a standard FGR calculation, one makes the approximation that the bath vibrational modes begin in thermal equilibrium. Thus, if the vibrational states for state 1 are labeled ${v}$⁠, the initial density matrix is as follows:

Note that $ĤB$ acts only on the vibrational components of any eigenstate, and ∑_v denotes a sum over all vibrational eigenstates of the bath. If one plugs Eq. (2.2) into Eq. (2.1), one finds that the correlation function reduces to a function dependent on only the difference t′ − t″, and the final result in the energy domain is Eq. (1.7). Note that evaluating either Eq. (1.7) or Eq. (2.1) can be difficult, although the task is made much easier when we assume that the bath is harmonic.

Now, in order to generalize the FGR approach above to nonequilibrium initial conditions in a tractable manner, the usual prescription^22,37 is to start with the equilibrium state above and then shift the eigenstates of each mode in space by a certain distance away from equilibrium. Mathematically, we can capture these initial conditions by applying to each bath eigenstate the following unitary translation operator:

This operation produces the nonequilibrium density matrix we will use as the initial state for subsequent calculations,

Equation (2.4) represents a mixed state where each eigenstate has been shifted in space along each mode α by a distance d_α. From here, following Refs. 37 and 22, one can plug Eq. (2.4) into Eq. (2.1) and use the same manipulations as in a typical FGR calculation to obtain a nonequilibrium expression. Within this class of nonequilibrium initial conditions, one can model photophysical and photochemical experiments, for which it is reasonable to assume that the bath was originally thermalized in the ground electronic state before an electronic transition is driven. Thereafter, system dynamics are launched along an excited electronic state. See Fig. 1.

Finally, let us address the question of convergence. For perturbative results (based on equilibrium or nonequilibrium initial conditions) to be valid, it is generally sufficient that the dynamics induced by the perturbation correspond to the slowest time scale of the system being analyzed. For spin-boson and spin-boson-like problems, this condition requires that the nuclear dynamics not be significantly affected by the slow population leakage from one electronic state to the other, a condition which is dictated by the size of the perturbation and the frequencies of the bath. The validity of the NEFGR is not directly dictated by the size of the shift d in Eq. (2.4).³⁸ 

Having reviewed the FGR/NEFGR formalisms, we will now apply the methods to the Hamiltonian with the coupling in Eq. (1.8).

We begin with equilibrium. We will now demonstrate that the rate dynamics show no dependence on the sign of the terms $Wα$ if the initial state of the bath is a mixed state of vibrational eigenstates on either diabat. In other words, if the initial bath subspace density matrix is diagonal, there will be no Berry force effects in the population dynamics. Note that this statement is very strong and applies to many initial conditions; after all, thermal equilibrium is just one special case of a diagonal bath density matrix.

To derive our result, we begin by applying a polaron transformation to the ESOC Hamiltonian, which is the typical approach in working with the standard spin-boson problem. The unitary operator that carries out the transformation is as follows:

In this new so-called polaron representation, the diagonal electronic energies are decoupled from the boson bath; all coupling to the boson bath is captured purely in the interstate coupling. The exact form of the Hamiltonian is as follows:

where $V$ is a constant that differs from V by a complex phase factor. The terms λ_α are the distances between the minima of states 1 and 2 along the alpha mode coordinates. Furthermore, we note that under the polaron transformation, eigenstates of either diabat are transformed to eigenstates of $ĤB$⁠. In other words, the transformed initial density matrix is

In Eqs. (3.2) and (3.3), we have added a superscript tilde over the energy $Ẽj$ and the density matrix $ρ̃̂$ and used stylized $H$ and $V$ so as to signify that these are quantities calculated after the initial polaron transformation.

Next, we introduce the unitary operator

which propagates each individual mode forward in time by an amount

We note that, as far as each mode (x_α, p_α) is concerned, $R̂$ is a function of the harmonic oscillator Hamiltonian, $p̂α22m+12mωα2x̂α2$⁠. Thus, $R̂$ preserves the eigenstates of $ĤB$⁠, merely multiplying them by a phase factor. Transforming $Ĥ$ under this operator leads to the (final) transformed Hamiltonian

At the same time, the initial density matrix is unchanged through this transformation,

Further details of this transformation (conjugation by $R̂$⁠) are given in Subsection 2 of the  Appendix.

In the end, the signs of the W_α terms completely disappear when we transform both (i) the Hamiltonian and (ii) the equilibrium initial state. Taking Eqs. (3.5) and (3.6) together, we find that (to any order) all vibronic dynamics depend only on the magnitudes of the terms W_α, not on their signs: in fact, the form of $Ĥ$ suggests that the system has a dressed reorganization energy that can be expressed as a sum of dressed single-mode reorganization energies,

The dressed reorganization energy $Er(tot)$ is composed of two distinct components: The first is the standard reorganization energy, which we will henceforth refer to as $Er(s)$⁠,

The second is a term arising from geometric phase effects, which we will term $Er(g)$⁠,

On the one hand, if one were to spectroscopically characterize the reorganization energy of a vibronic system with spin–orbit effects by measuring fluctuations in the energy gap, one would measure $Er(s)$⁠. On the other hand, if one were to conduct a series of experiments for a particular electron transfer system to construct a plot of ΔG vs rate, one would observe a shift in the peak corresponding to $Er(tot)$⁠. In other words, $Er(g)$ arises only when there is nuclear motion of interest. With this physical interpretation in mind, we will term $Er(tot)$ the total reorganization energy, $Er(s)$ the static reorganization energy, and $Er(g)$ the geometric reorganization energy. See Fig. 2.

Three key features of the Hamiltonian were necessary for the above proof: The first is that the bath modes were purely harmonic. The second is that the modes maintained their respective frequencies when there was a change in the electronic state. Third and finally, we relied on the fact that the initial state was a mixed state of vibrational eigenstates of one of the two system diabats [which allowed us to write Eq. (3.6)]. As mentioned in the Introduction, this state of affairs suggests that more distinct electron transfer rate effects may well emerge with either anharmonic bath models and/or nonequilibrium starting conditions.

In order to treat the nonequilibrium case using the formalism from above, the only change we must make is to modify the initial density matrix. We will use the density matrix defined in Eq. (2.4). As shown in Subsection 3 of the  Appendix, the resulting autocorrelation function is

Substituting this expression into the integral in Eq. (2.1), we obtain a general first-order population expression,

Here, P₂(t) describes the gradual development of population in state $2$ starting with P₂(0) = 0. One can take the time derivative of the equation above and obtain a time-dependent rate of transition from the $1$ state to the $2$ state. The formal (time-dependent) rate constant is given in Eq. (A24). In the high temperature Markovian limit, the rate expression can be explicitly evaluated (see Subsection 4 of the  Appendix), producing a time-dependent golden rule rate expression of the form

where $Er(tot)$ is the total reorganization energy that arises due to the combined effects of $Er(s)$ and $Er(g)$⁠. Note that Eq. (3.14) for the time-dependent rate constant has a form that is very similar to that of a classical Marcus theory time-independent rate expression. In this high temperature limit, the dressed reorganization energy $Er(tot)$ and dephasing function s(t) are enough to fully characterize the nonequilibrium rate dynamics. These parameters are, in turn, characterized by the bath density of modes, as well as the shift parameters λ_α, d_α, and W_α. The physical meaning of each of these terms is summarized in Table I and Fig. 3 as a reminder to the reader.

Let us now focus on the dephasing function s(t) in Eq. (3.15). This function contains all information regarding the transient behavior of the system, serving to capture the fluctuations in the interstate energy gap as the nuclear geometry evolves. Note, in particular, that the second term in s(t) changes when any single W_α changes sign, clearly indicating a spin-dependent rate effect. Additionally, in the case that the number of modes in the bath is very large, s(t) tends to 0 as t grows large, implying that the system tends to equilibrium dynamics. Therefore, s(t), indeed, contains information about the dephasing of the bath. As the bath dephases and returns to equilibrium, we expect, based on the result from Sec. III A, that all Berry force effects should be damped out.

Before concluding this section and presenting our results, let us note that, in order to simulate photoexcitation in a two-state system, one can simply set d_α = −λ_α for all α (see Fig. 3), and such a choice will manifest itself directly in the dynamics of the dephasing function in Eq. (3.15). Note, however, that the present model can also be applied to systems where the final acceptor state is not the same as the ground state, such as multistate models and/or intersystem crossings. In such cases, to retain full generality, one must not assume d_α = −λ_α.

In a typical analysis of a spin-boson system (i.e., in the Condon approximation), one assumes that the system–bath interaction can be characterized by a spectral density. Working from such an approximation, one can compute the reorganization energy, which fully characterizes relaxation dynamics. However, in the model we have investigated thus far, a single spectral density is insufficient to characterize the system–bath interaction; the system exchanges energy through both the system–bath coupling and spin–orbit effects so that we must fix two sets of parameters, λ_α and W_α.

To make such a parameterization as simple as possible, we will assume a density of modes $ρ(ω)=(ηω/ωC2)exp(−ω/ωC)$⁠. ω_C describes a cutoff frequency for the bath and η is a normalization factor that ensures that integration over the density function returns the correct number of modes in the bath. For the diabats, we will then further assume that the mass-weighted shifts between the two surfaces are uniform across all modes, i.e., $λαmα=λ$ for all α. Similarly, we set $Wα/mα$ equal to a constant W for all α. Note that W can be positive or negative. For the bath displacement parameters d_α, we set the mass-weighted displacements $dαmα$ equal to −κλ. As discussed earlier, κ = 1 corresponds to relaxation of a photoexcited system.

In order to characterize the dynamics, we need to compute the reorganization energy and the dephasing function. We begin with the reorganization energy. The total reorganization energy is the sum over all modes of the single-mode reorganization energies,

Below, it will sometimes be helpful to switch variables from (λ, W) to $(Er(tot),θ)$⁠,

Next, we calculate the dephasing function s(t),

In terms of $Er(tot)$ and θ, the dephasing function (plotted in Fig. 4) is

In Fig. 5, we show rate constants and population yields that were calculated for a photogenerated (κ = 1) distribution using the high temperature expression, Eq. (3.14). We fix the static reorganization energy $Er(s)$ and add in some amount of geometric reorganization energy $(Er(g)=W22m)$⁠; the ratio is quantified by θ in Eq. (4.2).

In Fig. 5(a), we plot the population for short and long times for different values of θ. At long times, as expected, we see that the dynamics for ±θ gradually approach the same equilibrium rate as time progresses. Here, we emphasize that the disagreement between the full NEFGR expression and the high temperature expression for θ = ±0.5 is not the result of any interesting dynamical effect that depends on θ. Instead, a simple numerical calculation shows that the full perturbational and high temperature expressions differ simply because of tunneling effects (whose importance depends on the total value of the reorganization energy $Er(tot)$⁠). Note that $Er(tot)(θ=0.5)=Er(tot)(θ=−0.5)≠Er(tot)(θ=0)$⁠.

The more interesting behavior occurs at short times-while the bath still has yet to decohere and where the population curves can depend sensitively on ±θ (i.e., ±W). To see these dynamics clearly, we zoom in within panels (b) and (c) of Fig. 5. Here, we find transient populations for θ = −0.5 that are 50% larger than those of θ = 0.5. In general, Berry force effects are larger when we calculate populations with the full (rather than the high temperature) expression, but overall, the high temperature expression clearly reflects the correct qualitative physics (especially if we compare only θ with −θ data).

Finally, in panels (d) and (e) of Fig. 5, we plot the derivative of the population data, i.e., the instantaneous rates. Here, we can clearly see differences in the features of the early time dynamics, including the magnitude and the relative timing of the populations spikes. Again, we find that the high temperature approximation is a reasonably good approximation for the full perturbational early time dynamics (and computationally cheaper to evaluate).

Finally, and most importantly, we would like to ask whether any rate differences between systems with opposite signs (±W) can persist for a significant period of time. In order to probe this question for a wide variety of θ values (while avoiding excessive computational expense), we will use the high temperature expression above. We can roughly quantify dynamical differences by focusing on rate dynamics at the time 1/ω_C, i.e., the characteristic dephasing time of the bath. For this dataset, we will fix the total reorganization energy $Er(tot)$⁠; if we were to fix $Er(s)$ as before, the variation in $Er(tot)$ would be very large and make it difficult to isolate effects arising from the geometric component. (Moreover, as noted above in Fig. 5, the quantitative accuracy of the high temperature approximation depends on the total value of $Er(tot)$ and whether one can safely ignore tunneling; by fixing $Er(tot)$⁠, we can thus reasonably assume a cancellation of error and compare the effect of exchanging $Er(g)$ with $Er(s)$⁠.)

Our results are plotted in Fig. 6. We consider two related quantities. First, we plot the spin polarization at time 1/ω_C,

where k(t) is as defined in Eq. (3.14) and k₋ is related to k₊ through a change in sign on the {W_α} terms. As shown in Fig. 6, the polarization can become quite large, reaching nearly 60% in photoexcitation processes in the 200–400 K range. Note the antisymmetry in the polarization plot across the line θ = 0. This antisymmetry arises naturally because changing the sign of θ corresponds to changing the signs on the {W_α} terms (which results in opposite spin polarization).

Second, in Fig. 6, we also plot

in order to quantify the absolute magnitude of the effect of the $Wα$ terms. We find that the nonequilibrium rates with W ≠ 0 can be very different from the rates with W = 0, sometimes by as much as a factor of 1000 for systems prepared far from equilibrium. Interestingly, the rate at time 1/ω_C is never significantly faster than the equilibrium rate; but there are certain regimes in which the transfer rate is much slower. This result suggests that the spin polarization observed is a result not of favoring one electronic spin, but rather disfavoring the other. Finally, we observe that as the temperature of the system is raised, polarization effects due to spin–orbit effects gradually disappear, reflecting the fact that geometric phase effects are highly dependent on bath coherences.

Although the high temperature expression is useful, it cannot provide a complete picture of the system behavior at low to intermediate temperatures. In Fig. 7, we plot the transient spin polarization as a function of temperature for θ = −0.5 using the full NEFGR expression (not the high temperature expression). We find a more nuanced picture; the spin polarization does not increase monotonically as one lowers the temperature, but rather has a maximum at some intermediate temperature. Interestingly, such non-monotonic temperature dependence also has been observed (experimentally and theoretically) by Das and co-workers in the context of spin polarization arising from conduction through DNA base pairs.³⁹ 

From the results presented above, it is clear that Berry forces can survive nuclear friction and that spin-dependent electron transfer dynamics are, indeed, possible. Although spin–orbit coupling effects (i.e., the effects of the sign of W) vanish in equilibrium and contribute only a small amount to the total reorganization energy, these effects can be quite important for systems prepared far from equilibrium. At the same time, however, the present work only scratches the surface of how electronic spin ties into electron transfer and opens up a slew of questions.

First, what is the most reasonable Hamiltonian that one should use when describing spin-dependent electron transfer? In the present article, we have simply assumed an exponential complex-valued diabatic coupling that changes according to $ei∑αWαx̂α$⁠. Where does such a term come from? In principle, we expect an interstate coupling in the electronic Hamiltonian to be of the form

where $a(x̂)$ is a standard, real-valued diabatic coupling and $ib(x̂)$ is the purely imaginary-valued spin–orbit coupling, as calculated in Subsection 1 of the  Appendix. Thus, the phase $d(x̂)$ must originate from a competition between different effects (e.g., the ratio between a position-dependent spin–orbit coupling and a position independent diabatic coupling or the ratio between a position-dependent diabatic coupling and a position independent spin–orbit coupling). If one assumes that $c(x̂)$ is a constant and one expands $d(x̂)≈∑Wαx̂α$ as a Taylor series, one recovers the ESOC Hamiltonian in Eq. (1.8). At present, our group is expending a great deal of effort running ab initio calculations searching for systems where |W| will be large and/or meaningful. One must wonder how large can |W| be, in practice (especially if SOC is small), and is it reasonable to make the exponentiation in Eq. (5.1) over the relevant volume of configuration space where two diabatic curves cross?

Questions about the validity of the ESOC model can be addressed in the future by investigating alternative forms for the interstate coupling, for instance, the more common Linear Vibronic Coupling (LVC) type model.^37,40,41 In the LVC model, the interstate coupling is taken to be linear in position for all bath modes, i.e., $V̂=∑αcαx̂α12+cα*x̂α21$⁠. The LVC model is standard for investigating dynamics of systems around conical intersections, but this model must also break down at some point because the off-diagonal couplings grow to infinity away from the origin. Nevertheless, if one works with an LVC (or an LVC-like model), a second question one can ask is what dynamics do we observe when a conical intersection is broken by a small amount of spin–orbit coupling? Note that, within an LVC model, adding spin–orbit coupling will be completely different from adding a real-valued constant coupling (of equal magnitude); a real-valued coupling would shift the LVC conical intersection, whereas an imaginary spin–orbit coupling would remove the conical intersection entirely (and lead to a large W). While such a removal might also not seem entirely realistic,^42–44 this scenario makes clear that electron transfer effects may emerge that are wholly unique to strongly spin-dependent systems.

A third question relates to the displaced harmonic oscillator model for nuclear motion. Above, the equivalence of equilibrium dynamics for $Ĥ$ and $Ĥ*$ arose almost accidentally and is clearly tied to the ESOC Hamiltonian. Although first-order FGR rates may be identical for $Ĥ$ and $Ĥ*$ for all rate processes [see Eq. (1.7)], there is no reason to expect (in general) that these equilibrium rates will be the same at second or higher order. For example, in Fig. 8, we plot dynamics for one single anharmonic mode, where dynamics are initialized in the lowest eigenstate of the upper diabat. Clearly, the dynamics show a small (but nonzero) spin-dependence. In general, will these differences be small (because they show up only at higher orders of perturbation theory) or will these differences be significant (because they can add together for many anharmonic modes)? One must also be very curious about the dynamics of a model with Duschinskii rotations⁴⁵—do such rotations (which introduce stronger directionality on the bath reorganization) amplify the effects of spin-dependent nuclear forces? One significant challenge will be to quantify the magnitude of spin-dependence for equilibrium dynamics in systems beyond the displaced oscillator approximation.

Finally, the goal of this research direction is to identify realistic systems that will display strong Berry force effects. Obviously, one would like to work with ab initio electronic structure Hamiltonians and run dynamics, but such a course of study is expensive. Even if, for a moment, we commit to studying model Hamiltonians, it is not wholly clear how to best map model dynamics to the real world. For instance, even though the standard spin-boson model has been used very successfully to model electron transfer over the years, it is clear that a more complicated parameter space of models will be necessary to discern Berry force effects. After all, within the ansatz of a two level system, Berry force arises by breaking the Condon approximation, and so we must require not one, but two spectral densities; one of these spectral densities mediates energy exchange between the electronic and vibrational systems (as in the standard spin-boson model) and the other characterizes how the vibrational states mediate the spin–orbit coupling between electronic states. Do these interactions proceed with the same natural set of vibrational states? How should these spectra be tuned to one another in order to maximize spin-related electron transfer effects? Can we disentangle these effects in a meaningful way? Finally, if we return to the question of electronic structure, one must ponder how best to map a realistic system to the model Hamiltonians discussed above: in particular, will we necessarily find that, within floppy molecular structures, metal ions with large spin–orbit couplings always experience bigger Berry force effects during an electronic transition or are the dynamics more complicated? For instance, some recent experiments suggest that the CISS effect can be sensitive to an ion with strong spin–orbit coupling,⁴⁶ while other experiments suggest that spin–orbit coupling in a substrate is not important to the overall CISS signal.⁴⁷ If we can broadly answer such questions, we will clearly gain a great deal of intuition as far as rationally designing viable molecules and pathways with exciting new chemistry and ideally gain the capacity to separate metal ions with large spin–orbit coupling from those with small spin–orbit coupling.

In conclusion, we have investigated a complex-valued spin-boson-like model Hamiltonian using nonequilibrium Fermi golden rule theory in order to provide a fundamental insight into how spin–orbit coupling might influence condensed phase electron transfer. Our investigations have shown that Berry force effects can survive nuclear friction for some (measurable) period of time, before the effect dissipates. The exact details of how long such an effect survives are complicated and sensitively depend on the exact nuclear motion—the ESOC model presented here is simple to treat theoretically, but difficult to map to an ab intio Hamiltonian. Further research will be necessary to better quantify the effect and gain intuition as to exactly how big these Berry forces can be (and when) and if/how nuclear vibrations either correlate with or cause a CISS effect.⁴⁸ Moreover, in the future, if we wish to run predictive simulations of realistic molecules with ab initio potentials, it will be necessary to develop efficient and inexpensive semiclassical algorithms^17,49,50 for treating nuclear curve-crossing dynamics in the presence of spin–orbit coupling.

This work was funded by the Center for Sustainable Separation of Metals, an NSF Center for Chemical Innovation (CCI), under Grant No. CHE-1925708.

The authors have no conflicts to disclose.

Data sharing is not applicable to this article as no new data were created or analyzed in this study.

In this appendix, we briefly discuss the connection between electronic spin and complex Hamiltonians and then proceed to derive some of the necessary equations presented above for the nonequilibrium golden rule calculations.

Consider a molecular system with two electronic states labeled $1$ and $2$ at a fixed geometry. In adding a spin–orbit effect, we must consider spin states $1,↑$⁠, $2,↑$⁠, $1,↓$⁠, and $2,↓$⁠. In principle, we can write the Hamiltonian as a sum

where $Ĥ0$ is spin-independent and ξ is the strength of the spin–orbit interaction. Since $Ĥ0$ is spin-independent, it can only couple states of identical spin. Therefore, in the spin-state basis, we can write H₀ as follows:

Now, we construct the $L̂⋅Ŝ$ matrix. We note that if $1$ and $2$ are molecular orbitals in an appropriate representation, they do not couple to themselves through any of the three component angular momentum operators. Thus, we can write the $L̂⋅Ŝ$ matrix as follows:

Noting that the angular momentum component operators are purely imaginary, we can write the matrix more concisely,

Note that α is purely imaginary, whereas β has a real component, so the upper right block is anti-Hermitian and traceless. Therefore, this 2 × 2 matrix is diagonalizable and has purely imaginary eigenvalues α′ and α′* = −α′. The lower left block is just the conjugate transpose of the (anti-Hermitian) upper right block, and so is diagonalized by the same change of basis. Furthermore, the eigenvalues of the lower left block are conjugate to the eigenvalues of the upper right block. Note that a change of spin basis does not affect $Ĥ0$⁠, which is the block diagonal. We can now write down the total fixed-geometry Hamiltonian in the new spin basis [recall that ξ is defined in Eq. (A1)],

Therefore, at any nuclear geometry, the four-dimensional Hilbert space above (for spin-electronic wavefunctions with two spatial orbitals) can always be partitioned into two distinct subspaces, each with a distinct Hamiltonian $Ĥ$ and $Ĥ*$ corresponding to one spin or the other opposite spin.

In this article, we have made the assumption that that such a partition will be valid and unchanged for all nuclear geometries. This strong assumption allows us to assume that we can entirely ignore all spin-flip processes.

To evaluate Eq. (3.5), we begin with the polaronic Hamiltonian in (3.2), restated here as

We will use $R̂$ as defined in Eq. (3.4) to perform a unitary transformation on the polaronic Hamiltonian. We need to compute $R̂ĤR̂†$⁠. $R̂$ obviously commutes with $ĤB$⁠. Thus, all we need to compute is $R̂ei∑αWαx̂α+λαp̂αR̂†$⁠. Using the fact that (for any operator $M̂$⁠)