Infinite-order perturbative treatment for quantum evolution with exchange

See allHide authors and affiliations

Science Advances  07 Aug 2020:
Vol. 6, no. 32, eabb6874
DOI: 10.1126/sciadv.abb6874


Many important applications in biochemistry, materials science, and catalysis sit squarely at the interface between quantum and statistical mechanics: Coherent evolution is interrupted by discrete events, such as binding of a substrate or isomerization. Theoretical models for such dynamics usually truncate the incorporation of these events to the linear response limit, thus requiring small step sizes. Here, we completely reassess the foundations of chemical exchange models and redesign a master equation treatment for exchange accurate to infinite order in perturbation theory. The net result is an astonishingly simple correction to the traditional picture, which vastly improves convergence with no increased computational cost. We demonstrate that this approach accurately and efficiently extracts physical parameters from complex experimental data, such as coherent hyperpolarization dynamics in magnetic resonance, and is applicable to a wide range of other systems.


Calculations of quantum evolution in dynamic systems, such as exchange or conversion between multiple discrete states, are important today in many disciplines (15). These calculations first became prominent in magnetic resonance more than 50 years ago with the McConnell equations (6), which were introduced first as a classical approximation of the spin dynamics in exchanging systems. These equations could readily describe the dynamic spectra of uncoupled spin-1/2 nuclei but were incapable of handling evolution under bilinear couplings. In contrast, the density matrix formalism (79) readily includes statistical averaging in the equilibrium state, and coherent evolution can be handled by unitary transformations involving calculation of a highly accessible propagator for spin systems.

As it would be computationally impossible to explicitly calculate, for instance, the dynamics of 1020 nuclear spins, one averages over each molecule to form a reduced density matrix, wherein the form of ensemble interactions is obfuscated. Therefore, dynamic exchange effects require a more careful treatment of the expression of the ensemble action in the reduced density picture, modifying the time evolution from the form given by pure quantum mechanics (tρ̂=i[ρ̂,Ĥ]). The exchange interaction has been historically derived as an analog to the case of Redfield relaxation theory (7, 10), but the ensemble dynamics that generate relaxation occur on a time scale far faster than the evolution of the quantum degrees of freedom (femtoseconds to picoseconds), effectively limiting the influence of these dynamics to the first observable moment. This is not valid for exchange, where it would be feasible for higher moments of the ensemble interaction to act on a time scale comparable to the coherent evolution.

Despite the maturity of models for exchange, there is still considerable motivation to develop new methods to efficiently and accurately explore dynamic effects in systems undergoing quantum evolution. On the forefront of magnetic resonance techniques are hyperpolarization methods (1115), which overcome the intrinsically low signal-to-noise limits by distilling spin order from an external source. Of particular interest over the past decade is signal amplification by reversible exchange (SABRE) (5, 1531), in which the singlet order of parahydrogen is converted into observable magnetization or more complex spin states on target ligands during transient interactions with an iridium catalyst (Fig. 1A). Optimization of this technique requires accurate modeling of the system, which, with the recent advent of coherently pumped SABRE experiments (Fig. 1B), has revealed bizarre and complex dynamics (5, 32). Accurately modeling the coherent hyperpolarization dynamics of systems like (15N,13C)-acetonitrile (Fig. 1C) and subsequently fitting the experimental data have been impossible within previous frameworks for exchange, given the multitudinous exchange interactions, such as coligand exchange events and number of coupled spins (21 total spins).

Fig. 1 SABRE provides an ideal system to challenge the limits of an exchange model, given the complexity of the underlying dynamics.

(A) SABRE transfers the singlet order of parahydrogen to a target ligand via reversible interactions with an iridium catalyst and exhibits nonlinear dynamics that are highly dependent on the relative concentrations of each species. (B) The coherent hyperpolarization dynamics can be probed by interleaving pulses at or near the SABRE resonance condition (red) with periods far off resonance (B = − 22.5 μT) to allow for exchange. (C) The (15N,13C)-acetonitrile SABRE system demonstrates rich dynamical information that varies with the resonant field as a result of a complex coupling network between 15N, 13C, and 1H. Lines are meant as visual guides only.

To that end, we completely reinterrogate the incorporation of dynamic exchange interactions in evolving quantum systems. We construct a reformulated dissipative master equation by recovering the traditional expression from the Dyson series and then continuing the derivation to infinite order in perturbation. The ramifications of extending the derivation of the dissipative master equation to all orders in the exchange interaction make a profound impact on the radius of convergence of exchange simulations with absolutely no additional computational cost by deriving a simple scaling factor that accounts for all moments in the ensemble motion. In addition to the most general case of exchange between distinguishable ensembles, we show solutions for pseudorotation generated by Abelian groups of order 2 and 3 as well as for quantum dynamical selection (QDS), where coherent degrees of freedom alter the exchange interaction. By coupling this new infinite-order treatment exchange to a formulation of the exchange operators that scales linearly with the number of distinguishable ensembles, we can easily model highly complex systems that would be untenable within alternate exchange formalisms (3, 25).


