## Abstract

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.

## INTRODUCTION

Calculations of quantum evolution in dynamic systems, such as exchange or conversion between multiple discrete states, are important today in many disciplines (*1*–*5*). 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 (*7*–*9*) 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 10^{20} 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 (*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 (*11*–*15*), 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*, *15*–*31*), 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 (^{15}N,^{13}C)-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).

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*).

## RESULTS

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-(^{15}N-pyridine) SABRE system, it would take approximately 300 years to run a single 60-s simulation at a modest exchange rate (*k*_{ex} = 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

Boosting Eq. 1 into the interaction frame gives

In this equation, *t* > *t*^{′}. At this juncture, we introduce *q* indexes through uncoupled exchange mechanisms,

Note that we have shifted to the commutation superoperator representation * _{pq}*. The δ-correlation imposes that

*) during a finite period of time,*

_{q}The term

Truncating the expansion here essentially imposes that exchange only acts linearly on *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 Δ*t* → *dt* 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

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 gives

Notice that δ* _{pq}* prevents any cross terms between

*p*and

*q*from arising in the summation, hence why integration generates a rate proportional to

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

Equation 11 may be simplified by establishing a more rigorous definition of the exchange operator

*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

Note that we have reindexed the alternating term for convenience. There is an intricate ramification of Eq. 13, in that the group G* _{q}* containing all powers of the exchange superoperators for a given system is isomorphic to its coset of linear power superoperators

*k*is simply the Maclaurin series for the exponential, Eq. 11 may be written in the Hilbert space as

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 Δ*t* → *dt*, 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 G^{2} and G^{3} 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

The operators ^{2} pseudorotations, which allows Eq. 15 to be reduced to

The form of the higher powers of the exchange interaction are given as^{2} pseudorotation in the Hilbert space

Again, the exponential term returns the canonical equation of motion for exchange in magnetic resonance when taken to the limit of analytical integration (Δ*t* → *dt*). 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

The case of G^{3} pseudorotations no longer permits the reduction of Eqs. 15 and 16, as

Plugging Eq. 19 into Eq. 11 gives the DMEx for G^{3} pseudorotation

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

The factorization of the three point correlator generates

As * _{pq}* is imposed to avoid violating the time ordering. For any (2

*n*+ 1)–point correlator, only the single term that averages

We may impose that* _{q}* is a constant that is the rate of the ensemble motion dictated by

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.

## DISCUSSION

### 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

In this equation, *t*/2 and then using Eq. 30 to generate the solution. Doing so shifts the actions of

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 ^{2}*J _{HH}* 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*-BuPCl

_{2}(Fig. 2B), which additionally exhibits a transition that is invariant under exchange and thus does not broaden (

*33*).

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

When calculating these spectra, we find that the traditional implementation and the DMEx converge to the same solution as Δ*t* → *dt*. 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 *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 ^{15}N polarization is calculated under the experimental conditions for SABRE–SHield Enables Alignment Transfer to Heteronuclei (SHEATH) (*16*).

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 ^{2}*J _{NH}* = −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 (

We begin with the dynamics of the free species. Exchange facilitates association of free species to the catalyst (*k _{N}*, 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 simply

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 (S* _{a}* and S

*), 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 manifolds*

_{b}Combining Eqs. 34 and 35 yields the equation of motion for the free species with exchange interactions

The bound species has two exchange interactions: one for the simultaneous exchange of a ligand and the hydrides occurring at rate of *k _{H}* and one for the exchange of target ligands at a rate of

*k*−

_{N}*k*. We will formulate the exchange operator for the bound species as a single entity,

_{H}*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

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 species

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 Δ*t* → *dt*. 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 σ

*≈ 0.5%.*

_{k}### 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-(^{15}N-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 ^{15}N 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 ^{1}H nuclei, with the largest approximation reported in literature using a single ^{1}H 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.

To emphasize the efficiency and flexibility of this framework, we used the DMExFR2 model to fit the coherent hyperpolarization dynamics of (^{15}N-^{13}C)-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 *k _{N}* = 14.5 ± 1.8 s

^{−1}and

*k*= 6.00 ± 0.75 s

_{H}^{−1}for the

*B*= − 1.65 μT data and

*k*= 15.0 ± 3.3 s

_{N}^{−1}and

*k*= 4.50 ± 0.98 s

_{H}^{−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 (^{15}N-^{13}C)-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.

## CONCLUSIONS AND OUTLOOK

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.

## MATERIALS AND METHODS

### Coherent SABRE-SHEATH experiments

A sample of (^{15}N,1-^{13}C)-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 MATERIALS

Supplementary material for this article is available at http://advances.sciencemag.org/cgi/content/full/6/32/eabb6874/DC1

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.

## REFERENCES AND NOTES

**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.

- Copyright © 2020 The Authors, some rights reserved; exclusive licensee American Association for the Advancement of Science. No claim to original U.S. Government Works. Distributed under a Creative Commons Attribution License 4.0 (CC BY).