Pursuing a traditional, master equation approach to chemical exchange requires the assumption that fundamentally discrete exchange events may be approximated as a continuous perturbation to the ensemble, shifting the model from assuming a Poisson process of a microcanonical ensemble to a Wiener process of a canonical ensemble. Quantum Monte Carlo (QMC) models discontinuously sample exchange events and are essentially exact, provided that one can iterate the solution to convergence (5). However, the cost of iterating a QMC solution to convergence would often make the calculation intractable on the time scale of experimental guidance. For instance, on the canonical 14-spin bis-(15N-pyridine) SABRE system, it would take approximately 300 years to run a single 60-s simulation at a modest exchange rate (kex = 50 s−1) to approximately 99% convergence.

The methods to describe dynamic evolution in quantum systems are well established in the case of spin relaxation in liquids, which relies on the perturbative Dyson series expansion of the interaction-frame propagator. In this case, as well as for the traditional case for exchange, the series is truncated to the leading term by assuming that the dynamic interaction is a small perturbation. The same result can be recovered by annealing the Liouville–von Neumann equation to the rate equations defining exchange by taking the tensor product between the quantum and chemical degrees of freedom. In both cases, this assumes that exchange acts linearly on the evolving quantum degrees of freedom, which is not well motivated. With our ansatz that discrete exchange events may be approximated as a continuous process, we recast the Dyson series for the case of exchange without any a priori assumptions as to the magnitude or order of the exchange interaction.

Reformulation of the Dyson series for chemical exchange

To begin, we partition the Hamiltonian into a stationary component Ĥ0 and a stochastically modulated exchange interaction Ĥ1(t), for which the equation of motion is given by (ℏ = 1)tρ̂=i[Ĥ0+Ĥ1(t),ρ̂](1)

Boosting Eq. 1 into the interaction frame givestσ̂=i[H˜̂1(t),σ̂](2)where we use the convention ρ̂σ̂ for distinguishability between the representations and H˜̂1(t)exp (iĤ0t)Ĥ1(t)exp (iĤ0t) is the interaction representation denoted by the tilde. Equation 2 may then be formally integrated and iteratively substituted back into the expression to generate the Dyson series. When doing so, we assume that the correlation time of the exchange interaction is much faster than the coherent evolution of the system, which allows us to write σ̂(t)=σ̂(t)t and extend the upper limit of integration to infinitytσˆ(t)=i[H˜ˆ1(t),σˆ(t)]T0dt[H˜ˆ1(t),[H˜ˆ1(t),σˆ(t)]]+(3)

In this equation, T is the Dyson time-ordering operator, which imposes t > t. At this juncture, we introduce H˜̂1(t) as the operator expansionH˜̂1(t)qF̂q(t)A˜̂q(4)where q indexes through uncoupled exchange mechanisms, F̂q(t) are real, stochastic operators describing the time evolution of exchange, and A˜̂q define the interaction of exchange with the evolving quantum system. We assume that the system is in a chemical steady state, and therefore, F̂q(t)=0 and is importantly not the first moment of the exchange rate. Furthermore, this has the repercussion that all odd-order terms in the expansion necessarily average to zero, ensuring that chemical exchange generates no complex phase rotations in σ̂. Substituting Eq. 4 into Eq. 3 and ensemble averaging gives the leading termtσˆ=Tpq0dtAˆˆpeiHˆ0(tt)AˆˆqeiHˆ0(tt)σˆFˆp(t)Fˆq(t)(5)

Note that we have shifted to the commutation superoperator representation ([Ô,]Ô̂) and have dropped the formal time dependence of σ̂ as well as the tilde notation for the interaction frame for legibility. The integrand is then the correlation function of the ensemble motion, which, for a stochastic, real-valued, and time-continuous (Wiener) process, is δ-correlated in time. To avoid violating the time-ordering operator, T, only self-correlated terms can give nonzero amplitudes upon integration, which may be accounted for by imposing δpq. The δ-correlation imposes that exp (iĤ0(tt)) is unity, making all interaction-frame superoperators for exchange identical to their lab space representation. Then, the time-ordered integral simply determines the rate of the exchange process, which is the probability of exchange given a characteristic lifetime (τq) during a finite period of time, (Δtτq1)/Δt. Upon integration, this givesT0dtFˆp(t)Fˆq(t)=0dtδ(tt)δpqΔtτq1Δt=1τq(6)

The term Â̂pÂ̂q can now be written as Â̂q2, which is immediately realized in its more common notation as the exchange superoperator K̂q, where the sign convention arises to ensure that exchange forms a completely positive map. Together with Eq. 6, this returns the canonical form of the chemical exchange master equation, which, in the Schrödinger representation, istρ̂=i[Ĥ0,ρ̂]+qK̂qρ̂τq(7)

Truncating the expansion here essentially imposes that exchange only acts linearly on ρ̂ during the finite period of time Δt (i.e., the simulation step size), as this was the condition under which Eq. 6 was defined. In only the simplest cases can Eq. 7 be homogenized and analytically integrated, wherein the finite period of time Δtdt and the assumption that an exchange interaction of any magnitude acts linearly on the quantum degrees of freedom are accurate. However, now that we have explicitly established the derivation of the chemical exchange master equation from the Dyson series, it is simple to continue the derivation to higher-order terms. Remembering that all odd-order terms necessarily go to zero upon ensemble averaging, the next nonzero interaction in the expansion arrives in the fourth-order term, which, after substitution of Eq. 4 and realizing that there can be only powers of Â̂q2 to give rise to the superoperator K̂q, istσˆ=pqAˆˆp2Aˆˆq2σˆT0dt0tdt0tdtFˆp(t)Fˆp(t)Fˆq(t)Fˆq(t)(8)

The integrand of Eq. 8 is a four-point correlator that may be factored into a sum of two-point correlator products by Isserlis’ theorem, where the form of each correlator is given by Eq. 6. There are (n − 1) !! identical terms for an n-point correlator after factorization when the process is δ-correlated, where n !! is the semifactorial of n, defined as the factorial using only integers of the same parity as n (5 !! = 5 × 3 × 1). In addition, given the time symmetry of the Wiener process, there will be (n − 1)! degenerate time orderings upon integration, accounted for with division of the correlator amplitude by the degeneracy. Integration then givesT0dt0tdt0tdtFˆp(t)Fˆp(t)Fˆq(t)Fˆq(t)=(41)!!(41)!(Δtτq)2/Δt=Δt2τq2(9)

Notice that δpq prevents any cross terms between p and q from arising in the summation, hence why integration generates a rate proportional to τq2. Equation 8 then becomestσ̂=q(Â̂q2)2Δt2τq2σ̂(10)

The fourth-order term describes the probability of two exchange events occurring during a finite period of time, as the exchange interaction can be expanded into successive applications of Â̂q2. We shall define the powers of the exchange interaction conditioned to specific cases and otherwise leave it in its more general form. Using the assumptions established here, it is then beneficial to rewrite the entire Dyson series for exchange astσ̂=q{1τqk=0(Â̂q2)k+11k!(Δt2τq)k}σ̂(11)

Equation 11 may be simplified by establishing a more rigorous definition of the exchange operator Â̂q2, which will show for a general case and more specific applications. Before doing so, it is pertinent to note that the general operator action of Eq. 11 can be written as(Aˆˆq2)k+1σˆ=(Aˆˆq2)kKˆqσˆ=γKˆqσˆ(12)where the first equality is inherent given the definition of K̂q and the second is possible if K̂qσ̂ are eigenfunctions of (Â̂q2)k, where γ is then a constant. Under this condition, the infinite sum in Eq. 11 would be independent of the interaction superoperator and evaluation of all moments of the exchange interaction would have no additional computational cost over a traditional formulation.

Exchange between distinguishable ensembles

The most general formulation for exchange is to form a composite vector space constructed from the direct sum of m discrete chemical configurations that form manifolds of quantum states, where exchange allows flow of populations between manifolds via projection operations. This formulation is unrestricted in the systems that it may describe, as it is trivial to project between systems of different sizes and projections may both encapsulate linear and nonlinear contributions to the evolution. When constructed carefully, this method grows linearly in cost with m, as projections need not act on the entire composite vector space.

We find the form of the exchange interaction for this case by recognizing that projections are idempotent operations (γ = 1 in Eq. 12), such that(Aˆˆq2)k+1=(1)k+1(Aˆˆq2)k+1=(1)k+1(Kˆq)k+1=(1)kKˆq(13)

Note that we have reindexed the alternating term for convenience. There is an intricate ramification of Eq. 13, in that the group Gq containing all powers of the exchange superoperators for a given system is isomorphic to its coset of linear power superoperators Sq(1), which equivalently form the kinetic equations for the chemical dynamics. Hence, considering all moments in the exchange interaction will only ever generate dynamics that are directly reflected in K̂q. Given Eq. 13 and that the summation over k is simply the Maclaurin series for the exponential, Eq. 11 may be written in the Hilbert space astρ̂=i[Ĥ0,ρ̂]+q{K̂qτq exp (Δt2τq)}ρ̂(14)

In this form, Eq. 14 is the exact, closed form solution of quantum evolution with exchange, which we call the exact dissipative master equation (DMEx). In its exact form, the only difference that arises in comparison to the traditional form is the exponential factor, where it is clear that in the limit where Δtdt, the equation converges back to Eq. 7. That arises as the impact of higher moments in exchange go away over small periods of time because it is impossible for multiple exchange events to occur simultaneously in the limit of infinitesimal step sizes. However, as Δt becomes larger, the higher-order terms account for moments in the dynamics of the ensemble that are not present when one assumes a linear coupling between the quantum and exchange degrees of freedom.

Pseudorotation by Abelian permutation groups

While Eq. 14 is valid for any exchanging quantum system, it would be inconvenient to expand a system into separate manifolds when the coherent evolution within those manifolds, or a subset of those manifolds, is identical. This is the case for exchange generated by pseudorotation, such as the canonical example of cyclohexane inversion in magnetic resonance. Hence, it is convenient to recast Eq. 11 for the cases of pseudorotation generated by two- and threefold symmetric permutation groups, which we shall call G2 and G3 pseudorotations, respectively. As we are recasting the DMEx for specific cases, we will drop the q-index and write the explicit form of the equation of motion. Note that all of the assumptions made to derive Eq. 11 remain valid for these conditions.

In the case of pseudorotations, which contain inherently coupled exchange processes, we define the first moment of the exchange interaction operator asÂ̂2K˜̂=12(R̂R̂1+R̂1R̂)Ê=12(K̂+K̂)(15)

The operators R̂±1R̂1 generate the forward (K̂) and backward (K̂) rotations, which are coupled with equivalent probability. We have also written the exchange superoperator K̂ in a traceless form for convenience instead of writing a separate term proportional to unity (Ê), as done in the general case. The forward and backward rotations are equivalent for G2 pseudorotations, which allows Eq. 15 to be reduced toK˜ˆ=RˆRˆ1Eˆ(16)

The form of the higher powers of the exchange interaction are given as(Aˆˆ2)k+1=(2)kK˜ˆ(17)and hence give the DMEx for G2 pseudorotation in the Hilbert spacetρ̂=i[Ĥ0,ρ̂]+R̂ρ̂R̂1ρ̂τ exp (Δtτ)(18)

Again, the exponential term returns the canonical equation of motion for exchange in magnetic resonance when taken to the limit of analytical integration (Δtdt). Note that in this case, the argument of the exponent is proportional to 1/τ, as opposed to the other cases presented here. This arises as a ramification of the definition that was used for K˜̂.

The case of G3 pseudorotations no longer permits the reduction of Eqs. 15 and 16, as KˆKˆ. It is pertinent to note the relation Rˆ2Rˆ2=Rˆ1Rˆ, which allows for any quadratic or higher rotation to be written in terms of the linear term. The higher powers of the exchange interaction then reduce to the simple form(Aˆˆ2)k+1=(1)kK˜ˆ(19)

Plugging Eq. 19 into Eq. 11 gives the DMEx for G3 pseudorotationtρ̂=i[Ĥ0,ρ̂]+K˜̂ρ̂τ exp (Δt2τ)(20)whereK˜̂ρ̂=12(R̂ρR̂1+R̂1ρR̂)ρ̂(21)

Again, the only difference between the exact treatment and the traditional implementation for exchange is the exponential term. Note that Eqs. 18 and 20 can be used alone or in conjunction with the more general DMEx formalism shown in Eq. 14 as a method of contracting the composite vector space, with an obvious example being to use Eq. 20 to compress the manifolds corresponding to the three orientations of methyl rotation with magnetically inequivalent interactions. This does require that the rate connecting the contracted manifolds to any other manifold be identical; otherwise, the contraction is invalid.

Quantum dynamical selection

An interesting case to examine is when the quantum and chemical dynamics are coupled, such as in QDS, where quantum evolution dictates the evolution of the exchange degrees of freedom. In this case, we must reassess the assumptions made to arrive at Eq. 11 pertaining to the stochastic motion of the ensemble. However, we will construct this case as exchange between distinguishable ensembles, leaving Eq. 13 intact. This case is restricted by the condition that[F̂q(t),σ̂]0(22)such that we may no longer ensemble average σ̂ and the ensemble motion operators F̂q separately. This gives rise to the leading term in the Dyson seriestσˆ=pqAˆˆpAˆˆqT0dtFˆp(t)Fˆq(t)σˆ(23)

The factorization of the three point correlator generatesF̂p(t)F̂q(t)σ̂=F̂p(t)F̂q(t)σ̂+F̂p(t)σ̂F̂q(t)+F̂q(t)σ̂F̂p(t)(24)

As F̂q=0 is retained with the assumption of a chemical steady state, only the first term in the factorization is retained. Furthermore, as the ensemble is still macroscopically described by a Wiener process, the δ correlation is retained and δpq is imposed to avoid violating the time ordering. For any (2n + 1)–point correlator, only the single term that averages σ̂ separately will be retained. If we then insert Eq. 24 into Eq. 23 and defineT0dtFˆq(t)Fˆq(t)Φˆˆq(25)where Φ̂̂q is an ensemble motion superoperator that determines the rate of the process as a function of the population in the quantum state coupled to the exchange process, then we obtain as the equation of motiontσˆ=qAˆˆq2Φˆˆqσˆ(26)

We may impose thatΦˆˆqσˆ=ζqσˆ(27)such that ζq is a constant that is the rate of the ensemble motion dictated by Â̂q, which states that σ̂ is an eigenstate of Φ̂̂q at all times. This relation is validated when δ is imposed, which prevents the Φ̂̂q from generating a tensor of rank 1 or higher. Hence, we may write the form of Eq. 25 asΦˆˆq=T0dtδ(tt)Δt/(Δtτq1)Pq=Pqτq(28)where Pq scales the rate by the projection of the quantum state, which couples to the exchange process and is a constant. As projection operators are idempotent, we may derive the DMEx for QDS directly from the derivation of the general case, with an additional factor of Pq that makes the rates time dependent; in the Hilbert space, this istρˆ=i[Hˆ0,ρˆ]+q{KˆqPqτq exp (ΔtPq2τq)}ρˆ(29)

Expressing the exchange rates as time-dependent quantities complicates interpretation of this case. However, if we choose to express the ensembles using both chemical and quantum degrees of freedom, then we recover Eq. 14 as the equation of motion for exchange. The exchange rates between the redefined manifolds are then time independent.


Performance of DMEx models

The chemical exchange master equation is only homogeneous and analytically integrable in the simplest of cases and can acquire nonlinearities when one considers reversible exchange between distinguishable ensembles. The DMEx equations of motion presented here are fully converged in the exchange interaction and handle that interaction exactly but are not exact for the full evolution, as there is no preservation of the time ordering between the quantum and exchange degrees of freedom that would be present in a full solution. This could be approached by calculation of higher-order time derivatives but at an additional computational cost. However, the next section focuses on the difference between the traditional master equation approach (Eq. 7) and the DMEx solutions, which shows marked improvements in permissible time steps, and then compares the DMEx to QMC simulations where these time orderings are preserved, but the computational overhead is very large. Therefore, we will evaluate the performance of DMEx methods using the solutionρˆ(t+Δt)=Uˆρˆ(t)Uˆ+qΔtτqexp (Δt2τq)(KˆqEˆ)(Uˆρˆ(t)Uˆ)(30)

In this equation, Uˆexp (iHˆ0Δt). This is an ideal computational method, as it only involves forward propagation of the solution, requires the fewest number of matrix operations, and produces linear evolution under the spin Hamiltonian and evolution to all orders in the exchange interaction. Equation 30 has a small intrinsic error associated with the solution, in that the first step only evolves quantum mechanically. A more accurate way to solve the equation of motion would be to evolve the initial density matrix backward in time by Δt/2 and then using Eq. 30 to generate the solution. Doing so shifts the actions of Uˆ and K̂q by a half step and corrects for this initial error. However, we have found that this makes little difference in the solution; thus, we retain Eq. 30 so to avoid generating a nonintegral number of steps. To isolate errors arising from exchange, we have constructed all of the following simulations in the Hilbert space of the system, where one can exactly evaluate quantum evolution in systems up to 15 coupled spins.

Dynamic nuclear magnetic resonance spectra under pseudorotation have been studied and understood for decades. Spectral features are well resolved in the limit of slow exchange, which broaden and coalesce as the exchange rate increases, and ultimately result in line narrowing in the fast exchange limit. This is reflected in spectra of s-trioxane (9) undergoing ring inversion (Fig. 2A), where the axial (blue) and equatorial (red) have different chemical shifts and the geminal 2JHH coupling is observable. As exchange increases, the spectrum collapses to a singlet, as the axial and equatorial positions become, on average, magnetically equivalent. Similar effects appear for the tert-butyl rotation in t-BuPCl2 (Fig. 2B), which additionally exhibits a transition that is invariant under exchange and thus does not broaden (33).

Fig. 2 Comparison of the traditional master equation solution (Eq. 7, black) and our DMEx models for G2 and G3 pseudorotations.

We use s-trioxane ring inversion to model G2 pseudorotation (A) and tert-butyl rotation in t-BuPCl2 to model G3 pseudorotation (B) The graphs at the bottom compare the root mean square deviation (RMSD) of the generated magnetization as a function of the time evolution such as in Fig. 1C using Δt = 10 μs (which is taken as the ground truth). The DMEx models retain good fidelity with no additional computational overhead, even with step sizes commonly being 10 times larger than were possible with traditional solutions.

For either of these systems, the pseudorotation matrices are generated by expressing a spin label permutation matrix in the appropriate basis. For convenience, we will use the Zeeman basis in this example. In the case of s-trioxane, where the axial and equatorial protons interchange, it is most convenient to setup the system such that axial protons have odd indices and equatorial protons have even indices. Then, the rotation R̂ is given byRˆ=Pˆ12Pˆ34Pˆ56(31)

P̂ij is the permutation matrix that interchanges the ∣αβ〉 and ∣βα〉 states, whereas the ∣αα〉 and ∣ββ〉 states are invariant under exchange. Using this method, it is trivial to arbitrarily reindex the entire system and is computationally efficient because the transformation from the original basis to the reindexed basis is unitary.

When calculating these spectra, we find that the traditional implementation and the DMEx converge to the same solution as Δtdt. However, the DMEx exhibits a substantially smaller error at any step size than the traditional implementation and only accrues an error on the order of 1% when the step size exceeds the average lifetime. In this limit, the traditional implementation loses stability and the trace of ρ̂ deviates from unity. This immediately provides the ability to take larger step sizes with the DMEx implementation. In the case of s-trioxane, an error in the solution of ≈1% requires Δt = 1 ms in the DMEx and Δt = 0.1 ms using the traditional implementation, thus requiring one to sample far fewer data points. In considering all moments of the exchange interaction, the radius of convergence of the Dyson expansion is far larger than it would be by assuming conditions similar to those used for exchange.

While these model systems provide illustrative examples of the performance of the DMEx model, they are far from the more challenging cases in dynamic systems. As noted previously, an interesting system that has gained much attention in the past decade is the hyperpolarization method SABRE, wherein large nonthermal nuclear magnetization is distilled from parahydrogen via reversible interactions with an iridium catalyst. Current efforts are focused on optimizing the extraction of spin order from parahydrogen, which requires accurate modeling of the quantum and exchange dynamics in realistic systems. For reference, an example simulation of the coupled coherent and exchange dynamics that drive SABRE hyperpolarization is shown in Fig. 3A, where the evolution of the 15N polarization is calculated under the experimental conditions for SABRE–SHield Enables Alignment Transfer to Heteronuclei (SHEATH) (16).

Fig. 3 Evaluation of model errors in 15N-SABRE simulations using our SABRE-specific DMExFR2 model, which adds free ligand effects, rebinding, and relaxation (FR2) to the DMEx.

An example of SABRE hyperpolarization dynamics is shown for reference, calculated on a six-spin 15N SABRE-SHEATH system (A). Comparing the DMEx and a QMC treatment, which is viewed as the “gold standard” but is computationally inefficient, indicates that there is a genuine but small difference of 0.142 ± 0.018%, on average, between the two solutions [red data, (B)]. The convergence error of the QMC is indicated by the black line. This error in the DMEx is attributed to the loss of the time orderings between the quantum and exchange degrees of freedom that are retained in the QMC. Even with nonlinear effects incorporated in the simulation, the DMExFR2 exhibits a larger radius of convergence over the traditional implementation by approximately a factor of 4 (C) and an improved self-consistency (parameterized by σk, the error in the predicted exchange rate) (D), which uses Δt = 10 μs solutions as ground truth.

In deriving the DMEx, we began with the ansatz that exchange could be considered as a time-continuous perturbation of the ensemble, but it is interesting to see when this assumption fails. The perturbation generated by exchange in the slow exchange limit is small, allowing the solution to be largely dictated by the quantum dynamics, and conversely, in the fast exchange limit, quantum evolution cannot generate large excursions from equilibrium when constantly disrupted by exchange. In the intermediate regime where SABRE exists, characterized by exchange rates on the order of the dominant couplings, it is no longer trivial to motivate that large excursions from equilibrium would not be impactful on the dynamics. To probe this, we compared our previous QMC model for SABRE (5) against the DMEx on a three-spin SABRE system (Fig. 3B) with a dominant coupling of 2JNH = −24 Hz. In this regime, there is a significant difference between the convergence error of the QMC solution (σQMC) and the DMEx solution; however, this error is, on average, only 0.142 ± 0.018%. Note that this analysis is limited to the smallest systems, given the large cost of iterating the QMC solution, and the error accrued by the DMEx is negligible on simulation time scale relevant to experimental guidance.

When modeling more complex systems, such as those often found in SABRE, it is critical for the cost of the DMEx to be augmented with an efficient method for exploring complex interactions, otherwise circumventing the benefits of an infinite-order treatment by excessive computational costs. In SABRE, these interactions include quantum evolution of multiple species, rebinding of previously polarized ligands to the activated complex, binding site competition with spin-inert coligands, and relaxation. We call this SABRE-specific model the “DMExFR2” to indicate that free ligand, rebinding, and relaxation effects are included. The most efficient way of accomplishing an efficient implementation of exchange is by expressing the interactions as block diagonal with respect to individual manifolds, which we call “manifold-diagonal” for simplicity and will motivate using the example of SABRE.

In SABRE, we primarily consider two different species: one in which the hyperpolarization target is bound to the iridium, which we call “bound species,” and another in solution, which we call “free species.” Coherent evolution in the manifolds is established by separately propagating a bound species density matrix (ρ̂bS) and the dissociated free species density matrix (ρ̂fS) under their respective nuclear spin Hamiltonianstρ̂fS=i[ĤfS,ρ̂fS](32)tρ̂bS=i[ĤbS,ρ̂bS](33)

We begin with the dynamics of the free species. Exchange facilitates association of free species to the catalyst (K̂a,fSρ̂fS), acting on ρ̂fS to remove free species from the manifold as they bind the complex, and allows for dissociation of bound ligand K̂d,fSρ̂bS, adding species back to the manifold. Both of these exchange processes happen at the exchange rate of the ligand, kN, with an action scaled proportional to the ratio of theconcentration of the iridium complex to the free ligand ([Ir]/[S]) to account for the inherent trace normalization of density matrices. The association operator is then simplyK̂a,bSρ̂fS=[Ir][S]ρ̂fS(34)

The dissociation operator then deposits an equivalent number of ligands from the bound species subsystem into the free species subsystem. For the case where both available binding sites in the iridium complex are exchanging with the target ligand, there are distinct subsets of the nuclear spins in the bound species (Sa and Sb), which may dissociate to join the free species with equal probability. We average these possibilities to generate the dissociation operator for the free species, remembering to apply the concentration scaling factor for exchange between manifoldsK̂d,fSρ̂bS=[Ir]2[S](Tr{H2Sa}ρ̂bS+Tr{H2Sb}ρ̂bS)(35)

Combining Eqs. 34 and 35 yields the equation of motion for the free species with exchange interactionstρ̂fS=i[ĤfS,ρ̂fS]+k˜N[Ir][S](K̂d,fSρ̂bSρ̂fS)(36)where kNkNexp (ΔtkN/2) and will be used as a notation for the DMEx rate going forward.

The bound species has two exchange interactions: one for the simultaneous exchange of a ligand and the hydrides occurring at rate of kH and one for the exchange of target ligands at a rate of kNkH. We will formulate the exchange operator for the bound species as a single entity, K̂ex, which takes multiple manifolds as arguments. Hydride exchange is restricted to occur only during ligand exchanges as the complex form a tetrahydride intermediate to facilitate this reaction. In the case where both parahydrogen and ligand exchange occur concurrently, we exchange the portion Δt(ka,H/kN) of ρ̂bS to reflect the new hydride population and new ligand population. This may be written as(k˜a,Hk˜N)ρ̂pH2ρ̂fSTr{H2,S}ρ̂bS(37)where ρ̂pH2 is the density matrix for pure singlet parahydrogen and Tr{H2, S} returns the density matrix for the ligand that remains bound. In the case where the hydrides do not exchange but the target ligand does, another portion of the density matrix Δt((1ka,H)/kN) must be reformulated to reflect the newly exchanged ligand(1k˜a,Hk˜N)Tr{S}ρ̂bSρ̂fS(38)where Tr{S}ρ̂bS is the density matrix for the subsystem of the remaining ligand and parahydrogen. This projection must be constructed carefully to ensure that the coherences are appropriately retained between the hydrides and remaining ligand. Note that while we are exchanging between the free and bound species subsystems, the scaling factor is not needed as the free species density matrix is, by definition, trace normalized. Therefore, one free ligand equivalent leaving the free species will look like one free ligand equivalent associating with the bound species. As this free ligand leaves the free species though, the appropriate reduction in the free species density matrix must be scaled by the concentration ratios. The full exchange operators can now be written as a combination of these two componentsK̂{S}(ρ̂bS,ρ̂fS)=(1k˜a,Hk˜N)Tr{S}ρ̂bSρ̂fS+(k˜a,Hk˜N)ρ̂pH2ρ̂fSTr{H2,S}ρ̂bS(39)

The two possible ligand exchanges from the two available binding sites, a and b, then average together to give the final exchange operator for the bound speciesK̂ex=12(K̂{Sa}+K̂{Sb})(40)

Note that Eqs. 35 and 40 contains terms that are quadratic in the magnetization density, arising from the effects of rebinding ligands that have already interacted with the species. Hence, this is a second-order nonlinear partial differential equation, which must be solved simultaneously with the equation of motion for the free species to define the full evolution of the system.

Furthermore, these nonlinearities are amplified as the ratio increases. It now becomes possible to efficiently represent the impact of concurrent evolution of the J-coupling networks in the free and bound species of the target ligand. In addition, we can now model the effects that various solution compositions will have on the polarization dynamics, given that rebinding of previously polarized ligand will significantly affect the evolution of the bound species under the nuclear spin Hamiltonian.

Even with the incorporation of the nonlinear terms to the DMEx, the solution convergence is still far faster than that of the traditional implementation (Fig. 3C), and the two models still converge in the limit when Δtdt. One can obtain the same error in the DMExFR2 with a step size that is four times larger than the traditional implementation. While we have focused on the accuracy of the simulation, its precision in reproducing input parameters, such as the exchange rate, are just as important, particularly as these models are used to extract physical parameters from experimental data. Under this condition, it is not only critical that the simulation is stable but also efficient, as large portions of phase space have to be searched to perform an experimental fit. To characterize the precision of the simulation, we introduce the parameter σk, which defines the relative shift in the predicted exchange rate in the simulation (Fig. 3D). Unexpectedly, there was essentially no exchange rate at which the traditional implementation provided a solution that was stably precise. In contrast to that, the DMEx model essentially perfectly reproduces the input exchange rate until k ≈ 300 s−1, and when nonlinearities are introduced to the simulation, the maximum deviation from the input exchange rate is only σk ≈ 0.5%.

Practical examples

As noted previously, guided in silico exploration of novel experimental methods that increase the hyperpolarization of SABRE is the focus of optimization efforts in the community. With the improved stability of the DMEx models, it is possible to explore realistic systems with complex coupling networks and reduce the calculation to an obtainable cost by using large simulation step sizes (Δt > 1 ms). The flexibility of this formulation to be expressed in either Hilbert or Liouville space additionally provides access to much larger spin systems than previously possible.

The case of the canonical bis-(15N-pyridine) SABRE-SHEATH system is particularly interesting, as it contains 14 strongly coupled spins in just the “iridium-bound” manifold with 22 total spins and is perhaps the most prevalent system in 15N SABRE. As the full system is far outside the scope of previous exchange models for SABRE, it has been traditionally acceptable to truncate the spin system to an approximate system, fully or partially removing ancillary 1H nuclei, with the largest approximation reported in literature using a single 1H per ligand (4). Even in this case, the dynamics of the truncated model diverge greatly from the actual system dynamics (Fig. 4A), the latter which can be explicitly calculated using the DMExFR2 model with either 2-ms (black) or 5-ms (red) step sizes with only minor deviation between the solutions. As a result, the truncated model optimizes to exchange rates that are false while retaining a deviation of ≈10% from the actual system dynamics when reoptimized to the erroneous rates (Fig. 4B). This means that any physical parameters extracted from experimental data by the model will be greatly confounded by the truncation errors inherent to the formulation.

Fig. 4 Importance of modeling SABRE systems using complete models.

(A) To reduce computational overhead, virtually all calculations to date remove ancillary spins from the system, such as artificially reducing the number of protons on the pyridines in the bis-(15N) SABRE complex from 10 to 2 (A). Doing so alters the hyperpolarization dynamics (blue) as compared to the explicit 14-spin calculation (black), which is stable using the DMExFR2 models for step sizes even up to 5 ms (red). (B) Fitting the 14-spin calculated time evolution with this smaller model produces incorrect values of the exchange rates. (C) Including all relevant exchange pathways when modeling SABRE systems is also crucial for predicting accurate exchange parameters. Here, we fit the experimental (15N,13C)-acetonitrile hyperpolarization dynamics to DMExFR2 models with (solid) and without (dashed) coligand exchange effects. When neglecting these exchange pathways, the predicted exchange rates differ from the correct values by 44 to 92%.

To emphasize the efficiency and flexibility of this framework, we used the DMExFR2 model to fit the coherent hyperpolarization dynamics of (15N-13C)-acetonitrile when exciting the sample with short (milisecond) pulses tuned to a field near the SABRE resonance condition, as described in our prior work (5, 32). Coherent evolution is then interrogated by varying the resonant pulse length, which encodes the dynamics in the final polarization detected. This is a multicomponent SABRE system containing 21 total spins and requires consideration of hyperpolarization-inactive coligand effects to accurately describe the dynamics. These effects allow for additional exchange pathways to influence the dynamics of the system. One of the most critical ramifications arising in allowing the hyperpolarizable ligand to exchange between positions on the complex and thus with which parahydrogen-derived hydride the ligand is coupled. In the limit of fast exchange, this makes the hydrides appear equivalent and would prevent the singlet order from being converted into observable magnetization. When coligand effects are included (solid lines), the experimental data can be reproduced with high fidelity to experiment at multiple field conditions (Fig. 4C), such as when the resonant pulse is B = − 1.65 (red) or −0.91 μT (blue). Furthermore, the extracted exchange rates for these datasets are kN = 14.5 ± 1.8 s−1 and kH = 6.00 ± 0.75 s−1 for the B = − 1.65 μT data and kN = 15.0 ± 3.3 s−1 and kH = 4.50 ± 0.98 s−1 for the B = − 0.91 μT data. When coligand effects are neglected, the predicted exchange rates can have errors of 44 to 92%.

Properly simulating this system requires two seven-spin manifolds for the two conformers of the iridium complex, a five-spin manifold for the free (15N-13C)-acetonitrile, and a two-spin manifold for parahydrogen. Fitting the experimental data would be intractable within any of the previous formalisms for SABRE dynamics (3, 25) as a function of the system size. However, when built using the DMExFR2 model in conjunction with the manifold-diagonal formalism for exchange introduced here, each simulation dataset, which consists of 32 simulations lasting 30 s using Δt = 2 ms, requires approximately 15 min to calculate, making a grid optimization fit possible within a day.


The foundations of exchange in dynamic quantum systems have been reassessed and derived in its exact form from infinite-order perturbation theory. In doing so, the DMEx formalism presented here accounts for higher moments in the ensemble action that are omitted from the traditional implementation. The speed and accuracy with which complex exchanging spin systems may be modeled using the DMEx formalism allows for extensive in silico experimentation and optimization in a way that has previously been inaccessible. In tandem, we have introduced a manifold-diagonal implementation for exchange, allowing the simulations of multicomponent systems with nonlinear exchange dynamics to scale linearly with the dimension of the composite space. Expressing individual manifolds in their Hilbert space representation affords efficient simulation of experimental data with high fidelity.

We anticipate that the results demonstrated here will have a radical impact on the simulation of complex dynamic systems. In the hyperpolarization community, the ability to accurately and efficiently simulate the entire SABRE system should greatly alter the optimization of the hyperpolarization efficiency. Annealing the DMEx formalism to state space reduction techniques has the potential to introduce efficient simulation of systems as large as biomolecules, offering the possibility to markedly reduce computational time and improve simulation accuracy using larger time steps. While these results are presented within the framework of SABRE, they are immediately applicable to the dynamics of any quantum system undergoing exchange, provided that the chemical dynamics are at equilibrium. More broadly, these results demonstrate a method to encapsulate the coupling between distinguishable ensembles past the perturbative limit. This has the potential to have a transformative impact in the field of open quantum system dynamics for applications like quantum computing, where quantum information loss to a bath is critical for understanding the limits of the system.


Coherent SABRE-SHEATH experiments

A sample of (15N,1-13C)-acetonitrile (100 mM), pyridine (33 mM), and IrIMesCODCl (5 mM; IMes, 1,3-bis(2,4,6-trimethylphenyl)-imidazol-2-ylidine; COD, 1,5-cyclooctadiene) was prepared in methanol-d4 and activated by bubbling parahydrogen gas (7 bars) through the solution. This sample was then bubbled inside of a solenoid positioned within a triple-layer micrometal shield, which was stroboscopically pumped between the resonant field condition (B = − 1.65 or −0.91 μT) for a given pulse length interleaved with a nonresonant field of B ≈ − 22.5 μT for 250 ms; this was repeated for 60 s to generate SABRE hyperpolarization. The coherent dynamics are interrogated by varying the resonant pulse length.


Supplementary material for this article is available at

This is an open-access article distributed under the terms of the Creative Commons Attribution license, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.


Acknowledgments: Funding: The presented research was funded by the NSF grant CHE-1665090. Author contributions: J.R.L. derived the form of the infinite-order correction to the chemical exchange master equation, here called DMEx. S.L.E. developed the methods for higher-order effects, such as free ligand evolution and rebinding. Together, J.R.L. and S.L.E. wrote the simulation code and performed error analysis of the various models. C.P.N.T. collected the experimental data. J.R.L., S.L.E., and W.S.W. prepared the manuscript. Competing interests: The authors declare that they have no competing interests. Data and materials availability: All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Additional data related to this paper may be requested from the authors.

Stay Connected to Science Advances

Navigate This Article