Scalable modeling of multi-spin ensembles in SABRE hyperpolarization: a symmetry-based framework for zero and ultralow fields

This work presents a theoretical framework for quantitative, scalable modeling of signal amplification by reversible exchange (SABRE) experiments under zero- and ultralow-field (ZULF) conditions. SABRE exploits the singlet spin order of parahydrogen to hyperpolarize nuclear spins of substrates without chemical modification, enhancing NMR signals. In the ZULF SABRE method, polarization transfer occurs in ultralow magnetic fields where Zeeman interactions are comparable to or weaker than scalar couplings, enabling coherent mixing of spin states and revealing interactions often suppressed at high fields. Our approach captures the full quantum dynamics of SABRE, including coherent evolution, chemical exchange, and relaxation, within a Liouville space formalism. We demonstrate that the Hamiltonian, relaxation, and exchange superoperators possess symmetry with respect to the total spin, allowing the dynamics to be rigorously restricted to the zero-quantum coherence subspace. This symmetry-based reduction yields a scalable framework for efficient simulation of multi-spin SABRE systems, allowing the treatment of arbitrary spin ensembles, including those containing ¹⁵N, ¹³C, ¹H, and other nuclei. The approach is validated against full Liouville space calculations for small systems and is further applied to a 14-spin SABRE complex, demonstrating its ability to treat spin systems of a complexity well beyond the reach of conventional full Liouville space simulations. The framework thus provides a predictive tool for optimal polarization fields, ZULF NMR spectra, and the design of novel hyperpolarization experiments.

1 Introduction

Nuclear magnetic resonance (NMR) is a powerful spectroscopic technique that provides detailed information about molecular structure and dynamics. Its sensitivity, however, is fundamentally limited by the low thermal polarization of nuclear spins, a limitation that is particularly severe for heteronuclei whose low gyromagnetic ratios and natural abundances result in intrinsically weak NMR signals (Levitt, 2008). Enhancing the polarization levels, i.e. hyperpolarization, of the magnetic nuclei is therefore of central importance for extending the applicability of NMR to chemically and biologically relevant systems (Eills et al., 2023; Petersen et al., 2025; McBride et al., 2025; Kuhn et al., 2023; Zachrdla et al., 2025; Cavallari et al., 2018; Angelovski et al., 2023; Keshari and Wilson, 2014; Boutin et al., 2011).

Among the available hyperpolarization techniques, signal amplification by reversible exchange (SABRE) has emerged as a versatile and experimentally accessible method for boosting NMR sensitivity (Adams et al., 2009; Rayner and Duckett, 2018; Buntkowsky et al., 2022; Pravdivtsev et al., 2021; Atkinson et al., 2009). In SABRE, the spin order of parahydrogen is transferred to a target substrate through reversible binding to a transient metal–ligand complex, without permanent chemical modification of the molecule, as shown in Fig. 1A (Barskiy et al., 2019a; Pravdivtsev et al., 2013b; Shchepin et al., 2019). Importantly, SABRE has been shown to efficiently hyperpolarize a variety of heteronuclear spin systems, including ¹³C-, ¹⁵N-, ¹⁹F-, ³¹P-, and ⁷⁷Se-containing compounds (Kozinenko et al., 2025; Vaneeckhaute et al., 2023; Svyatova et al., 2019; Schmidt et al., 2023; Kiryutin et al., 2025; Markelov et al., 2024; Kiryutin et al., 2018; Burns et al., 2015; Olaru et al., 2018; Markelov et al., 2025; TomHon et al., 2025; Zheng et al., 2026). In biological NMR, heteronuclear detection is especially advantageous, as proton-based experiments are limited by the strong water background and short relaxation times (Mishra et al., 2016; Cho et al., 2019; Kim et al., 2025; Dey et al., 2020; Park and Wang, 2022).

Figure 1 Conceptual overview of the symmetry-based framework for SABRE under zero- and ultralow-field (ZULF) NMR conditions. (A) The full Liouville space has a dimension of 4^N, where N is the number of spins, rendering simulations of SABRE that include coherent spin evolution, relaxation, and chemical exchange computationally demanding. (B) Two symmetry reductions are employed: (1) permutation symmetry of equivalent spins, which is exploited via the introduction of effective spins K, and (2) conservation of the total spin projection F_z, which allows the dynamics to be rigorously restricted to the zero-quantum coherence (ZQC) subspace. (C) As a result, the problem is reduced to a substantially smaller ZQC Liouville space representation that still fully captures the relevant spin dynamics, thereby enabling quantitative simulations of multi-spin SABRE systems.

Importantly, the efficiency of SABRE hyperpolarization is governed by the coherent spin dynamics within the coupled spin network formed by the parahydrogen-derived hydrides and the nuclear spins of the substrate in the transient metal complex (Pravdivtsev et al., 2013b, 2015; Ivanov et al., 2014; Pravdivtsev et al., 2014a, b, c; Kiryutin et al., 2013). Despite growing experimental interest in SABRE, the number of rigorous theoretical models describing its spin dynamics remains relatively limited. In the present work, we adopt a fully rigorous description of SABRE spin dynamics based on a Liouville space formulation, in which coherent evolution, relaxation, and chemical exchange are treated on equal footing through a master equation. Although alternative approaches exist (Lindale et al., 2020), our focus here is specifically on a framework that captures the complete quantum dynamics of the system (Knecht et al., 2020; Knecht and Ivanov, 2019; Knecht et al., 2016). This approach is firmly rooted in the standard theory of NMR spin dynamics (Limbach, 1991; Abergel and Palmer, 2005; Cuperlovic et al., 2000; Zaiss and Bachert, 2013; Schwartz et al., 1982; Jameson and Bruschweiler, 2021) and provides a conceptually transparent description of SABRE polarization transfer. However, the dimensionality of Liouville space grows exponentially with the number of spins involved, making direct simulations computationally intractable for large heteronuclear SABRE systems; see Fig. 1A. Consequently, the development of a general and scalable framework for such multi-spin ensembles remains a significant challenge. Motivated by this, we present in this work a complementary approach to recent efforts (Mamone et al., 2025) that leverages symmetry to render the full Liouville space treatment computationally tractable for complex SABRE systems.

In parallel with the development of hyperpolarization methods, zero- and ultralow-field (ZULF) NMR has been gaining increasing attention as a rapidly developing complement to conventional high-field NMR that enables access to spin dynamics and interactions inaccessible at high magnetic fields (Blanchard and Budker, 2016; Blanchard et al., 2020; Ledbetter et al., 2011, 2009; Put et al., 2021; Kiryutin et al., 2021; Barskiy et al., 2025, 2019b; Xu et al., 2026; Sheberstov et al., 2021). In this regime, experiments are performed in magnetic fields where Zeeman interactions are comparable to or weaker than scalar J couplings. ZULF NMR therefore offers several distinctive advantages that have driven its increasing adoption. The absence of a strong static magnetic field eliminates line broadening due to magnetic field inhomogeneities and enables exceptionally high spectral resolution (Blanchard et al., 2013; Picazo-Frutos et al., 2024; Sjolander et al., 2020). Moreover, ZULF NMR is inherently compatible with compact, magnetically shielded setups and non-inductive detection schemes based on optically pumped magnetometers, allowing high sensitivity without the need for superconducting magnets (Ledbetter et al., 2009; Jiang et al., 2019; Hong et al., 2025). For heteronuclear systems in particular, ZULF NMR provides access to spin interactions and dynamical effects that are strongly suppressed or effectively truncated at high magnetic fields (Blanchard et al., 2015; Teleanu et al., 2025; Pravdivtsev et al., 2013a). In the ZULF regime, heteronuclear spin evolution is governed by the full, non-secular spin Hamiltonian, allowing otherwise negligible coupling terms and coherent mixing of spin states to influence the observed spectra (Kiryutin et al., 2021; Zhukov et al., 2020, 2021). This makes ZULF NMR especially well-suited for probing heteronuclear spin dynamics beyond the high-field approximation.

The combination of the parahydrogen-based hyperpolarization with zero- and ultralow-field NMR provides a powerful approach for enhancing the sensitivity of NMR detection (Theis et al., 2012, 2011; Buckenmaier et al., 2025; Rodriguez et al., 2025; Van Dyke et al., 2022; Put et al., 2023). By leveraging hyperpolarization in the ultralow-field regime, enhanced signals can be observed from collective spin states comprising ¹H, ¹³C, ¹⁵N, and other nuclei. The theoretical consideration of SABRE under the ZULF conditions, or ZULF SABRE, is helpful for guiding the design of experiments, estimating optimal fields for polarization buildup, and predicting ZULF NMR spectra, thereby supporting the development of new polarization transfer strategies and their adaptation to the unique conditions of SABRE systems.

Theoretical and computational studies of SABRE spin dynamics have attracted significant attention in recent years. In a number of works, it has been shown that the key contributions to polarization transfer in SABRE can be understood in terms of zero-quantum operators, highlighting the central role of the zero-quantum coherence (ZQC) subspace (Eriksson et al., 2022; Snadin et al., 2024; Markelov et al., 2021; Li et al., 2022). In Liouville space, this subspace is defined by operators with zero eigenvalue of the superoperator of the z projection of the total spin, ${\widehat{\widehat{\mathbf{F}}}}_z$, where the quantization axis is set by the direction of the ultralow external magnetic field. In particular, Mamone et al. (2025) demonstrated that the description in terms of zero-quantum operators, including exchange and relaxation, successfully captures SABRE dynamics while reducing the dimensionality of the problem.

Here, we place these observations on a rigorous footing. We formally show (Fig. 1B) that the Hamiltonian, relaxation, and chemical exchange superoperators in Liouville space possess a well-defined symmetry with respect to ${\widehat{\widehat{\mathbf{F}}}}_z$. This symmetry guarantees that the dynamics is restricted to the ZQC subspace and that only these components contribute to polarization transfer in ZULF SABRE. Exploiting this symmetry, we develop a general and scalable computational framework based on systematic restriction to the ZQC subspace (Fig. 1C). The approach is validated against full Liouville space simulations for small spin systems, reproducing exact results while providing computational speedups of 30–50 times. For larger systems, such as a 14-spin SABRE complex, the method enables simulations that would be computationally infeasible in the full Liouville space.

2 Theory

2.1 Block-diagonal decomposition of Hilbert space

We begin by considering the low-field spin Hamiltonian governing the coherent dynamics of a coupled N-spin system (in units of ℏ):

$$\begin{array}{}\text{(1)}& \widehat{H}\left(t\right)=-B_{\mathrm{0}}\left(t\right)\sum _{l=\mathrm{1}}^N{\mathit{\gamma }}^l{\widehat{I}}_z^l+\mathrm{2}\mathit{\pi }\sum _{l<m}^N{J}^{lm}({\widehat{\mathbf{I}}}^l\cdot {\widehat{\mathbf{I}}}^m),\end{array}$$

where N is the total number of nuclei in the spin system, B₀(t) is an external magnetic field, γ^l is the gyromagnetic ratio of nucleus l, J^lm is the J-coupling constant between nuclei l and m, and $({\widehat{\mathbf{I}}}^l\cdot {\widehat{\mathbf{I}}}^m)$ denotes the scalar product of the corresponding spin operators. At ultralow magnetic fields, chemical shift terms are neglected. However, the framework presented below is not limited to this regime; chemical shifts can be incorporated straightforwardly without modifying the theoretical description.

In this section, we analyze the symmetry properties of the Hamiltonian in Eq. (1) and show that they give rise to a block-diagonal structure in an appropriate basis. As a result, the spin dynamics can be rigorously restricted to invariant subspaces of reduced dimensionality. This symmetry-based reduction provides a practical route for treating large, coupled spin systems and underpins the SABRE applications presented in the subsequent sections.

Importantly, throughout this work we treat the general and experimentally relevant case of a non-zero magnetic field, B₀≠0. In practice, a strictly zero field is never realized due to residual magnetic fields, which makes the present approach directly applicable to SABRE polarization transfer at ultralow fields and to the simulation of ZULF NMR spectra across different coupling regimes, ranging from the J-coupling-dominated to the Zeeman-dominated limit (Buckenmaier et al., 2025). A separate and careful symmetry analysis would be required only in the idealized limit B₀=0, which lies beyond the scope of the present study.

2.1.1 The group of magnetically equivalent nuclei

Assume that we have a group G_P of P magnetically equivalent spins, reflecting the underlying molecular symmetry, for instance, the three protons of a methyl (CH₃) group or the two protons of a methylene (CH₂) group. Therefore, for each external nucleus m∉G_P, the J coupling between m and any spin p∈G_P is the same: $J^{mp}=J^{m{G}_P}$. Thus, the J-coupling constant is independent of the index of the nucleus within G_P and depends only on the outer nucleus m. The magnetic equivalence enables us to treat the group G_P as a pseudo-nucleus with total spin K. The Hilbert space of P identical spin-$\mathrm{1}/\mathrm{2}$ nuclei can be decomposed into invariant subspaces corresponding to definite total spin (Messiah, 1962):

$$\begin{array}{}\text{(2)}& \underset{i=\mathrm{1}}{\overset{P}{\otimes }}{D}_{\mathrm{1}/\mathrm{2}}^i=\underset{K=K_{min}}{\overset{P/\mathrm{2}}{\oplus }}{\mathit{\mu }}_K{D}_K,\end{array}$$

where $D_{\mathrm{1}/\mathrm{2}}^i$ denotes the spin-$\mathrm{1}/\mathrm{2}$ representation for spin i, D_K is the irreducible representation of the collective spin system corresponding to total spin K, and

$$\begin{array}{}\text{(3)}& {\mathit{\mu }}_K={\displaystyle \frac{\mathrm{2}K+\mathrm{1}}{P/\mathrm{2}+K+\mathrm{1}}}\left(\genfrac{}{}{0ex}{}{P}{P/\mathrm{2}-K}\right)\end{array}$$

is the multiplicity of spin K, where $\left(\genfrac{}{}{0ex}{}{n}{k}\right)=n\mathrm{!}/\left[\right(n-k\left)\mathrm{!}k\mathrm{!}\right]$ is the binomial coefficient. This commutation relation, $[{\widehat{K}}^{\mathrm{2}},\widehat{H}(t\left)\right]=\mathrm{0}$, allows the Hamiltonian in Eq. (1) to be block diagonalized according to the total spin K of the group G_P (Barskiy and Pravdivtsev, 2025):

$$\begin{array}{}\text{(4)}& \widehat{H}\left(t\right)=\underset{K=K_{min}}{\overset{P/\mathrm{2}}{\oplus }}\langle K\left|\widehat{H}\right(t\left)\right|K\rangle =\underset{K=K_{min}}{\overset{P/\mathrm{2}}{\oplus }}{\widehat{H}}_K\left(t\right),\end{array}$$

where

$$\begin{array}{}\text{(5)}& K_{min}=\left\{\begin{array}{ll}\mathrm{0}& \text{if }P\text{ is even},\\ \mathrm{1}/\mathrm{2}& \text{if }P\text{ is odd}.\end{array}\right.\end{array}$$

Similarly, when the density matrix of the full spin ensemble commutes with the total spin operator, $[{\widehat{K}}^{\mathrm{2}},\widehat{\mathit{\rho }}(t\left)\right]=\mathrm{0}$, it can be decomposed into the corresponding invariant subspaces:

$$\begin{array}{}\text{(6)}& \widehat{\mathit{\rho }}\left(t\right)=\underset{K=K_{min}}{\overset{P/\mathrm{2}}{\oplus }}{g}_K{\widehat{\mathit{\rho }}}_K\left(t\right),\end{array}$$

where

$$\begin{array}{}\text{(7)}& g_K={\displaystyle \frac{(\mathrm{2}K+\mathrm{1}){\mathit{\mu }}_K}{{\mathrm{2}}^P}}\end{array}$$

is a statistical weight of the K subspace with ${\sum }_{K=K_{min}}^{P/\mathrm{2}}{g}_K=\mathrm{1}$. Each block of the density matrix is normalized in the usual way, $\text{Tr}\mathit{\left\{}{\widehat{\mathit{\rho }}}_K\right(t\left)\mathit{\right\}}=\mathrm{1}$. For an observable $\widehat{O}$ that commutes with the total spin operator, $[{\widehat{K}}^{\mathrm{2}},\widehat{O}]=\mathrm{0}$, the expectation value can be calculated block-wise. First, the operator is decomposed into the following corresponding blocks:

$$\begin{array}{}\text{(8)}& \widehat{O}=\underset{K=K_{min}}{\overset{P/\mathrm{2}}{\oplus }}{\widehat{O}}_K,\end{array}$$

and then the ensemble average is obtained as

$$\begin{array}{}\text{(9)}& \begin{array}{rl}\langle \widehat{O}\rangle \left(t\right)& =\text{Tr}\mathit{\left\{}\widehat{O}\widehat{\mathit{\rho }}\right(t\left)\mathit{\right\}}=\sum _{K=K_{min}}^{P/\mathrm{2}}{g}_K\text{Tr}\mathit{\left\{}{\widehat{O}}_K{\widehat{\mathit{\rho }}}_K\right(t\left)\mathit{\right\}}\\ & =\sum _{K=K_{min}}^{P/\mathrm{2}}{g}_K\langle {\widehat{O}}_K\rangle \left(t\right).\end{array}\end{array}$$

2.1.2 The z projection of total spin

Each block ${\widehat{H}}_K\left(t\right)$ corresponding to the total spin K of the group G_P can be further decomposed according to the z projection of the total spin. Specifically, consider a single K block:

$$\begin{array}{}\text{(10)}& \begin{array}{rl}{\widehat{H}}_K\left(t\right)& =-B_{\mathrm{0}}\left(t\right)\left(\sum _{l=\mathrm{1}}^{N-P}{\mathit{\gamma }}^l{\widehat{I}}_z^l+{\mathit{\gamma }}^P{\widehat{K}}_z\right)\\ & +\mathrm{2}\mathit{\pi }\sum _{l<m}^{N-P}{J}^{lm}({\widehat{\mathbf{I}}}^l\cdot {\widehat{\mathbf{I}}}^m)+\mathrm{2}\mathit{\pi }\sum _{l=\mathrm{1}}^{N-P}{J}^{l{G}_P}({\widehat{\mathbf{I}}}^l\cdot \widehat{\mathbf{K}}),\end{array}\end{array}$$

where γ^P is the gyromagnetic ratio of a nucleus in the group G_P of P magnetically equivalent nuclei, and $J^{l{G}_P}$ denotes the J-coupling constant between an external nucleus l∉G_P and the nuclei within the group G_P. J couplings within the group G_P are omitted as they do not affect the dynamics within the K block. The Hamiltonian ${\widehat{H}}_K\left(t\right)$ commutes with the z projection of the total spin of all nuclei, $[{\widehat{F}}_z,{\widehat{H}}_K(t\left)\right]=\mathrm{0}$, where

$$\begin{array}{}\text{(11)}& {\widehat{F}}_z=\sum _{l=\mathrm{1}}^{N-P}{\widehat{I}}_z^l+{\widehat{K}}_z,\end{array}$$

and the first sum runs over a set of non-equivalent spin-$\mathrm{1}/\mathrm{2}$ nuclei outside the group G_P. Consequently, ${\widehat{H}}_K\left(t\right)$ is block diagonal with respect to the eigenvalues of ${\widehat{F}}_z$, which takes values $F_z\in \left\{\frac{N-P}{\mathrm{2}}+K,\frac{N-P}{\mathrm{2}}+K-\mathrm{1},\mathrm{\dots },-\frac{N-P}{\mathrm{2}}-K\right\}$. In other words, ${\widehat{H}}_K\left(t\right)$ can be decomposed as $(\mathrm{2}K+\mathrm{1}+N-P)$ blocks:

$$\begin{array}{}\text{(12)}& \begin{array}{rl}{\widehat{H}}_K\left(t\right)& =\underset{F_z=-(N-P)/\mathrm{2}-K}{\overset{(N-P)/\mathrm{2}+K}{\oplus }}\langle F_z\left|{\widehat{H}}_K\right(t\left)\right|F_z\rangle \\ & =\underset{F_z=-(N-P)/\mathrm{2}-K}{\overset{(N-P)/\mathrm{2}+K}{\oplus }}{\widehat{H}}_{K,F_z}\left(t\right),\end{array}\end{array}$$

with each block ${\widehat{H}}_{K,F_z}\left(t\right)$ corresponding to specific values of K and F_z.

2.1.3 Coherent spin dynamics

The time evolution of the density matrix $\widehat{\mathit{\rho }}\left(t\right)$ is governed by the Liouville–von Neumann (LvN) equation:

$$\begin{array}{}\text{(13)}& {\displaystyle \frac{\partial \widehat{\mathit{\rho }}\left(t\right)}{\partial t}}=-i\left[\widehat{H}\right(t),\widehat{\mathit{\rho }}(t\left)\right],\end{array}$$

with an initial state $\widehat{\mathit{\rho }}\left(\mathrm{0}\right)$. If the initial state commutes with the symmetry operators, $[{\widehat{K}}^{\mathrm{2}},\widehat{\mathit{\rho }}(\mathrm{0}\left)\right]=[{\widehat{F}}_z,\widehat{\mathit{\rho }}(\mathrm{0}\left)\right]=\mathrm{0}$, then the block-diagonal decomposition of $\widehat{\mathit{\rho }}\left(\mathrm{0}\right)$ is valid:

$$\begin{array}{}\text{(14)}& \begin{array}{rl}\widehat{\mathit{\rho }}\left(\mathrm{0}\right)& =\underset{K=K_{min}}{\overset{P/\mathrm{2}}{\oplus }}{g}_K{\widehat{\mathit{\rho }}}_K\left(\mathrm{0}\right)\\ & =\underset{K=K_{min}}{\overset{P/\mathrm{2}}{\oplus }}\underset{F_z=-(N-P)/\mathrm{2}-K}{\overset{(N-P)/\mathrm{2}+K}{\oplus }}{g}_K{\widehat{\mathit{\rho }}}_{K,F_z}\left(\mathrm{0}\right),\end{array}\end{array}$$

where we introduced the reduced density matrices ${\widehat{\mathit{\rho }}}_{K,F_z}\left(t\right)$. The dimension of the F_z block in the Hilbert space ℋ, denoted as $\text{dim}\left({\mathcal{H}}_{F_z}\right)$, is given by

$$\begin{array}{}\text{(15)}& \text{dim}\left({\mathcal{H}}_{F_z}\right)=\sum _{m_K=-K}^K\left(\genfrac{}{}{0ex}{}{N-P}{F_z+(N-P)/\mathrm{2}-m_K}\right).\end{array}$$

Moreover, if the Hamiltonian commutes with the same symmetry operators,

$$\begin{array}{}\text{(16)}& \begin{array}{rl}\text{If }[{\widehat{K}}^{\mathrm{2}},& \widehat{H}\left(t\right)]=[{\widehat{F}}_z,\widehat{H}\left(t\right)]=\mathrm{0},\\ \text{then }& [{\widehat{K}}^{\mathrm{2}},\widehat{\mathit{\rho }}(t\left)\right]=[{\widehat{F}}_z,\widehat{\mathit{\rho }}(t\left)\right]=\mathrm{0}\text{ for any }t.\end{array}\end{array}$$

This means that the block-diagonal structure of the density matrix is conserved for all instants of time:

$$\begin{array}{}\text{(17)}& \widehat{\mathit{\rho }}\left(t\right)=\underset{K=K_{min}}{\overset{P/\mathrm{2}}{\oplus }}\underset{F_z=-(N-P)/\mathrm{2}-K}{\overset{(N-P)/\mathrm{2}+K}{\oplus }}{g}_K{\widehat{\mathit{\rho }}}_{K,F_z}\left(t\right)\text{ for any }t.\end{array}$$

The reduced density matrices ${\widehat{\mathit{\rho }}}_{K,F_z}\left(t\right)$ obey the LvN equation independently within each block:

$$\begin{array}{}\text{(18)}& {\displaystyle \frac{\partial {\widehat{\mathit{\rho }}}_{K,F_z}\left(t\right)}{\partial t}}=-i\left[{\widehat{H}}_{K,F_z}\right(t),{\widehat{\mathit{\rho }}}_{K,F_z}(t\left)\right].\end{array}$$

Thus, the block-diagonal decomposition reduces the computational complexity by considering smaller blocks with fixed values of K and F_z independently.

2.2 Block-diagonal decomposition of Liouville space

2.2.1 Coherent spin dynamics

In Liouville space, an operator $\widehat{A}$ from Hilbert space is mapped to a commutation superoperator $\widehat{\widehat{\mathbf{A}}}$ acting on an operator $\widehat{O}$ according to

$$\begin{array}{}\text{(19)}& \widehat{\widehat{\mathbf{A}}}\widehat{O}=[\widehat{A},\widehat{O}].\end{array}$$

At the same time, the density matrix in Liouville space is represented as a vector; see Fig. 2A. Throughout this paper, operators are denoted by hatted symbols, whereas superoperators are denoted by bold double-hatted symbols. Under this mapping, the Hamiltonian $\widehat{H}\left(t\right)$ is represented by the Liouville space superoperator $\widehat{\widehat{\mathbf{H}}}\left(t\right)$, and similarly, the effective-spin operator ${\widehat{K}}^{\mathrm{2}}$ is mapped to the superoperator ${\widehat{\widehat{\mathbf{K}}}}^{\mathrm{2}}$. In Liouville space, the coherent spin dynamics is described by the LvN equation:

$$\begin{array}{}\text{(20)}& {\displaystyle \frac{\partial \widehat{\mathit{\rho }}\left(t\right)}{\partial t}}=-i\widehat{\widehat{\mathbf{H}}}\left(t\right)\widehat{\mathit{\rho }}\left(t\right).\end{array}$$

If the initial density matrix commutes with the effective-spin operator ${\widehat{K}}^{\mathrm{2}}$, i.e.,

$$\begin{array}{}\text{(21)}& \text{If }[{\widehat{K}}^{\mathrm{2}},\widehat{\mathit{\rho }}(\mathrm{0}\left)\right]=\mathrm{0},\text{ then }{\widehat{\widehat{\mathbf{K}}}}^{\mathrm{2}}\widehat{\mathit{\rho }}\left(\mathrm{0}\right)=\mathrm{0}.\end{array}$$

Moreover, Liouville space conserves commutation relations so that

$$\begin{array}{}\text{(22)}& \text{If }[{\widehat{K}}^{\mathrm{2}},\widehat{H}(t\left)\right]=\mathrm{0},\text{ then }[{\widehat{\widehat{\mathbf{K}}}}^{\mathrm{2}},\widehat{\widehat{\mathbf{H}}}(t\left)\right]=\mathrm{0}.\end{array}$$

Consequently, the block-diagonal structure of the density matrix is conserved during the evolution

$$\begin{array}{}\text{(23)}& \text{so that }{\widehat{\widehat{\mathbf{K}}}}^{\mathrm{2}}\widehat{\mathit{\rho }}\left(t\right)=\mathrm{0},\text{ meaning }[{\widehat{K}}^{\mathrm{2}},\widehat{\mathit{\rho }}(t\left)\right]=\mathrm{0}.\end{array}$$

Therefore, the blocks corresponding to different total spin values K remain uncoupled. Analogous to the Hilbert space treatment, in Liouville space we can also treat the group G_P of P magnetically equivalent nuclear spins as a pseudo-nucleus with total spin $K\in \left\{\frac{P}{\mathrm{2}},\frac{P}{\mathrm{2}}-\mathrm{1},\mathrm{\dots },K_{min}\right\}.$

Figure 2 Schematic illustration of the zero-quantum coherence (ZQC) reduction in Liouville space. (A) Transformation of the density matrix $\widehat{\mathit{\rho }}\left(t\right)$ from Hilbert space to Liouville space. (B) Decomposition of the Hamiltonian and the relaxation superoperators according to the coherence order, i.e., the eigenvalues of ${\widehat{\widehat{\mathbf{F}}}}_z$: +2 (double-quantum coherence, +DQC), +1 (single-quantum coherence, +SQC), 0 (zero-quantum coherence, ZQC), and the corresponding negative coherence orders. (C) The initial state of the system is assumed to satisfy ${\widehat{\widehat{\mathbf{F}}}}_z\widehat{\mathit{\rho }}\left(\mathrm{0}\right)=\mathrm{0}$ so that only ZQC operators are initially populated. (D) Only the ZQC subspace is relevant for spin evolution under the total Liouvillian $\widehat{\widehat{\mathbf{L}}}=-i\widehat{\widehat{\mathbf{H}}}+\widehat{\widehat{\mathbf{\Gamma }}}$, implying that ${\widehat{\widehat{\mathbf{F}}}}_z\widehat{\mathit{\rho }}\left(t\right)=\mathrm{0}$, and $\widehat{\mathit{\rho }}\left(t\right)$ has no components outside the ZQC subspace. Consequently, all other coherence order subspaces can be omitted without affecting the spin dynamics.

Consequently, the dynamics can be considered independently for each K subspace, and the expectation value of an observable is obtained by averaging over the blocks according to Eqs. (8) and (9), as in Hilbert space. The time evolution of the reduced density matrix ${\widehat{\mathit{\rho }}}_K\left(t\right)$ in Liouville space is then given by

$$\begin{array}{}\text{(24)}& {\displaystyle \frac{\partial {\widehat{\mathit{\rho }}}_K\left(t\right)}{\partial t}}=-i{\widehat{\widehat{\mathbf{H}}}}_K\left(t\right){\widehat{\mathit{\rho }}}_K\left(t\right),\end{array}$$

where we introduced the Hamiltonian superoperator reduced onto the K subspace as ${\widehat{\widehat{\mathbf{H}}}}_K\left(t\right)•=\left[{\widehat{H}}_K\right(t),•]$, and the density matrix is normalized as $\text{Tr}\mathit{\left\{}{\widehat{\mathit{\rho }}}_K\right(t\left)\mathit{\right\}}=\mathrm{1}$.

The symmetry associated with ${\widehat{F}}_z$ in Hilbert space, see Eq. (11), is directly inherited in Liouville space. To analyze this, the operator ${\widehat{F}}_z$ is mapped to the corresponding superoperator ${\widehat{\widehat{\mathbf{F}}}}_z$ using Eq. (19). The eigenvalues of ${\widehat{\widehat{\mathbf{F}}}}_z$ are commonly referred to as the coherence orders, which are well-known in standard high-field NMR; see Appendix A for details. A key property of the Liouville space evolution of Eq. (24) is

$$\begin{array}{}\text{(25)}& \begin{array}{rl}\text{If }{\widehat{\widehat{\mathbf{F}}}}_z{\widehat{\mathit{\rho }}}_K\left(\mathrm{0}\right)& =\mathrm{0}\text{ and }[{\widehat{\widehat{\mathbf{F}}}}_z,{\widehat{\widehat{\mathbf{H}}}}_K(t\left)\right]=\mathrm{0},\\ & \text{then }{\widehat{\widehat{\mathbf{F}}}}_z{\widehat{\mathit{\rho }}}_K\left(t\right)=\mathrm{0}\text{ for any }t.\end{array}\end{array}$$

Equation (25) is schematically illustrated in Fig. 2B–D. In other words, the coherent dynamics in Liouville space conserves evolution strictly within the subspace of zero eigenvalue of ${\widehat{\widehat{\mathbf{F}}}}_z$, i.e., the zero-quantum coherence (ZQC) block. The ZQC subspace is spanned by operators of the following form, as discussed in Appendix A:

$$\begin{array}{}\text{(26)}& \text{ZQC}=\text{span}\mathit{\left\{}\right|F_z,\mathit{\xi }\rangle \langle F_z,\mathit{\lambda }\left|\mathit{\right\}},\end{array}$$

with $F_z\in \left\{\frac{N-P}{\mathrm{2}}+K,\frac{N-P}{\mathrm{2}}+K-\mathrm{1},\mathrm{\dots },-\frac{N-P}{\mathrm{2}}-K\right\}$ and ξ,λ run over all possible magnetic quantum numbers of the individual nuclear spins for a given value of the z projection of total spin F_z.

The coherent evolution within the ZQC block is then governed by the standard LvN equation; see Eq. (24):

$$\begin{array}{}\text{(27)}& {\displaystyle \frac{\partial {\widehat{\mathit{\rho }}}_K\left(t\right)}{\partial t}}=-i{\widehat{\widehat{\mathbf{H}}}}_K^{\text{ZQC}}\left(t\right){\widehat{\mathit{\rho }}}_K\left(t\right),\end{array}$$

where the dynamics is restricted to the ZQC subspace, as illustrated in Fig. 2D. Here, we also introduced ${\widehat{\widehat{\mathbf{H}}}}_K^{\text{ZQC}}\left(t\right)$ as the Hamiltonian superoperator projected onto the ZQC subspace. The subscript ZQC of the density matrix is omitted since ${\widehat{\mathit{\rho }}}_K\left(t\right)\equiv \mathrm{0}$ outside the ZQC subspace; see Fig. 2D.

We also emphasize that the consideration of the ZQC subspace drastically reduces the dimensionality of the problem. The full Liouville space has a dimension of 4^N, whereas the ZQC subspace with the effective-spin reduction has a smaller dimension:

$$\begin{array}{}\text{(28)}& \begin{array}{rl}\text{dim}& \left(\text{ZQC}\right)\\ & =\sum _{F_z=-(N-P)/\mathrm{2}-K}^{(N-P)/\mathrm{2}+K}{\left\{\text{dim}\left({\mathcal{H}}_{F_z}\right)\right\}}^{\mathrm{2}}\\ & =\sum _{F_z=-(N-P)/\mathrm{2}-K}^{(N-P)/\mathrm{2}+K}{\left\{\sum _{m_K=-K}^K\left(\genfrac{}{}{0ex}{}{N-P}{F_z+(N-P)/\mathrm{2}-m_K}\right)\right\}}^{\mathrm{2}}\\ & <{\mathrm{4}}^N,\end{array}\end{array}$$

where Eq. (15) was used to calculate $\text{dim}\left({\mathcal{H}}_{F_z}\right)$. This restriction to the ZQC subspace allows efficient treatment of spin dynamics by working with a smaller independent block rather than the full Liouville space.

2.2.2 Relaxational spin dynamics

Relaxational spin dynamics is described within the Redfield formalism, which leads to the following equation of motion for the density matrix $\widehat{\mathit{\rho }}\left(t\right)$:

$$\begin{array}{}\text{(29)}& {\displaystyle \frac{\partial \widehat{\mathit{\rho }}\left(t\right)}{\partial t}}=(-i\widehat{\widehat{\mathbf{H}}}(t)+\widehat{\widehat{\mathbf{\Gamma }}})\widehat{\mathit{\rho }}\left(t\right)=\widehat{\widehat{\mathbf{L}}}\left(t\right)\widehat{\mathit{\rho }}\left(t\right),\end{array}$$

where $\widehat{\widehat{\mathbf{H}}}\left(t\right)$ is the coherent Hamiltonian superoperator, $\widehat{\widehat{\mathbf{\Gamma }}}$ is the relaxational superoperator describing stochastically modulated dynamics, and $\widehat{\widehat{\mathbf{L}}}\left(t\right)=-i\widehat{\widehat{\mathbf{H}}}\left(t\right)+\widehat{\widehat{\mathbf{\Gamma }}}$ is the total Liouvillian. In this work, we used the random fluctuating field (RFF) mechanism in the fast-motion limit, also known as the extreme narrowing regime. Within this approximation, the relaxation superoperator takes the following form (Pileio, 2010):

$$\begin{array}{}\text{(30)}& \widehat{\widehat{\mathbf{\Gamma }}}=-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\sum _{n,l}^N{\displaystyle \frac{C^{nl}}{\sqrt{T_{\mathrm{1}}^n{T}_{\mathrm{1}}^l}}}\sum _{m=-\mathrm{1}}^{\mathrm{1}}(-\mathrm{1}{)}^m{\widehat{\widehat{\mathbf{T}}}}_{\mathrm{1},-m}^n{\widehat{\widehat{\mathbf{T}}}}_{\mathrm{1},m}^l,\end{array}$$

where we assumed that the local fields are isotropic and introduced the following quantities: $T_{\mathrm{1}}^n$ is the longitudinal T₁-relaxation time of the nucleus n, $C^{nl}\in [\mathrm{0},\mathrm{1}]$ is a correlation constant between the local fluctuating magnetic fields at the spatial positions of nuclear spins n and l (Freeman et al., 1970), and ${\widehat{\widehat{\mathbf{T}}}}_{\mathrm{1},m}^n\mathbf{•}=[{\widehat{T}}_{\mathrm{1},m}^n,\mathbf{•}]$ is the rank-1 irreducible spherical superoperator acting on nucleus n. The corresponding spherical tensor operators are defined as

$$\begin{array}{}\text{(31)}& {\widehat{T}}_{\mathrm{1},\mathrm{1}}^n=-{\displaystyle \frac{{\widehat{I}}_x^n+i{\widehat{I}}_y^n}{\sqrt{\mathrm{2}}}},{\widehat{T}}_{\mathrm{1},\mathrm{0}}^n={\widehat{I}}_z^n,{\widehat{T}}_{\mathrm{1},-\mathrm{1}}^n={\displaystyle \frac{{\widehat{I}}_x^n-i{\widehat{I}}_y^n}{\sqrt{\mathrm{2}}}}.\end{array}$$

For the spin systems containing a group G_P of P magnetically equivalent nuclei (e.g. CH₂ or CH₃ groups), Eq. (30) can be rewritten in a more convenient form by separating contributions from equivalent and non-equivalent spins:

$$\begin{array}{}\text{(32)}& \begin{array}{rl}\widehat{\widehat{\mathbf{\Gamma }}}& =-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\sum _{n\notin G_P}{\displaystyle \frac{\mathrm{1}}{T_{\mathrm{1}}^n}}\sum _{m=-\mathrm{1}}^{\mathrm{1}}(-\mathrm{1}{)}^m{\widehat{\widehat{\mathbf{T}}}}_{\mathrm{1},-m}^n{\widehat{\widehat{\mathbf{T}}}}_{\mathrm{1},m}^l\\ & -{\displaystyle \frac{\mathrm{1}}{\mathrm{2}{T}_{\mathrm{1}}^P}}\sum _{n,l\in G_P}\sum _{m=-\mathrm{1}}^{\mathrm{1}}(-\mathrm{1}{)}^m{\widehat{\widehat{\mathbf{T}}}}_{\mathrm{1},-m}^n{\widehat{\widehat{\mathbf{T}}}}_{\mathrm{1},m}^l,\end{array}\end{array}$$

where we have explicitly accounted for the equivalence of the nuclei in group G_P. Specifically, all nuclei in G_P are assumed to have identical longitudinal relaxation times $T_{\mathrm{1}}^P$, and the local fluctuating magnetic fields at their spatial positions are taken to be fully correlated, such that C^nl=1 for $n,l\in G_P$. All other local fields are assumed to be uncorrelated; i.e., C^nl=0 if either n∉G_P or l∉G_P.

The second term in Eq. (32), and therefore the full relaxation superoperator $\widehat{\widehat{\mathbf{\Gamma }}}$, commutes with the superoperator of the total spin of the group G_P, ${\widehat{\widehat{\mathbf{K}}}}^{\mathrm{2}}$. This follows from the identity

$$\begin{array}{}\text{(33)}& \sum _{n,l\in G_P}\sum _{m=-\mathrm{1}}^{\mathrm{1}}(-\mathrm{1}{)}^m{\widehat{\widehat{\mathbf{T}}}}_{\mathrm{1},-m}^n{\widehat{\widehat{\mathbf{T}}}}_{\mathrm{1},m}^l=\sum _{n,l\in G_P}({\widehat{\widehat{\mathbf{T}}}}_{\mathrm{1}}^n\cdot {\widehat{\widehat{\mathbf{T}}}}_{\mathrm{1}}^l),\end{array}$$

which represents a scalar product of rank-1 spherical superoperators in Liouville space (Pyper, 1971). This is fully analogous to the Hilbert space, where the scalar product of equivalent nuclear spins, ${\sum }_{n,l\in G_P}({\widehat{\mathbf{I}}}^n\cdot {\widehat{\mathbf{I}}}^l)$, commutes with ${\widehat{K}}^{\mathrm{2}}$. As a consequence, $[{\widehat{\widehat{\mathbf{K}}}}^{\mathrm{2}},\widehat{\widehat{\mathbf{\Gamma }}}]=\mathrm{0}$, and the relaxation superoperator can be decomposed into independent blocks corresponding to different values of ${\widehat{\widehat{\mathbf{K}}}}^{\mathrm{2}}$. Importantly, we are interested only in the zero-eigenvalue subspace of ${\widehat{\widehat{\mathbf{K}}}}^{\mathrm{2}}$ since

$$\begin{array}{}\text{(34)}& \begin{array}{rl}\text{If }{\widehat{\widehat{\mathbf{K}}}}^{\mathrm{2}}\widehat{\mathit{\rho }}\left(\mathrm{0}\right)=\mathrm{0}\text{ and }& [{\widehat{\widehat{\mathbf{K}}}}^{\mathrm{2}},\widehat{\widehat{\mathbf{\Gamma }}}]=[{\widehat{\widehat{\mathbf{K}}}}^{\mathrm{2}},\widehat{\widehat{\mathbf{H}}}(t\left)\right]=\mathrm{0},\\ & \text{then }{\widehat{\widehat{\mathbf{K}}}}^{\mathrm{2}}\widehat{\mathit{\rho }}\left(t\right)=\mathrm{0}\text{ for any }t.\end{array}\end{array}$$

Therefore, relaxational dynamics also conserves the density matrix $\widehat{\mathit{\rho }}\left(t\right)$ within the zero-eigenvalue subspace of ${\widehat{\widehat{\mathbf{K}}}}^{\mathrm{2}}$; i.e. ${\widehat{\widehat{\mathbf{K}}}}^{\mathrm{2}}\widehat{\mathit{\rho }}\left(t\right)=[{\widehat{K}}^{\mathrm{2}},\widehat{\mathit{\rho }}(t\left)\right]=\mathrm{0}$. As a consequence, the blocks of the density matrix corresponding to different values of the total spin K of the magnetically equivalent nuclei in the group G_P are not mixed by relaxation. This allows us to treat relaxation independently for each value of the nuclear spin $K\in \left\{\frac{P}{\mathrm{2}},\frac{P}{\mathrm{2}}-\mathrm{1},\mathrm{\dots },K_{min}\right\}$ of the group G_P (Doronin et al., 2011) and introduce the reduced relaxation superoperator

$$\begin{array}{}\text{(35)}& \begin{array}{rl}{\widehat{\widehat{\mathbf{\Gamma }}}}_K& =-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\sum _{n\notin G_P}{\displaystyle \frac{\mathrm{1}}{T_{\mathrm{1}}^n}}\sum _{m=-\mathrm{1}}^{\mathrm{1}}(-\mathrm{1}{)}^m{\widehat{\widehat{\mathbf{T}}}}_{\mathrm{1},-m}^n{\widehat{\widehat{\mathbf{T}}}}_{\mathrm{1},m}^n\\ & -{\displaystyle \frac{\mathrm{1}}{\mathrm{2}{T}_{\mathrm{1}}^P}}\sum _{m=-\mathrm{1}}^{\mathrm{1}}(-\mathrm{1}{)}^m{\widehat{\widehat{\mathbf{T}}}}_{\mathrm{1}\left(K\right),-m}{\widehat{\widehat{\mathbf{T}}}}_{\mathrm{1}\left(K\right),m},\end{array}\end{array}$$

where ${\widehat{\widehat{\mathbf{T}}}}_{\mathrm{1}\left(K\right),m}$ is a rank-1 irreducible spherical superoperator constructed for a spin K according to Eq. (31). A crucial property of the superoperator ${\widehat{\widehat{\mathbf{\Gamma }}}}_K$ is that it commutes with the superoperator ${\widehat{\widehat{\mathbf{F}}}}_z$:

$$\begin{array}{}\text{(36)}& \left[{\widehat{\widehat{\mathbf{F}}}}_z,{\widehat{\widehat{\mathbf{\Gamma }}}}_K\right]=\mathrm{0}.\end{array}$$

This follows from the identities

$$\begin{array}{}\text{(37)}& \begin{array}{rl}& \sum _{m=-\mathrm{1}}^{\mathrm{1}}(-\mathrm{1}{)}^m{\widehat{\widehat{\mathbf{T}}}}_{\mathrm{1},-m}^n{\widehat{\widehat{\mathbf{T}}}}_{\mathrm{1},m}^n=\left({\widehat{\widehat{\mathbf{T}}}}_{\mathrm{1}}^n\cdot {\widehat{\widehat{\mathbf{T}}}}_{\mathrm{1}}^n\right),\\ & \sum _{m=-\mathrm{1}}^{\mathrm{1}}(-\mathrm{1}{)}^m{\widehat{\widehat{\mathbf{T}}}}_{\mathrm{1}\left(K\right),-m}{\widehat{\widehat{\mathbf{T}}}}_{\mathrm{1}\left(K\right),m}=\left({\widehat{\widehat{\mathbf{T}}}}_{\mathrm{1}\left(K\right)}\cdot {\widehat{\widehat{\mathbf{T}}}}_{\mathrm{1}\left(K\right)}\right),\end{array}\end{array}$$

which represent scalar products of rank-1 spherical superoperators in Liouville space. Such scalar products are rotationally invariant and therefore commute with ${\widehat{\widehat{\mathbf{F}}}}_z$ (this can also be verified by direct evaluation of the commutators). Combining Eqs. (25) and (36), we obtain ${\widehat{\widehat{\mathbf{F}}}}_z{\widehat{\mathit{\rho }}}_K\left(t\right)=\mathrm{0}$ for any t.

Thus, relaxational dynamics driven by ${\widehat{\widehat{\mathbf{\Gamma }}}}_K$ and, therefore, the full Liouville space dynamics governed by Eq. (29) conserve the evolution within the ZQC subspace; see Fig. 2B–D. Therefore, the evolution equation takes the following form:

$$\begin{array}{}\text{(38)}& \begin{array}{rl}{\displaystyle \frac{\partial {\widehat{\mathit{\rho }}}_K\left(t\right)}{\partial t}}& =(-i{\widehat{\widehat{\mathbf{H}}}}_K^{\text{ZQC}}(t)+{\widehat{\widehat{\mathbf{\Gamma }}}}_K^{\text{ZQC}})\widehat{\mathit{\rho }}\left(t\right)\\ & ={\widehat{\widehat{\mathbf{L}}}}_K^{\text{ZQC}}\left(t\right)\widehat{\mathit{\rho }}\left(t\right),\end{array}\end{array}$$

where we introduced ${\widehat{\widehat{\mathbf{L}}}}_K^{\text{ZQC}}\left(t\right)$ as the Liouvillian superoperator projected onto the ZQC subspace. We emphasize that all components of the density matrix outside the ZQC subspace are identically zero, as shown in Fig. 2C and D.

It should be noted that a wide class of relaxation mechanisms possesses the same symmetry properties: the effective-spin reduction and conservation of a coherence order (Kowalewski and Mäler, 2006):

$$\begin{array}{}\text{(39)}& [{\widehat{\widehat{\mathbf{K}}}}^{\mathrm{2}},\widehat{\widehat{\mathbf{\Gamma }}}]=\left[{\widehat{\widehat{\mathbf{F}}}}_z,\widehat{\widehat{\mathbf{\Gamma }}}\right]=\mathrm{0},\end{array}$$

including dipole–dipole relaxation, chemical shift anisotropy (CSA) with axial symmetry, and cross-correlated relaxation between these mechanisms. Nevertheless, for the majority of SABRE experiments, the random fluctuating field (RFF) model provides a sufficiently accurate description of relaxation.

2.3 Block-diagonal decomposition of SABRE master equations

The SABRE hyperpolarization process arises from an intricate interplay between spin and chemical exchange dynamics. Polarization is transferred from a parahydrogen molecule to the substrate during its transient binding within the Ir-based polarization transfer complex, as illustrated in Fig. 3. Inside the complex, the nuclear spins of the parahydrogen-derived hydrides and the nuclear spins of the substrate become strongly coupled, leading to the formation of entangled spin states. This coherent spin evolution enables efficient transfer of spin order from parahydrogen to the substrate nuclei. Upon dissociation of the substrate from the metal complex, the accumulated spin polarization is released into the free substrate, resulting in its observable hyperpolarization of the nuclear spins.

Figure 3 Schematic representation of the SABRE polarization transfer process from parahydrogen (pH₂) to the substrate (S) mediated by the transient polarization transfer complex (C). Spin dynamics within the complex is governed by the Liouvillian superoperator ${\widehat{\widehat{\mathbf{L}}}}_C$, while spin evolution in the pool of free substrate molecules is described by ${\widehat{\widehat{\mathbf{L}}}}_S$. In our model, parahydrogen supply and exchange with the complex are assumed to be fast such that the overall exchange is effectively determined by the substrate dissociation rate constant k_d and association rate W_a. Here, L and S₁ denote the axial and equatorial ligand of the complex, respectively, and IMes is a heterocyclic carbene ligand (1,3-bis(2,4,6-trimethylphenyl)-1,3-dihydro-2H-imidazol-2-ylidene).

Thus, spin dynamics governs the polarization transfer processes occurring within the metal complex, whereas chemical dynamics determines how this polarization is transferred to and accumulated in the ensemble of free substrate molecules. Both coherent spin evolution and stochastic chemical exchange are therefore essential ingredients of the SABRE mechanism.

To rigorously consider this sophisticated SABRE dynamics, we used the SABRE master equation (Knecht et al., 2016):

$$\begin{array}{}\text{(40)}& \left\{\begin{array}{rl}{\displaystyle \frac{\mathrm{d}{\widehat{\mathit{\rho }}}_S\left(t\right)}{\mathrm{d}t}}& {\displaystyle }={\widehat{\widehat{\mathbf{L}}}}_S\left(t\right){\widehat{\mathit{\rho }}}_S\left(t\right)-W_{\mathrm{a}}{\widehat{\mathit{\rho }}}_S\left(t\right)+k_{\mathrm{d}}{\text{Tr}}_{{\mathrm{H}}_{\mathrm{2}}}\mathit{\left\{}{\widehat{\mathit{\rho }}}_C\right(t\left)\mathit{\right\}},\\ {\displaystyle \frac{\mathrm{d}{\widehat{\mathit{\rho }}}_C\left(t\right)}{\mathrm{d}t}}& {\displaystyle }={\widehat{\widehat{\mathbf{L}}}}_C\left(t\right){\widehat{\mathit{\rho }}}_C\left(t\right)-k_{\mathrm{d}}{\widehat{\mathit{\rho }}}_C\left(t\right)\\ {\displaystyle }& {\displaystyle }+W_{\mathrm{a}}\mathit{\left\{}{\widehat{\mathit{\rho }}}_S\right(t)\otimes {\widehat{\mathit{\rho }}}_{{\text{pH}}_{\mathrm{2}}}\mathit{\}},\end{array}\right.\end{array}$$

where we introduced the following quantities: ${\widehat{\mathit{\rho }}}_{S,C}\left(t\right)$ values are the density matrices of the free substrate and the SABRE complex, respectively; ${\widehat{\mathit{\rho }}}_{{\text{pH}}_{\mathrm{2}}}$ is the density matrix of parahydrogen (singlet spin state of two nuclei); ${\widehat{\widehat{\mathbf{L}}}}_{S,C}\left(t\right)=-i{\widehat{\widehat{\mathbf{H}}}}_{S,C}\left(t\right)+{\widehat{\widehat{\mathbf{\Gamma }}}}_{S,C}$ values are the Liouvillian superoperators of the substrate and the complex; k_d is the substrate dissociation rate constant; and $W_{\mathrm{a}}=k_{\mathrm{d}}\left[C\right]/\left[S\right]$ is the substrate association rate constant. Importantly, in Eq. (40) the density matrices are normalized by the following concentrations: $\text{Tr}\mathit{\left\{}{\widehat{\mathit{\rho }}}_S\right(t\left)\mathit{\right\}}=\left[S\right]/\left(\right[C]+[S\left]\right)$ and $\text{Tr}\mathit{\left\{}{\widehat{\mathit{\rho }}}_C\right(t\left)\mathit{\right\}}=\left[C\right]/\left(\right[C]+[S\left]\right)$. Equation (40) arises as the rigorous Markovian limit of the integral encounter theory, a first-principle approach to chemical kinetics in the condensed phase (Doktorov et al., 2002; Feskov et al., 2005; Ivanov et al., 2004). SABRE Eq. (40) is linear and valid in the regime of fast parahydrogen supply and fast dihydrogen exchange with the complex. Although the full SABRE kinetics is generally non-linear (Pravdivtsev and Hövener, 2019), the linear approximation is adequate for most practical purposes.

Equation (40) can be transformed to the superoperator form

$$\begin{array}{}\text{(41)}& \left\{\begin{array}{rl}{\displaystyle \frac{\mathrm{d}{\widehat{\mathit{\rho }}}_S\left(t\right)}{\mathrm{d}t}}& {\displaystyle }={\widehat{\widehat{\mathbf{A}}}}_S\left(t\right){\widehat{\mathit{\rho }}}_S\left(t\right)+k_{\mathrm{d}}{\widehat{\widehat{\mathbf{S}}}}_{{\text{TrH}}_{\mathrm{2}}}{\widehat{\mathit{\rho }}}_C\left(t\right),\\ {\displaystyle \frac{\mathrm{d}{\widehat{\mathit{\rho }}}_C\left(t\right)}{\mathrm{d}t}}& {\displaystyle }={\widehat{\widehat{\mathbf{A}}}}_C\left(t\right){\widehat{\mathit{\rho }}}_C\left(t\right)+W_{\mathrm{a}}{\widehat{\widehat{\mathbf{S}}}}_{\text{Kron}}{\widehat{\mathit{\rho }}}_S\left(t\right).\end{array}\right.\end{array}$$

Here, we introduce the following auxiliary superoperators:

$$\begin{array}{}\text{(42)}& \begin{array}{rl}{\widehat{\widehat{\mathbf{A}}}}_S\left(t\right)& ={\widehat{\widehat{\mathbf{L}}}}_S\left(t\right)-W_{\mathrm{a}}{\widehat{\widehat{\mathbf{1}}}}_S,\\ {\widehat{\widehat{\mathbf{A}}}}_C\left(t\right)& ={\widehat{\widehat{\mathbf{L}}}}_C\left(t\right)-k_{\mathrm{d}}{\widehat{\widehat{\mathbf{1}}}}_C,\end{array}\end{array}$$

where ${\widehat{\widehat{\mathbf{1}}}}_{S,C}$ values are the identity superoperators acting on the density matrices according to ${\widehat{\widehat{\mathbf{1}}}}_S{\widehat{\mathit{\rho }}}_S\left(t\right)={\widehat{\mathit{\rho }}}_S\left(t\right)$ and ${\widehat{\widehat{\mathbf{1}}}}_C{\widehat{\mathit{\rho }}}_C\left(t\right)={\widehat{\mathit{\rho }}}_C\left(t\right)$, ${\widehat{\widehat{\mathbf{S}}}}_{{\text{TrH}}_{\mathrm{2}}}{\widehat{\mathit{\rho }}}_C\left(t\right)={\text{Tr}}_{{\mathrm{H}}_{\mathrm{2}}}\mathit{\left\{}{\widehat{\mathit{\rho }}}_C\right(t\left)\mathit{\right\}}$ is the partial trace superoperator over the states of the hydrides in the complex, and ${\widehat{\widehat{\mathbf{S}}}}_{\text{Kron}}{\widehat{\mathit{\rho }}}_S\left(t\right)=\mathit{\left\{}{\widehat{\mathit{\rho }}}_S\right(t)\otimes {\widehat{\mathit{\rho }}}_{{\text{pH}}_{\mathrm{2}}}\mathit{\}}$ is the Kronecker product superoperator.

The initial state of the complex and the free substrate is typically assumed to be completely non-polarized:

$$\begin{array}{}\text{(43)}& \begin{array}{rl}{\widehat{\mathit{\rho }}}_S\left(\mathrm{0}\right)& ={\displaystyle \frac{\left[S\right]}{\left[C\right]+\left[S\right]}}{\displaystyle \frac{{\widehat{\mathrm{1}}}_S}{\text{dim}\left({\mathcal{H}}_S\right)}},\\ {\widehat{\mathit{\rho }}}_C\left(\mathrm{0}\right)& ={\displaystyle \frac{\left[C\right]}{\left[C\right]+\left[S\right]}}{\displaystyle \frac{{\widehat{\mathrm{1}}}_C}{\text{dim}\left({\mathcal{H}}_C\right)}},\end{array}\end{array}$$

where ${\widehat{\mathrm{1}}}_{S,C}$ values are the identity operators, and dim(ℋ_S,C) values are the dimensions of the Hilbert spaces of the free substrate and the complex, respectively.

2.3.1 The group of magnetically equivalent nuclei

In what follows, we demonstrate how block-diagonal decomposition can be employed to reduce the dimensionality of Eq. (41).

We note that the initial density matrices in Eq. (43) commute with any operator in the Hilbert space, as they are proportional to the identity operator. Consequently, if the substrate contains a group of magnetically equivalent nuclei G_P, it can be treated as a single effective nucleus with total spin K; see Eqs. (4) and (5). This reduction is justified because Eq. (34) holds for both the free substrate and the SABRE complex. In this case, where the magnetic equivalence is conserved in the free substrate and in the complex, the chemical exchange superoperators do not mix subspaces with different values of K.

As a consequence, the first step in the dimensionality reduction is to solve Eq. (41) independently for each value of the total spin K of the group of magnetically equivalent nuclei. To avoid unnecessary notational complexity, the index K will be omitted in what follows, with the understanding that the density matrices ${\widehat{\mathit{\rho }}}_S\left(t\right)$ and ${\widehat{\mathit{\rho }}}_C\left(t\right)$ correspond to a fixed value of K.

2.3.2 The z projection of total spin

The next step is to exploit the symmetry associated with the z projection of the total spin of the substrate, ${\widehat{\widehat{\mathbf{F}}}}_{Sz}$, and of the complex, ${\widehat{\widehat{\mathbf{F}}}}_{Cz}$. In particular, the following relations hold:

$$\begin{array}{}\text{(44)}& \left\{\begin{array}{rl}{\displaystyle }{\widehat{\widehat{\mathbf{F}}}}_{Sz}{\widehat{\mathit{\rho }}}_S\left(\mathrm{0}\right)& {\displaystyle }=\mathrm{0},\\ {\displaystyle }\left[{\widehat{\widehat{\mathbf{F}}}}_{Sz},{\widehat{\widehat{\mathbf{A}}}}_S\right]& {\displaystyle }=\mathrm{0}.\end{array}\right.\left\{\begin{array}{rl}{\displaystyle }{\widehat{\widehat{\mathbf{F}}}}_{Cz}{\widehat{\mathit{\rho }}}_C\left(\mathrm{0}\right)& {\displaystyle }=\mathrm{0},\\ {\displaystyle }\left[{\widehat{\widehat{\mathbf{F}}}}_{Cz},{\widehat{\widehat{\mathbf{A}}}}_C\right]& {\displaystyle }=\mathrm{0}.\end{array}\right.\end{array}$$

The last two identities follow from the fact that the relaxation and coherent dynamics conserve the coherence order:

$$\begin{array}{}\text{(45)}& \begin{array}{rl}& \left[{\widehat{\widehat{\mathbf{F}}}}_{Sz},{\widehat{\widehat{\mathbf{L}}}}_S\right]=-i\underset{\mathrm{0}}{\underbrace{\left[{\widehat{\widehat{\mathbf{F}}}}_{Sz},{\widehat{\widehat{\mathbf{H}}}}_S\right]}}+\underset{\mathrm{0}}{\underbrace{\left[{\widehat{\widehat{\mathbf{F}}}}_{Sz},{\widehat{\widehat{\mathbf{\Gamma }}}}_S\right]}}=\mathrm{0},\\ & \left[{\widehat{\widehat{\mathbf{F}}}}_{Cz},{\widehat{\widehat{\mathbf{L}}}}_C\right]=-i\underset{\mathrm{0}}{\underbrace{\left[{\widehat{\widehat{\mathbf{F}}}}_{Cz},{\widehat{\widehat{\mathbf{H}}}}_C\right]}}+\underset{\mathrm{0}}{\underbrace{\left[{\widehat{\widehat{\mathbf{F}}}}_{Cz},{\widehat{\widehat{\mathbf{\Gamma }}}}_C\right]}}=\mathrm{0}\end{array}\end{array}$$

and the definition of the auxiliary superoperators; see Eq. (42) (${\widehat{\widehat{\mathbf{1}}}}_{S,C}$ commutes with any superoperator). As a result, the first terms in Eq. (41), governed by the auxiliary superoperators ${\widehat{\widehat{\mathbf{A}}}}_{S,C}$, conserve the dynamics within the zero-quantum coherence (ZQC) subspaces of the Liouville spaces of the complex and the free substrate. These subspaces are denoted as ZQC_C and ZQC_S, respectively.

Additionally, in SABRE the situation is significantly complicated by the presence of reversible chemical exchange. Special attention must therefore be paid to the exchange superoperators ${\widehat{\widehat{\mathbf{S}}}}_{{\text{TrH}}_{\mathrm{2}}}$ and ${\widehat{\widehat{\mathbf{S}}}}_{\text{Kron}}$, introduced in Eq. (41), since they explicitly couple the Liouville spaces of the complex and the free substrate. In Appendices B and C, we demonstrate that these chemical exchange superoperators conserve the zero coherence order. Consequently, ${\widehat{\widehat{\mathbf{S}}}}_{\text{Kron}}:{\text{ZQC}}_S\to {\text{ZQC}}_C,{\widehat{\widehat{\mathbf{S}}}}_{{\text{TrH}}_{\mathrm{2}}}:{\text{ZQC}}_C\to {\text{ZQC}}_S$, and the following relations hold:

$$\begin{array}{}\text{(46)}& {\widehat{\widehat{\mathbf{F}}}}_{Sz}{\widehat{\mathit{\rho }}}_S\left(t\right)=\mathrm{0}\text{ and }{\widehat{\widehat{\mathbf{F}}}}_{Cz}{\widehat{\mathit{\rho }}}_C\left(t\right)=\mathrm{0}\text{ for any }t.\end{array}$$

Thus, the dynamics can be rigorously restricted to the ZQC_S and ZQC_C subspaces, as illustrated in Fig. 4. Importantly, all other coherence order subspaces can be omitted without affecting the description of SABRE spin dynamics, since the density matrices are identically zero outside the ZQC subspaces (Fig. 4).

Figure 4 Block-diagonal decomposition for SABRE master equations. Only zero-quantum coherence (ZQC) subspaces of the free substrate (ZQC_S) and the complex (ZQC_C) are relevant in the SABRE matrix. Chemical exchange superoperators ${\widehat{\widehat{\mathbf{S}}}}_{{\text{TrH}}_{\mathrm{2}}}$ and ${\widehat{\widehat{\mathbf{S}}}}_{\text{Kron}}$ couple these ZQC subspaces. Outside the ZQC subspaces, the density matrices of the free substrate ${\widehat{\mathit{\rho }}}_S$ and the complex ${\widehat{\mathit{\rho }}}_C$ are zero. The density matrices are represented in Liouville space as vectors.

Thus, the dimensionality of the problem is significantly reduced by considering that Eq. (41) is projected onto the ZQC subspaces instead of the full Liouville spaces:

$$\begin{array}{}\text{(47)}& \left\{\begin{array}{rl}{\displaystyle \frac{\mathrm{d}{\widehat{\mathit{\rho }}}_S\left(t\right)}{\mathrm{d}t}}& {\displaystyle }={\widehat{\widehat{\mathbf{A}}}}_S^{\text{ZQC}}\left(t\right){\widehat{\mathit{\rho }}}_S\left(t\right)+k_{\mathrm{d}}{\widehat{\widehat{\mathbf{S}}}}_{{\text{TrH}}_{\mathrm{2}}}^{\text{ZQC}}{\widehat{\mathit{\rho }}}_C\left(t\right),\\ {\displaystyle \frac{\mathrm{d}{\widehat{\mathit{\rho }}}_C\left(t\right)}{\mathrm{d}t}}& {\displaystyle }={\widehat{\widehat{\mathbf{A}}}}_C^{\text{ZQC}}\left(t\right){\widehat{\mathit{\rho }}}_C\left(t\right)+W_{\mathrm{a}}{\widehat{\widehat{\mathbf{S}}}}_{\text{Kron}}^{\text{ZQC}}{\widehat{\mathit{\rho }}}_S\left(t\right),\end{array}\right.\end{array}$$

where ${\widehat{\widehat{\mathbf{A}}}}_S^{\text{ZQC}}\left(t\right)$ and ${\widehat{\widehat{\mathbf{A}}}}_C^{\text{ZQC}}\left(t\right)$ are the superoperators ${\widehat{\widehat{\mathbf{A}}}}_S\left(t\right)$ and ${\widehat{\widehat{\mathbf{A}}}}_C\left(t\right)$ reduced onto the ZQC_S and ZQC_C subspaces, respectively, and ${\widehat{\widehat{\mathbf{S}}}}_{{\text{TrH}}_{\mathrm{2}}}^{\text{ZQC}}$ and ${\widehat{\widehat{\mathbf{S}}}}_{\text{Kron}}^{\text{ZQC}}$ are the exchange superoperators calculated between the corresponding ZQC blocks. Importantly, these equations provide a rigorous, complete description of the full SABRE dynamics without any additional approximations while fully exploiting the block-diagonal structure of the density matrices to reduce computational complexity. In our notation, the ZQC subscript for the density matrices in Eq. (47) is omitted, as they are identically zero outside the zero-quantum coherence subspaces.

2.3.3 The dimension of the SABRE matrix

In this work, we introduce three different SABRE matrices for comparison: the full SABRE matrix (without any reduction), the K-SABRE matrix (with effective-spin reduction for the magnetically equivalent nuclei), and the ZQC SABRE matrix (with both effective-spin and ZQC reductions). The dimensions of the matrices are simply the sum of the dimensions of the Liouville spaces for the free substrate and the complex:

$$\begin{array}{}\text{(48)}& \begin{array}{rl}& \text{dim}\left(\text{Full SABRE}\right)={\mathrm{4}}^N+{\mathrm{4}}^{N+\mathrm{2}},\\ & \text{dim}\left(\text{K-SABRE}\right)=(\mathrm{2}K+\mathrm{1}{)}^{\mathrm{2}}({\mathrm{4}}^{N-P}+{\mathrm{4}}^{N-P+\mathrm{2}}),\end{array}\end{array}$$

where N is the number of substrate spins, N+2 is the number of spins in the complex (see Fig. 4), N−P is the number of magnetically non-equivalent nuclear spins in the substrate, $N-P+\mathrm{2}$ is the number of magnetically non-equivalent nuclear spins in the complex, and K is the total spin of the group G_P of magnetically equivalent nuclei. In turn, the dimension of the ZQC SABRE matrix of Eq. (47) is given by the sum of dimensions of the ZQC subspaces:

$$\begin{array}{}\text{(49)}& \begin{array}{rl}& \text{dim}\left(\text{ZQC SABRE}\right)=\\ & \sum _{F_z=-(N-P)/\mathrm{2}-K}^{(N-P)/\mathrm{2}+K}{\left\{\sum _{m_k=-K}^K\left(\genfrac{}{}{0ex}{}{N-P}{F_z+(N-P)/\mathrm{2}-m_k}\right)\right\}}^{\mathrm{2}}\\ & +\sum _{F_z=-(N+\mathrm{2}-P)/\mathrm{2}-K}^{(N+\mathrm{2}-P)/\mathrm{2}+K}{\left\{\sum _{m_k=-K}^K\left(\genfrac{}{}{0ex}{}{N+\mathrm{2}-P}{F_z+(N+\mathrm{2}-P)/\mathrm{2}-m_k}\right)\right\}}^{\mathrm{2}},\end{array}\end{array}$$

where we used Eq. (28). The dimension of the ZQC-reduced SABRE matrix is significantly smaller than the dimension of the K-SABRE and the full SABRE matrices. The asymptotic behavior at $\sqrt{N-P}\gg K$ gives the following result (see Eq. 48 and Appendix D):

$$\begin{array}{}\text{(50)}& \begin{array}{rl}& \text{dim}\left(\text{Full SABRE}\right)\sim {\mathrm{4}}^{N+\mathrm{2}},\\ & \text{dim}\left(\text{K-SABRE}\right)\sim (\mathrm{2}K+\mathrm{1}{)}^{\mathrm{2}}{\mathrm{4}}^{N-P+\mathrm{2}},\\ & \text{dim}\left(\text{ZQC SABRE}\right)\sim {\displaystyle \frac{(\mathrm{2}K+\mathrm{1}{)}^{\mathrm{2}}}{\sqrt{\mathit{\pi }(N-P)}}}{\mathrm{4}}^{N-P+\mathrm{2}}.\end{array}\end{array}$$

Thus, the K-SABRE matrix only uses reduction based on the effective-spin treatment, which allows us to reduce the dimension by the factor $\frac{{\mathrm{4}}^P}{(\mathrm{2}K+\mathrm{1}{)}^{\mathrm{2}}}$. The following application of the ZQC reduction allows us to achieve better asymptotics for large spin systems since the dimension of the matrix is reduced by a greater factor of $\frac{\text{dim}\left(\text{Full SABRE}\right)}{\text{dim}\left(\text{ZQC SABRE}\right)}=\frac{\sqrt{\mathit{\pi }(N-P)}{\mathrm{4}}^P}{(\mathrm{2}K+\mathrm{1}{)}^{\mathrm{2}}}$ ≫ 1.

2.4 On the validity of the zero-quantum coherence reduction

Our framework with the ZQC reduction is valid if the following conditions are met.

1. The Hamiltonian

   The Hamiltonian of the spin system $\widehat{H}\left(t\right)$ is assumed to commute with the z projection of the total spin ${\widehat{F}}_z$. This assumption holds as long as no transverse (x or y) radiofrequency (RF) pulses are applied. Consequently, the proposed framework is particularly valuable for low and ultralow fields, where nontrivial spin dynamics occurs “spontaneously”, i.e., in the absence of RF pulses. Specialized state-space reduction techniques are also employed in high-field NMR simulations, notably in the Spinach software framework, which enables large-scale simulations of spin dynamics under arbitrary RF pulse sequences (Kuprov et al., 2007; Kuprov, 2018; Hogben et al., 2011).

2. The relaxation superoperator

   The relaxation superoperator $\widehat{\widehat{\mathbf{\Gamma }}}$ is assumed to commute with ${\widehat{\widehat{\mathbf{F}}}}_z$. This approximation is valid for a variety of relevant NMR relaxation mechanisms in the extreme narrowing regime, including dipolar relaxation, chemical shift anisotropy, and their cross-correlations. A possible exception is chemical shift anisotropy in the absence of axial symmetry, which represents a rather uncommon case.

3. The initial condition

   The initial density matrices of the free substrate, ${\widehat{\mathit{\rho }}}_S\left(\mathrm{0}\right)$, and the SABRE complex, ${\widehat{\mathit{\rho }}}_C\left(\mathrm{0}\right)$, are assumed to commute with their respective z projection of total spin operators; i.e., $[{\widehat{F}}_{Sz},{\widehat{\mathit{\rho }}}_S(\mathrm{0}\left)\right]=\mathrm{0}$ and $[{\widehat{F}}_{Cz},{\widehat{\mathit{\rho }}}_C(\mathrm{0}\left)\right]=\mathrm{0}$. This condition is satisfied for the standard SABRE initial condition given by Eq. (43).

It should be noted that, if only relaxation and coherent evolution are considered without chemical exchange, the dynamics can be reduced onto subspaces with fixed magnetic quantum numbers F_C and F_S without invoking the full ZQC subspaces. Strictly speaking, Eq. (46) ensures that the density matrices ${\widehat{\mathit{\rho }}}_C\left(t\right)$ and ${\widehat{\mathit{\rho }}}_S\left(t\right)$ commute with ${\widehat{F}}_{Cz}$ and ${\widehat{F}}_{Sz}$, implying their block-diagonal structure in F_C and F_S, respectively. However, the chemical exchange superoperator ${\widehat{\widehat{\mathbf{S}}}}_{{\text{TrH}}_{\mathrm{2}}}$, which describes the decay of correlations between the substrate and hydrides in the complex, couples the subspaces of the free substrate and the complex with different quantum numbers. Consequently, although the density matrices remain block diagonal in magnetic quantum numbers F_C and F_S, the SABRE equations themselves are not block diagonal in F_C and F_S, meaning that they cannot be reduced to a single block with fixed F_C and F_S. Therefore, in SABRE, the full ZQC subspaces of the substrate and the complex must be considered.

3 Methods

3.1 The algorithm for numerical simulations

This section summarizes the theoretical framework and outlines the complete algorithm for numerical simulation.

1. Identify permutation-equivalent spin groups. Determine the group G_P of magnetically equivalent spin-$\mathrm{1}/\mathrm{2}$ nuclei within the substrate molecule. This typically includes protons in chemically equivalent groups, such as $-{\mathrm{CH}}_{\mathrm{2}}-$ or −CH₃. If the molecule contains several distinct sets of magnetically equivalent nuclei, each is treated separately as groups $G_{P_{\mathrm{1}}},G_{P_{\mathrm{2}}},\mathrm{\dots }$, containing $P_{\mathrm{1}},P_{\mathrm{2}},\mathrm{\dots }$ nuclei, respectively.

2. Replace each equivalent group with an effective spin. Each group $G_{P_i}$ of magnetically equivalent spins is replaced by a single effective spin K_i, whose quantum number can take values in the range

   $$\begin{array}{}\text{(51)}& K_i\in [K_{i,\text{min}},P_i/\mathrm{2}],\end{array}$$

   where $K_{i,\text{min}}=\mathrm{0}$ if P_i is even, and $K_{i,\text{min}}=\mathrm{1}/\mathrm{2}$ if P_i is odd.

3. Construct the effective spin system. The original system of N nuclei is thereby reduced to a system consisting of the following:

   • the $(N-P_{\mathrm{1}}-P_{\mathrm{2}}-\mathrm{\dots })$ non-equivalent spin-$\mathrm{1}/\mathrm{2}$ nuclei,

   • a single effective spin K₁ (representing group $G_{P_{\mathrm{1}}}$),

   • a single effective spin K₂ (representing group $G_{P_{\mathrm{2}}}$), and so on.

4. Build the Hamiltonian superoperators for the complex and free substrate. The Hamiltonian operator for the effective-spin system reads as follows (see Eq. 10):

   $$\begin{array}{}\text{(52)}& \begin{array}{rl}& {\widehat{H}}_{K_{\mathrm{1}},K_{\mathrm{2}},\mathrm{\dots }}\left(t\right)=\\ & -B_{\mathrm{0}}\left(t\right)\left(\sum _{l=\mathrm{1}}^{N-P_{\mathrm{1}}-P_{\mathrm{2}}-\mathrm{\dots }}{\mathit{\gamma }}^l{\widehat{I}}_z^l+{\mathit{\gamma }}^{P_{\mathrm{1}}}{\widehat{K}}_{\mathrm{1}z}+{\mathit{\gamma }}^{P_{\mathrm{2}}}{\widehat{K}}_{\mathrm{2}z}+\mathrm{\dots }\right)\\ & +\mathrm{2}\mathit{\pi }\sum _{l<m}^{N-P_{\mathrm{1}}-P_{\mathrm{2}}-\mathrm{\dots }}{J}^{lm}({\widehat{\mathbf{I}}}^l\cdot {\widehat{\mathbf{I}}}^m)\\ & +\mathrm{2}\mathit{\pi }\sum _{l=\mathrm{1}}^{N-P_{\mathrm{1}}-P_{\mathrm{2}}-\mathrm{\dots }}\left(J^{l{G}_{P_{\mathrm{1}}}}\right({\widehat{\mathbf{I}}}^l\cdot {\widehat{\mathbf{K}}}_{\mathrm{1}})+J^{l{G}_{P_{\mathrm{2}}}}({\widehat{\mathbf{I}}}^l\cdot {\widehat{\mathbf{K}}}_{\mathrm{2}})+\mathrm{\dots })\\ & +\mathrm{2}\mathit{\pi }\left(J^{G_{P_{\mathrm{1}}}{G}_{P_{\mathrm{2}}}}\right({\widehat{\mathbf{K}}}_{\mathrm{1}}\cdot {\widehat{\mathbf{K}}}_{\mathrm{2}})+J^{G_{P_{\mathrm{1}}}{G}_{P_{\mathrm{3}}}}({\widehat{\mathbf{K}}}_{\mathrm{1}}\cdot {\widehat{\mathbf{K}}}_{\mathrm{3}})+\mathrm{\dots }),\end{array}\end{array}$$

   where the sums run over all magnetically non-equivalent spin-$\mathrm{1}/\mathrm{2}$ nuclei, and the remaining terms include the scalar J couplings within and between the different spin groups of the equivalent nuclei. The Hamiltonian is then mapped to Liouville space as the commutation superoperator; see Eq. (19):

   $$\begin{array}{}\text{(53)}& \begin{array}{rl}\widehat{\widehat{\mathbf{H}}}& {}_{K_{\mathrm{1}},K_{\mathrm{2}},\mathrm{\dots }}\left(t\right)=\\ & {\widehat{H}}_{K_{\mathrm{1}},K_{\mathrm{2}},\mathrm{\dots }}\left(t\right)\otimes \widehat{\mathrm{1}}-\widehat{\mathrm{1}}\otimes {\widehat{H}}_{K_{\mathrm{1}},K_{\mathrm{2}},\mathrm{\dots }}(t{)}^{\top },\end{array}\end{array}$$

   where $\widehat{\mathrm{1}}$ is the identity operator, and ⊤ denotes the matrix transpose. The same construction is applied separately to the SABRE complex and the free substrate, using their respective spin systems and parameters.

5. Build the relaxation superoperators for the complex and free substrate. The relaxation superoperator for the effective-spin system is given by (see Eq. 35)

   $$\begin{array}{}\text{(54)}& \begin{array}{rl}{\widehat{\widehat{\mathbf{\Gamma }}}}_{K_{\mathrm{1}},K_{\mathrm{2}},\mathrm{\dots }}& =-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\sum _{n=\mathrm{1}}^{N-P_{\mathrm{1}}-P_{\mathrm{2}}-\mathrm{\dots }}{\displaystyle \frac{\mathrm{1}}{T_{\mathrm{1}}^n}}\sum _{m=-\mathrm{1}}^{\mathrm{1}}(-\mathrm{1}{)}^m{\widehat{\widehat{\mathbf{T}}}}_{\mathrm{1},-m}^n{\widehat{\widehat{\mathbf{T}}}}_{\mathrm{1},m}^n\\ & -{\displaystyle \frac{\mathrm{1}}{\mathrm{2}{T}_{\mathrm{1}}^{P_{\mathrm{1}}}}}\sum _{m=-\mathrm{1}}^{\mathrm{1}}(-\mathrm{1}{)}^m{\widehat{\widehat{\mathbf{T}}}}_{\mathrm{1}\left(K_{\mathrm{1}}\right),-m}{\widehat{\widehat{\mathbf{T}}}}_{\mathrm{1}\left(K_{\mathrm{1}}\right),m}\\ & -\mathrm{\dots },\end{array}\end{array}$$

   where the first sum runs over all $(N-P_{\mathrm{1}}-P_{\mathrm{2}}-\mathrm{\dots })$ magnetically non-equivalent spin-$\mathrm{1}/\mathrm{2}$ nuclei, and other terms represent a contribution for a given group of magnetically equivalent spins with longitudinal relaxation time T₁. The irreducible spherical superoperators are given by

   $$\begin{array}{}\text{(55)}& \begin{array}{rl}& {\widehat{\widehat{\mathbf{T}}}}_{\mathrm{1},m}^n={\widehat{T}}_{\mathrm{1},m}^n\otimes \widehat{\mathrm{1}}-\widehat{\mathrm{1}}\otimes {\widehat{T}}_{\mathrm{1},m}^{n\top },\\ & {\widehat{\widehat{\mathbf{T}}}}_{\mathrm{1}\left(K_i\right),m}={\widehat{T}}_{\mathrm{1}\left(K_i\right),m}\otimes \widehat{\mathrm{1}}-\widehat{\mathrm{1}}\otimes {\widehat{T}}_{\mathrm{1}\left(K_i\right),m}^{\top }.\end{array}\end{array}$$

   The corresponding irreducible spherical tensor operators ${\widehat{T}}_{\mathrm{1},m}^n$ (for a spin-$\mathrm{1}/\mathrm{2}$ nucleus) and ${\widehat{T}}_{\mathrm{1}\left(K_i\right),m}$ (for an effective spin K_i) are constructed according to the standard definition in Eq. (31), using the appropriate spin quantum number ($\mathrm{1}/\mathrm{2}$ or K_i, respectively). This construction is performed independently for the SABRE complex and the free substrate, using their respective spin systems and relaxation times T₁.

6. Construct the projected Liouvillians in the zero-quantum coherence subspaces. The total Liouvillian superoperators for the free substrate (S) and the SABRE complex (C) are constructed by combining the respective Hamiltonian and relaxation superoperators from steps 4 and 5:

   $$\begin{array}{}\text{(56)}& {\widehat{\widehat{\mathbf{L}}}}_{S,C}\left(t\right)=-i{\widehat{\widehat{\mathbf{H}}}}_{S,C}\left(t\right)+{\widehat{\widehat{\mathbf{\Gamma }}}}_{S,C},\end{array}$$

   where ${\widehat{\widehat{\mathbf{H}}}}_{S,C}\left(t\right)$ is defined by Eq. (53) and ${\widehat{\widehat{\mathbf{\Gamma }}}}_{S,C}$ by Eq. (54). Following Eq. (41), we then form the auxiliary superoperators ${\widehat{\widehat{\mathbf{A}}}}_{S,C}\left(t\right)$ for each system. Finally, the superoperators ${\widehat{\widehat{\mathbf{A}}}}_{S,C}\left(t\right)$, see Eq. (42), are projected onto the zero-quantum coherence (ZQC) subspaces of the free substrate (ZQC_S) and the complex (ZQC_C), respectively:

   $$\begin{array}{}\text{(57)}& {\widehat{\widehat{\mathbf{A}}}}_S^{\text{ZQC}}\left(t\right)=\left({\widehat{s}}_i\left|{\widehat{\widehat{\mathbf{A}}}}_S\left(t\right)\right|{\widehat{s}}_j\right),{\widehat{s}}_i,{\widehat{s}}_j\in {\text{ZQC}}_S,\end{array}$$

   where the brackets denote the scalar product in Liouville space, ${\widehat{s}}_{i,j}$ values are the basis operators in the ZQC_S subspace (see Eq. 26), and indices i,j run over all possible values in the ZQC_S subspace. In the same manner, we introduce

   $$\begin{array}{}\text{(58)}& {\widehat{\widehat{\mathbf{A}}}}_C^{\text{ZQC}}\left(t\right)=\left({\widehat{c}}_i\left|{\widehat{\widehat{\mathbf{A}}}}_C\left(t\right)\right|{\widehat{c}}_j\right),{\widehat{c}}_i,{\widehat{c}}_j\in {\text{ZQC}}_C,\end{array}$$

   with ${\widehat{c}}_{i,j}$ being the basis operators in the ZQC_C subspace.

7. Construct the chemical exchange superoperators. The chemical exchange superoperators ${\widehat{\widehat{\mathbf{S}}}}_{{\text{TrH}}_{\mathrm{2}}}^{\text{ZQC}}$ and ${\widehat{\widehat{\mathbf{S}}}}_{\text{Kron}}^{\text{ZQC}}$, which describe the reversible binding dynamics between the free substrate and the complex, are constructed following the formalism detailed in the supporting information of Knecht et al. (2016). These exchange superoperators are then projected to act between the previously defined zero-quantum coherence subspaces of the free substrate (ZQC_S) and the complex (ZQC_C). This yields the following final projected exchange superoperators:

   $$\begin{array}{}\text{(59)}& \begin{array}{rl}& {\widehat{\widehat{\mathbf{S}}}}_{{\text{TrH}}_{\mathrm{2}}}^{\text{ZQC}}=\left({\widehat{s}}_i\left|{\widehat{\widehat{\mathbf{S}}}}_{{\text{TrH}}_{\mathrm{2}}}\right|{\widehat{c}}_j\right),\\ & {\widehat{\widehat{\mathbf{S}}}}_{\text{Kron}}^{\text{ZQC}}=\left({\widehat{c}}_j\left|{\widehat{\widehat{\mathbf{S}}}}_{\text{Kron}}\right|{\widehat{s}}_i\right),\\ & \text{with }{\widehat{s}}_i\in {\text{ZQC}}_S,{\widehat{c}}_j\in {\text{ZQC}}_C,\end{array}\end{array}$$

   where the basis operators of the ZQC subspaces, ${\widehat{s}}_i,{\widehat{c}}_j$, are introduced in Eqs. (57) and (58). These superoperators correctly couple the reduced-dimensionality state representations of the two chemical species.

8. Build and solve the SABRE master equation in the ZQC subspaces. The time-dependent master equation, Eq. (47), is then solved by numerical integration, using the constructed ZQC SABRE matrix and appropriate initial conditions; see Eq. (43). The solution yields the time-evolving density matrices of the free substrate and the complex, confined to their respective ZQC subspaces: ${\widehat{\mathit{\rho }}}_S\left(t\right)\in {\text{ZQC}}_S$, ${\widehat{\mathit{\rho }}}_C\left(t\right)\in {\text{ZQC}}_C$. A key feature of this representation is that all matrix elements of the density matrices outside these subspaces are identically zero.

9. Calculate the physical observables. The primary quantities of interest are the expectation values of spin operators for the free substrate. These are calculated using the reduced density matrix ${\widehat{\mathit{\rho }}}_S\left(t\right)$ obtained in step 8. More precisely, for a given set of effective high-spin quantum numbers $(K_{\mathrm{1}},K_{\mathrm{2}},\mathrm{\dots })$, we operate with the corresponding density matrix ${\widehat{\mathit{\rho }}}_{(S,K_{\mathrm{1}},K_{\mathrm{2}},\mathrm{\dots })}\left(t\right)$. The expectation value of an observable operator ${\widehat{O}}_{(K_{\mathrm{1}},K_{\mathrm{2}},\mathrm{\dots })}$ (e.g., magnetization along a specific axis) is computed as

   $$\begin{array}{}\text{(60)}& \langle {\widehat{O}}_{K_{\mathrm{1}},K_{\mathrm{2}},\mathrm{\dots }}\rangle \left(t\right)={\displaystyle \frac{\text{Tr}\mathit{\{}{\widehat{O}}_{K_{\mathrm{1}},K_{\mathrm{2}},\mathrm{\dots }}\cdot {\widehat{\mathit{\rho }}}_{S,K_{\mathrm{1}},K_{\mathrm{2}},\mathrm{\dots }}(t\left)\mathit{\right\}}}{\text{Tr}\left\{{\widehat{\mathit{\rho }}}_{S,K_{\mathrm{1}},K_{\mathrm{2}},\mathrm{\dots }}\left(t\right)\right\}}},\end{array}$$

   where the trace is taken over the Hilbert space of the free substrate's spin system. The effective-spin description represents a statistical mixture of different total spin states for each group of magnetically equivalent nuclei. Therefore, the final, physically observable $\widehat{O}$ expectation value is obtained by summing the results from Eq. (60) over all possible combinations $(K_{\mathrm{1}},K_{\mathrm{2}},\mathrm{\dots })$, weighted by their respective statistical weights:

   $$\begin{array}{}\text{(61)}& \langle \widehat{O}\rangle \left(t\right)=\sum _{K_{\mathrm{1}},K_{\mathrm{2}},\mathrm{\dots }}{g}_{K_{\mathrm{1}},K_{\mathrm{2}},\mathrm{\dots }}\langle {\widehat{O}}_{K_{\mathrm{1}},K_{\mathrm{2}},\mathrm{\dots }}\rangle \left(t\right),\end{array}$$

   where the total statistical weights are defined as $g_{K_{\mathrm{1}},K_{\mathrm{2}},\mathrm{\dots }}=g_{K_{\mathrm{1}}}\cdot g_{K_{\mathrm{2}}}\mathrm{\cdots }$. The individual factors $g_{K_i}$ for a group of P_i magnetically equivalent spin-$\mathrm{1}/\mathrm{2}$ nuclei are defined by Eqs. (3) and (7).

3.2 Computational resources

All numerical simulations and matrix manipulations were performed using the in-house Python code. The calculations were executed on a standard desktop workstation with the following specifications: CPU 2 × Intel Xeon E5-2620 v3 at 2.40 GHz (12 cores, 24 threads) and RAM 256 GB DDR4.

The most memory-intensive step of our calculations was the storage of the ZQC SABRE matrix in the compressed sparse row (CSR) format for a system of N=14 nuclear spins, represented by a matrix of dimension $\sim \mathrm{3.3}\times {\mathrm{10}}^{\mathrm{6}}\times \mathrm{3.3}\times {\mathrm{10}}^{\mathrm{6}}$. This required approximately 7 GB of RAM. The total computation time per magnetic field with this spin system was about 3.5 h. The total computation time of the ZULF SABRE spectrum was about 100 h.

All scalar J-coupling constants and longitudinal T₁-relaxation times used in the simulations are provided in Appendix E. In all cases, we considered the complex to be a single substrate molecule and two hydride ligands.

4 Results and discussion

This section validates the ZQC subspace reduction method for computationally tractable, small spin systems by demonstrating that it yields results identical to those of the full, exact calculation while providing a dramatic speedup (∼30–50 fold). This rigorous benchmark establishes the accuracy of the reduction. Furthermore, for a system of fully non-equivalent spins (representing the worst-case scenario for the reduction), the dimension of the Liouville space matrices is reduced by a factor of $\mathcal{O}\left(\sqrt{N}\right)$, and the computational time is reduced by a factor of 𝒪(N) compared to full Liouville space simulations, where N is the number of substrate spins, underscoring the scalability of the proposed approach.

For the large spin systems, where full simulations are computationally prohibitive, the validated ZQC method becomes an essential tool. Its performance gain unlocks the modeling of multi-spin systems, as illustrated in Fig. 5: panel A shows the zero-to-ultralow field (ZULF) SABRE experimental protocol, while panel B presents the two studied substrates: the small system ([¹⁵N,¹³C₂]acetonitrile, 6 spins) and the large system ([¹⁵N,¹³C₄]butyronitrile, 12 spins).

Figure 5 Schematic of the ZULF SABRE experiment and molecular systems. (A) The two-stage zero-to-ultralow field (ZULF) SABRE protocol: (1) polarization generation at a magnetic field B₀∼1 µT for a duration t_pol and (2) ZULF NMR signal detection during the interval t_ZULF at a residual near-zero field B_ZULF∼0.01 µT. The ZULF NMR signal is obtained via Fourier transform of ZULF free induction decay (ZULF FID). Parahydrogen bubbling is maintained continuously throughout both stages. (B) Chemical structures of the isotopically labeled substrates used in the simulations: [¹⁵N,¹³C₂]acetonitrile (6 spins) and [¹⁵N,¹³C₄]butyronitrile (12 spins). The numbers in parentheses indicate the total number of spin-$\mathrm{1}/\mathrm{2}$ nuclei in each system.

Using this approach, we address two central aspects of ZULF SABRE (see Fig. 5A): (1) calculating the magnetic field dependence of hyperpolarization to identify the optimal polarization field B₀ and (2) simulating the resulting ZULF NMR spectra for direct comparison with the experiment.

4.1 Validation of the ZQC reduction method: [¹⁵N,¹³C₂]acetonitrile

The [¹⁵N,¹³C₂]acetonitrile molecule constitutes a relatively small spin system (6 spins in the substrate, 8 spins in the SABRE complex), allowing for a direct numerical benchmark between the exact, unreduced calculation and the ZQC reduction approach.

First, we simulated the first stage of the ZULF SABRE experiment (Fig. 5A), the magnetic field dependence of the polarization in [¹⁵N,¹³C₂]acetonitrile, by varying the polarization field B₀ over a range of microtesla values. We specifically calculated the net polarization of each nucleus, corresponding to the observable operator ${\widehat{I}}_z$. The resulting magnetic field dependences for the different nuclei of acetonitrile are presented in Fig. 6: 1-¹⁵N (Fig. 6A), 2-¹³C (Fig. 6B), 3-¹³C (Fig. 6C), and the methyl protons (Fig. 6D).

For each nucleus, the ZQC-reduced calculation was performed for two fixed effective spin states of the three equivalent protons in the −CH₃ group: $K=\mathrm{1}/\mathrm{2}$ (top row) and $K=\mathrm{3}/\mathrm{2}$ (middle row). The final, physically observable dependence (bottom row) was obtained by statistically averaging these two results with equal weights, $g_{\mathrm{1}/\mathrm{2}}=g_{\mathrm{3}/\mathrm{2}}=\mathrm{1}/\mathrm{2}$, as prescribed by Eqs. (3), (7), and (61). For validation, the results of the exact, full SABRE matrix calculation (without any reduction) are included in the bottom row of Fig. 6 as dashed lines. They are virtually indistinguishable from the results obtained with the ZQC reduction, demonstrating a negligible relative numerical residual of the order of 10⁻⁵.

Figure 6 Magnetic field dependence of the hyperpolarization in free [¹⁵N,¹³C₂]acetonitrile computed via the ZQC reduction method. The panels correspond to (A) 1-¹⁵N, (B) 2-¹³C, (C) 3-¹³C, and (D) ¹H nuclei. For each nucleus, results are shown for three representations of the methyl (CH₃) group: top row – fixed effective spin $K=\mathrm{1}/\mathrm{2}$ for the three equivalent protons, middle row – fixed effective spin $K=\mathrm{3}/\mathrm{2}$, and bottom row – weighted statistical average over $K=\mathrm{1}/\mathrm{2}$ and $K=\mathrm{3}/\mathrm{2}$ states with weights $g_{\mathrm{1}/\mathrm{2}}=g_{\mathrm{3}/\mathrm{2}}=\mathrm{1}/\mathrm{2}$. In the bottom-row panel, the dashed white lines show the result of the exact, full Liouville space calculation (without any reduction) for the complete spin system, while the solid lines correspond to the ZQC-reduced calculation. The simulation parameters are the dissociation rate constant k_d=10 s⁻¹, catalyst to substrate ratio $\left[C\right]/\left[S\right]=\mathrm{0.1}$, and polarization time t_pol=1 s.

Crucially, a direct benchmark comparison of the different computational approaches is summarized in Table 1. The ZQC reduction method achieves a 34-fold reduction in computation time compared to the full calculation. For reference, Table 1 also includes the intermediate case using only effective-spin reduction (without ZQC projection), which provides a smaller, 6-fold speedup. Notably, the observed speedup roughly scales with the ratio of the matrix dimensions, confirming that the ZQC projection effectively reduces the computational cost in proportion to the reduced size of the Liouville space block.

Table 1 Computation times per magnetic field point for three approaches: (1) the full SABRE matrix, (2) the SABRE matrix with effective-spin reduction of the equivalent nuclei in the methyl or methylene groups (denoted as the K-SABRE matrix), and (3) the ZQC-reduced SABRE matrix (denoted by boldface).

^a The maximum dimension is the largest matrix size for the system, corresponding to the highest possible effective spins of the equivalent groups: $K=\mathrm{3}/\mathrm{2}$ (CH₃) for acetonitrile and K₁=1, K₂=1, and $K_{\mathrm{3}}=\mathrm{3}/\mathrm{2}$ (two CH₂ and one CH₃) for butyronitrile.
^b For the 8-spin system, the ZQC SABRE matrix reduces the computation time ∼ 34-fold.
^c For the 14-spin system, the ZQC SABRE matrix reduces the computation time ∼ 86-fold.

The obtained magnetic field dependences exhibit broad profiles (extending up to ∼3 µT) and a complex structure featuring multiple extrema and sign changes. This nontrivial behavior stems from the coherent interplay of several types of nuclei (¹⁵N, ¹³C, ¹H) within the molecule, coupled through the network of scalar interactions. Motivated by this nontrivial magnetic field dependence, we extended the simulations to a broader range (up to 10 µT) and examined the effect of the dissociation rate constant k_d. The resulting two-dimensional maps of hyperpolarization are shown in Fig. 7 for each nucleus: 1-¹⁵N (Fig. 7A), 2-¹³C (Fig. 7B), 3-¹³C (Fig. 7C), and the methyl protons (Fig. 7D). Critically, the ZQC reduction renders such two-dimensional scans computationally feasible, requiring only ∼1 h despite the complexity of the simulated system, which includes 8 magnetic nuclei in the SABRE complex. Thus, the maximum polarization for all nuclei occurs at approximately 0.5 µT across the entire range of dissociation rates k_d. Beyond this common optimum, several distinct features emerge: (1) the width of the 1-¹⁵N polarization profile increases systematically with k_d (see Fig. 7A), (2) the 2-¹³C nucleus exhibits two extrema – a primary positive maximum near 0.5 µT and a secondary negative extremum at higher fields (2–3 µT) (see Fig. 7B), and (3) the 3-¹³C nucleus displays an unusually complex and broad profile with multiple extrema persisting up to 10 µT. This pronounced multi-peak structure is non-typical in SABRE and directly reflects the intricate network of heteronuclear couplings within the multi-spin molecule.

Figure 7 Two-dimensional dependence of the hyperpolarization in free [¹⁵N,¹³C₂]acetonitrile on the dissociation rate constant k_d and the polarization transfer field B₀, computed via the ZQC reduction method. The panels correspond to (A) 1-¹⁵N, (B) 2-¹³C, (C) 3-¹³C, and (D) ¹H nuclei. The simulation parameters are catalyst to substrate ratio $\left[C\right]/\left[S\right]=\mathrm{0.1}$ and polarization time t_pol=1 s.

As the next step, we simulated the complete ZULF SABRE spectrum using the optimized two-stage protocol: (1) pre-polarization at the optimal field B₀=0.5 µT (identified above), followed by (2) ZULF FID signal detection at a near-zero field B_ZULF=0.01 µT. The detected ZULF signal is proportional to the total z magnetization of all nuclei, represented by the operator

$$\begin{array}{}\text{(62)}& {\widehat{O}}_{\text{ZULF}}=\sum _l{\mathit{\gamma }}^l{\widehat{I}}_z^l,\end{array}$$

where the sum runs over all spins in the substrate. A detailed procedure for simulating such ZULF NMR spectra is described in Stern and Sheberstov (2023). In brief, we computed the time-dependent expectation value $\langle {\widehat{O}}_{\text{ZULF}}\rangle \left(t\right)$, see Eqs. (60) and (61), and then obtained the frequency-domain spectrum via a standard Fourier transform.

The resulting simulated ZULF SABRE spectrum of free [¹⁵N,¹³C₂]acetonitrile is shown in Fig. 8. The solid line shows the spectrum obtained with the ZQC reduction method, while the dashed line corresponds to the exact, full calculation without any reduction. The two curves are visually indistinguishable, with a negligible relative numerical difference of the order of 10⁻⁵, providing a final validation of the accuracy of the ZQC approach. The computational speedup in this case is ∼54-fold with the ZQC reduction.

Figure 8 Simulated ZULF SABRE spectrum of free [¹⁵N,¹³C₂]acetonitrile. The right part of the spectrum is multiplied by a factor of 8. The dashed line shows the result of the exact, full-space calculation (without ZQC or effective spin reduction) for the complete 6-spin system. The solid line corresponds to the spectrum obtained via the ZQC reduction method with effective spin treatment. The vertical dashed lines mark the characteristic spectral frequencies of acetonitrile. The high-frequency domain at ¹J_CH = 136 Hz and 2×¹J_CH = 272 Hz and the low-frequency domain at ¹J_CN, ${\mathit{\nu }}_{\mathrm{1}}=\frac{\mathrm{5}}{\mathrm{8}}\left(J_{\mathrm{CC}}+\mathrm{3}{\times }^{\mathrm{2}}{J}_{\mathrm{CH}}\right)$ = 19.4 Hz, ${\mathit{\nu }}_{\mathrm{2}}=\frac{\mathrm{3}}{\mathrm{4}}\left(J_{\mathrm{CC}}+{}^{\mathrm{2}}{J}_{\mathrm{CH}}\right)$ = 36.8 Hz, and ${\mathit{\nu }}_{\mathrm{3}}=\frac{\mathrm{3}}{\mathrm{8}}\left(J_{\mathrm{CC}}-\mathrm{5}{\times }^{\mathrm{2}}{J}_{\mathrm{CH}}\right)$ = 38.7 Hz. The simulation parameters are polarization field B₀ = 0.5 µT, detection field B_ZULF = 0.01 µT, polarization time t_pol = 10 s, acquisition time t_ZULF = 5 s, dissociation rate constant k_d = 100 s⁻¹, and catalyst to substrate ratio $\left[C\right]/\left[S\right]$ = 0.1. The computation times were 1415 s for the full matrix and 26 s for the ZQC-reduced matrix, demonstrating a ∼ 54-fold speedup.

The simulated ZULF NMR spectrum exhibits characteristic fingerprints of different scalar couplings in acetonitrile. A fully analytical description of the spectrum is nontrivial; however, qualitative insight can be obtained by considering the ¹³C–¹³CH₃ moiety as an effective B(XA₃) spin system, as analyzed previously (Barskiy et al., 2025). Within this simplified picture, prominent high-frequency features are expected near the couplings ¹J_CH=136 Hz and 2×¹J_CH = 272 Hz (see Appendix E for the J couplings). In addition, several low-frequency transitions are predicted in the first-order perturbation theory at ${\mathit{\nu }}_{\mathrm{1}}=\frac{\mathrm{5}}{\mathrm{8}}\left(J_{\mathrm{CC}}+\mathrm{3}{\times }^{\mathrm{2}}{J}_{\mathrm{CH}}\right)$ = 19.4 Hz, ${\mathit{\nu }}_{\mathrm{2}}=\frac{\mathrm{3}}{\mathrm{4}}\left(J_{\mathrm{CC}}+{}^{\mathrm{2}}{J}_{\mathrm{CH}}\right)$ = 36.8 Hz, and ${\mathit{\nu }}_{\mathrm{3}}=\frac{\mathrm{3}}{\mathrm{8}}\left(J_{\mathrm{CC}}-\mathrm{5}{\times }^{\mathrm{2}}{J}_{\mathrm{CH}}\right)$ = 38.7 Hz. The inclusion of the ¹⁵N nucleus, which is coupled to ¹³C via the scalar coupling ¹J_CN = −18 Hz, calls for a more detailed analysis beyond the reduced-spin model considered above. Although additional transitions arise upon inclusion of ¹⁵N, the overall spectral structure remains in qualitative agreement with the simplified description. We note that a ZULF NMR spectrum of [¹⁵N,¹³C₂]acetonitrile with thermal pre-polarization at 1.8 T has been reported previously in Ledbetter et al. (2011). In contrast to the thermal case, the relative line intensities in the SABRE-enhanced ZULF spectrum differ substantially, reflecting the non-Boltzmann spin orders generated via SABRE.

4.2 Simulation of a large spin ensemble: [¹⁵N,¹³C₄butyronitrile]

[¹⁵N,¹³C₄]butyronitrile provides a representative example of a multi-spin system comprising 12 magnetically active nuclei so that the corresponding SABRE complex contains a total of 14 spins. Direct simulations of such systems rapidly become computationally demanding on a standard desktop computer. As shown in Table 1, the computation time required for a single magnetic field point is approximately 300 h for the full Liouville space treatment and about 24 h when only the effective-spin reduction is applied, resulting in impractically long simulation times. The use of the ZQC reduction decreases the dimensionality of the problem by a factor of ∼88 (Table 1), reducing the computation time to approximately 3.5 h per field point and enabling systematic simulations of such multi-spin SABRE systems. Importantly, in our treatment we considered the $-{\mathrm{CH}}_{\mathrm{2}}-$ and −CH₃ groups of butyronitrile as the corresponding effective spins and then averaged the results with different statistical weights shown in Table 2.

Table 2 Effective-spin configurations of the $-{\mathrm{CH}}_{\mathrm{2}}-$ and −CH₃ groups of [¹⁵N–¹³C₄]butyronitrile and their corresponding statistical weights.

We first calculated the magnetic field dependences of the ¹⁵N and ¹³C polarization during the pre-polarization stage of the ZULF SABRE experiment, as shown in Fig. 9A and B, respectively. The simulated results capture the main qualitative features: (1) there is a very broad dependence of the 1-¹⁵N polarization on the pre-polarization field B₀, extending up to ∼2 µT, reflecting the complexity of the underlying multi-spin system, and (2) the polarization levels of the different ¹³C sites decrease along the chain. Importantly, our calculations can incorporate the actual sweeping profile of the magnetic field ramp, which is essential for capturing adiabatic effects, although doing so will increase the computational time. It should also be noted that the 3-¹³C and 4-¹³C sites are close to chemical equivalence; i.e. the difference in their chemical shifts is smaller than their J-coupling constants with other nuclei. Accordingly, we report their total magnetization (3-¹³C + 4-¹³C) as the measured observable.

Figure 9 ZULF SABRE of free [¹⁵N,¹³C₄]butyronitrile. (A, B) Simulated magnetic field dependences of the ¹⁵N (A) and ¹³C (B) polarization. (C) Simulated ZULF NMR spectrum with the ZQC reduction obtained after pre-polarization at B₀ = 0.5 µT. (Top) Full ZULF NMR spectrum, in which only the low-frequency features are visible. (Bottom) The same spectrum scaled by a factor of 25, revealing characteristic high-frequency features associated with the scalar couplings ¹J_CH = 136 Hz, $\mathrm{3}/\mathrm{2}{\times }^{\mathrm{1}}{J}_{\mathrm{CH}}$ = 204 Hz, and 2×¹J_CH = 272 Hz, indicated by dashed lines. The total computation time of the ZULF NMR spectrum was ∼ 100 h. The simulation parameters are detection field B_ZULF = 0.01 µT, polarization time t_pol = 10 s, acquisition time t_ZULF = 5 s, dissociation rate constant k_d = 50 s⁻¹, and catalyst to substrate ratio $\left[C\right]/\left[S\right]$ = 0.37.

The simulated ZULF NMR spectrum is shown in Fig. 9C. The most intense signal appears in the low-frequency domain, centered around ∼10 Hz, and broadens as the number of coupled magnetic nuclei increases. In addition to this dominant low-frequency feature, characteristic fingerprints of the $-{\mathrm{CH}}_{\mathrm{2}}-$ and −CH₃ groups are visible: frequencies ¹J_CH=136 Hz and 2×¹J_CH = 272 Hz from the −CH₃ group and $\mathrm{3}/\mathrm{2}{\times }^{\mathrm{1}}{J}_{\mathrm{CH}}$ = 204 Hz from the $-{\mathrm{CH}}_{\mathrm{2}}-$ group. These signals are broadened due to splitting by the J couplings present in the multi-spin system.

4.3 Benchmarking and scalability analysis

In this section, we analyze the scalability of the proposed ZQC reduction approach in comparison to full Liouville space simulations as a function of the number of spins in the substrate molecule. To demonstrate its efficiency, we consider the worst-case scenario in which the substrate contains only N non-equivalent nuclei (the presence of magnetically equivalent nuclei would only further exponentially improve the scalability; see Eq. 50). The scaling benchmark as a function of N is shown in Fig. 10.

Figure 10 Scaling benchmark of full and ZQC-reduced Liouville space matrices for a substrate molecule with N non-equivalent spins. (A) Dimensionality of the full SABRE (open squares) and ZQC-reduced SABRE (filled circles) Liouville space matrices as a function of N. The solid lines represent the asymptotic scaling for $\sqrt{N}\gg \mathrm{1}$. Note that the SABRE complex contains N+2 spins. (B) Dimensionality reduction achieved with the ZQC approach. The predicted asymptotic scaling is $\sqrt{\mathit{\pi }N}$. (C) Computation time per magnetic field B₀ for the full SABRE (squares) and ZQC-reduced SABRE (circles) Liouville space matrices as a function of N. (D) Computational speedup achieved with the ZQC approach, which scales as 𝒪(N).

We first examine the dimensionality of the matrices, shown in Fig. 10A. For full Liouville space simulations, the dimensionality is given exactly by $\text{dim(Full SABRE)}={\mathrm{4}}^N+{\mathrm{4}}^{N+\mathrm{2}}$ (the SABRE complex contains N+2 spins). For the ZQC-reduced matrices, a compact asymptotic expression valid for $\sqrt{N}\gg \mathrm{1}$ is given by Eq. (50): $\text{dim(ZQC SABRE)}\sim {\mathrm{4}}^{N+\mathrm{2}}/\sqrt{\mathit{\pi }N}$. As shown in Fig. 10A, this asymptotic approximation (solid line) almost perfectly matches the dimensionality of the ZQC-reduced SABRE matrix. Although the dimensionality still grows exponentially with N, it is reduced by a factor of $\sqrt{\mathit{\pi }N}$ compared to the full Liouville space (see Fig. 10B).

Similarly, we investigated the dependence of the computation time per single magnetic field B₀ on the number of spins N, as shown in Fig. 10C. For small $N=\mathrm{1},\mathrm{2}$, the computation time of the ZQC approach is nearly identical for both matrices. For larger N, however, it scales linearly, i.e., as 𝒪(N) (Fig. 10D). This demonstrates that the reduction in computation time scales quadratically with the reduction in matrix dimension.

5 Conclusions

In this work, we have shown that the Hamiltonian, relaxation, and chemical exchange superoperators in SABRE systems under zero-to-ultralow-field conditions all possess well-defined symmetry with respect to the z projection of the total spin, where the z axis is defined by the external ultralow magnetic field. This approach is valid across different coupling regimes encountered in ZULF NMR, from the J-coupling-dominated to the Zeeman-dominated limits, making the approach general and broadly applicable. Importantly, this symmetry ensures that chemical exchange conserves the coherence order, allowing the dynamics to be rigorously restricted to the zero-quantum coherence (ZQC) subspace. As only operators within this subspace contribute to polarization transfer in SABRE, focusing on ZQC drastically reduces the effective dimensionality of the Liouville space while retaining a fully rigorous description of the spin dynamics.

As a first validation, we applied this approach to [¹⁵N,¹³C₂]acetonitrile (8 spins in the SABRE complex), reproducing full Liouville space results with excellent agreement while accelerating computations by a factor of ∼30 for magnetic field dependence calculations and ∼54 for ZULF NMR spectra. To demonstrate the scalability of the method, we further applied it to [¹⁵N,¹³C₄]butyronitrile (14 spins in the SABRE complex), a system that would be practically impossible to simulate with full Liouville space calculations within a reasonable time frame. Using the ZQC reduction, these simulations became feasible, providing insight into both field dependence of polarization and ZULF NMR spectra. These results confirm the benchmarking analysis for a system of N non-equivalent substrate spins: the matrix dimension is reduced roughly by a factor of $\sqrt{\mathit{\pi }N}$, while the computation time is accelerated roughly by a factor of N.

Overall, the present framework provides a predictive and scalable tool for analyzing complex SABRE spin dynamics in the absence of transverse oscillating fields, which is particularly valuable for capturing the rich structure of field-dependent polarization curves and the intricate ZULF NMR spectra of multi-spin molecules. This approach therefore offers a practical pathway for designing experiments and optimizing hyperpolarization strategies in large (up to 15 spins) and chemically diverse spin systems.

Appendix A: The coherence orders

The operator ${\widehat{Q}}_p$ is said to have a coherence order p if it satisfies the following relation with the operator of the z projection of total spin ${\widehat{F}}_z$:

$$\begin{array}{}\text{(A1)}& [{\widehat{F}}_z,{\widehat{Q}}_p]=p{\widehat{Q}}_p.\end{array}$$

In terms of Liouville space, where the operator ${\widehat{F}}_z$ is mapped to the superoperator ${\widehat{\widehat{\mathbf{F}}}}_z$, Eq. (A1) transforms to the eigenvalue problem:

$$\begin{array}{}\text{(A2)}& [{\widehat{F}}_z,{\widehat{Q}}_p]={\widehat{\widehat{\mathbf{F}}}}_z{\widehat{Q}}_p=p{\widehat{Q}}_p.\end{array}$$

Thus, ${\widehat{Q}}_p$ is an eigenoperator of ${\widehat{\widehat{\mathbf{F}}}}_z$ with eigenvalue p.

As a demonstrative example, we consider two spin-$\mathrm{1}/\mathrm{2}$ particles and define ${\widehat{F}}_z={\widehat{I}}_{\mathrm{1}z}+{\widehat{I}}_{\mathrm{2}z}$. As an example, we determine the coherence order of the operator $\widehat{Q}=|\mathit{\alpha }\mathit{\alpha }\rangle \langle \mathit{\beta }\mathit{\beta }|$, where |α〉 and |β〉 denote the standard Zeeman “spin-up” and “spin-down” states, respectively. In this case, we substitute $\widehat{Q}$ into Eq. (A2):

$$\begin{array}{}\text{(A3)}& \begin{array}{rl}{\widehat{\widehat{\mathbf{F}}}}_z|\mathit{\alpha }\mathit{\alpha }\rangle \langle \mathit{\beta }\mathit{\beta }|& =[{\widehat{F}}_z,|\mathit{\alpha }\mathit{\alpha }\rangle \langle \mathit{\beta }\mathit{\beta }\left|\right]=\left[\right({\widehat{I}}_{\mathrm{1}z}+{\widehat{I}}_{\mathrm{2}z}),|\mathit{\alpha }\mathit{\alpha }\rangle \langle \mathit{\beta }\mathit{\beta }\left|\right]\\ & =({\widehat{I}}_{\mathrm{1}z}+{\widehat{I}}_{\mathrm{2}z})|\mathit{\alpha }\mathit{\alpha }\rangle \langle \mathit{\beta }\mathit{\beta }|-|\mathit{\alpha }\mathit{\alpha }\rangle \langle \mathit{\beta }\mathit{\beta }|({\widehat{I}}_{\mathrm{1}z}+{\widehat{I}}_{\mathrm{2}z})\\ & =(\frac{\mathrm{1}}{\mathrm{2}}+\frac{\mathrm{1}}{\mathrm{2}})|\mathit{\alpha }\mathit{\alpha }\rangle \langle \mathit{\beta }\mathit{\beta }|-|\mathit{\alpha }\mathit{\alpha }\rangle \langle \mathit{\beta }\mathit{\beta }|(-\frac{\mathrm{1}}{\mathrm{2}}-\frac{\mathrm{1}}{\mathrm{2}})\\ & =+\mathrm{2}|\mathit{\alpha }\mathit{\alpha }\rangle \langle \mathit{\beta }\mathit{\beta }|.\end{array}\end{array}$$

Therefore, $|\mathit{\alpha }\mathit{\alpha }\rangle \langle \mathit{\beta }\mathit{\beta }|$ is a +Double-Quantum Coherence (+DQC) operator. In the same manner, one can find that the populations are the zero-quantum coherence (ZQC) operators:

$$\begin{array}{}\text{(A4)}& \begin{array}{rl}{Q}_{\mathrm{0}}^{\left(\mathrm{1}\right)}& =|\mathit{\alpha }\mathit{\alpha }\rangle \langle \mathit{\alpha }\mathit{\alpha }|,Q_{\mathrm{0}}^{\left(\mathrm{2}\right)}=|\mathit{\alpha }\mathit{\beta }\rangle \langle \mathit{\alpha }\mathit{\beta }|,\\ Q_{\mathrm{0}}^{\left(\mathrm{3}\right)}& =|\mathit{\beta }\mathit{\alpha }\rangle \langle \mathit{\beta }\mathit{\alpha }|,Q_{\mathrm{0}}^{\left(\mathrm{4}\right)}=|\mathit{\beta }\mathit{\beta }\rangle \langle \mathit{\beta }\mathit{\beta }|.\end{array}\end{array}$$

Additionally, two more ZQC operators can be found:

$$\begin{array}{}\text{(A5)}& Q_{\mathrm{0}}^{\left(\mathrm{5}\right)}=|\mathit{\alpha }\mathit{\beta }\rangle \langle \mathit{\beta }\mathit{\alpha }|,Q_{\mathrm{0}}^{\left(\mathrm{6}\right)}=|\mathit{\beta }\mathit{\alpha }\rangle \langle \mathit{\alpha }\mathit{\beta }|.\end{array}$$

Generally, for an arbitrary spin system, a ZQC operator is constructed as

$$\begin{array}{}\text{(A6)}& {\widehat{Q}}_{\mathrm{0}}=|F_z,\mathit{\xi }\rangle \langle F_z,\mathit{\lambda }|,\end{array}$$

where F_z, which runs over all possible values for the considered spin system, and ξ,λ denote different spin states with identical values of F_z. This construction ensures that ${\widehat{\widehat{\mathbf{F}}}}_z{\widehat{Q}}_{\mathrm{0}}=\mathrm{0}$.

For example, the ZQC operators in Eqs. (A4) and (A5) can be expressed using the notations introduced in Eq. (A6):

$$\begin{array}{}\text{(A7)}& \begin{array}{rl}& Q_{\mathrm{0}}^{\left(\mathrm{1}\right)}=\left|\mathrm{1}\right(\mathrm{1}/\mathrm{2},\mathrm{1}/\mathrm{2})\rangle \langle \mathrm{1}(\mathrm{1}/\mathrm{2},\mathrm{1}/\mathrm{2}\left)\right|,\\ & Q_{\mathrm{0}}^{\left(\mathrm{2}\right)}=\left|\mathrm{0}\right(\mathrm{1}/\mathrm{2},-\mathrm{1}/\mathrm{2})\rangle \langle \mathrm{0}(\mathrm{1}/\mathrm{2},-\mathrm{1}/\mathrm{2}\left)\right|,\\ & Q_{\mathrm{0}}^{\left(\mathrm{3}\right)}=\left|\mathrm{0}\right(-\mathrm{1}/\mathrm{2},\mathrm{1}/\mathrm{2})\rangle \langle \mathrm{0}(-\mathrm{1}/\mathrm{2},\mathrm{1}/\mathrm{2}\left)\right|,\\ & Q_{\mathrm{0}}^{\left(\mathrm{4}\right)}=|-\mathrm{1}(-\mathrm{1}/\mathrm{2},-\mathrm{1}/\mathrm{2})\rangle \langle -\mathrm{1}(-\mathrm{1}/\mathrm{2},-\mathrm{1}/\mathrm{2}\left)\right|,\\ & Q_{\mathrm{0}}^{\left(\mathrm{5}\right)}=\left|\mathrm{0}\right(\mathrm{1}/\mathrm{2},-\mathrm{1}/\mathrm{2})\rangle \langle \mathrm{0}(-\mathrm{1}/\mathrm{2},\mathrm{1}/\mathrm{2}\left)\right|,\\ & Q_{\mathrm{0}}^{\left(\mathrm{6}\right)}=\left|\mathrm{0}\right(-\mathrm{1}/\mathrm{2},\mathrm{1}/\mathrm{2})\rangle \langle \mathrm{0}(\mathrm{1}/\mathrm{2},-\mathrm{1}/\mathrm{2}\left)\right|.\end{array}\end{array}$$

Here, the individual magnetic quantum numbers m_z of two $\mathrm{1}/\mathrm{2}$ particles are shown in brackets.

It is therefore straightforward to calculate the full number of independent ZQC operators. The total number $\text{dim}\left({\mathcal{H}}_{F_z}\right)$ of different spin states with a fixed F_z is given in Eq. (15). The number of ZQC operators for a fixed F_z is given by ${\left\{\text{dim}\left({\mathcal{H}}_{F_z}\right)\right\}}^{\mathrm{2}}$, as follows directly from the “ket-bra” construction of Eq. (A6). Finally, the total number of ZQC operators is obtained by summing over all possible F_z, yielding Eq. (28).

In a broader context, the role of the coherence orders in modern NMR techniques is discussed in Ivanov et al. (2023).

Appendix B: Theorem 1. Conservation of zero coherence order in SABRE chemical reaction

B1 Dissociation from the complex

$$\begin{array}{}\text{(B1)}& \text{If }{\widehat{\widehat{\mathbf{F}}}}_{Cz}{\widehat{\mathit{\rho }}}_C\left(t\right)=\mathrm{0}\text{, then }{\widehat{\widehat{\mathbf{F}}}}_{Sz}{\text{Tr}}_{{\mathrm{H}}_{\mathrm{2}}}\mathit{\left\{}{\widehat{\mathit{\rho }}}_C\right(t\left)\mathit{\right\}}=\mathrm{0}\end{array}$$

Proof

The superoperator ${\widehat{\widehat{\mathbf{F}}}}_z$ acts on the density matrix according to the following rule:

$$\begin{array}{}\text{(B2)}& {\widehat{\widehat{\mathbf{F}}}}_z\widehat{\mathit{\rho }}\left(t\right)=[{\widehat{F}}_z,\widehat{\mathit{\rho }}(t\left)\right].\end{array}$$

Consider the following equation:

$$\begin{array}{}\text{(B3)}& \begin{array}{rl}{\text{Tr}}_{{\mathrm{H}}_{\mathrm{2}}}& \left\{{\widehat{\widehat{\mathbf{F}}}}_{Cz}{\widehat{\mathit{\rho }}}_C\left(t\right)\right\}\\ & ={\text{Tr}}_{{\mathrm{H}}_{\mathrm{2}}}\left\{\left({\widehat{\widehat{\mathbf{F}}}}_{Sz}+{\widehat{\widehat{\mathbf{F}}}}_{{\mathrm{H}}_{\mathrm{2}}z}\right){\widehat{\mathit{\rho }}}_C\left(t\right)\right\}\\ & ={\text{Tr}}_{{\mathrm{H}}_{\mathrm{2}}}\left\{{\widehat{\widehat{\mathbf{F}}}}_{Sz}{\widehat{\mathit{\rho }}}_C\left(t\right)\right\}+{\text{Tr}}_{{\mathrm{H}}_{\mathrm{2}}}\left\{{\widehat{\widehat{\mathbf{F}}}}_{{\mathrm{H}}_{\mathrm{2}}z}{\widehat{\mathit{\rho }}}_C\left(t\right)\right\},\end{array}\end{array}$$

where we accounted for ${\widehat{\widehat{\mathbf{F}}}}_{Cz}={\widehat{\widehat{\mathbf{F}}}}_{Sz}+{\widehat{\widehat{\mathbf{F}}}}_{{\mathrm{H}}_{\mathrm{2}}z}$. Here, the first superoperator acts only on the substrate’s spins, whereas the second superoperator acts only on the spins of the hydrides.

At first, we prove the following relation:

$$\begin{array}{}\text{(B4)}& {\text{Tr}}_{{\mathrm{H}}_{\mathrm{2}}}\left\{{\widehat{\widehat{\mathbf{F}}}}_{{\mathrm{H}}_{\mathrm{2}}z}{\widehat{\mathit{\rho }}}_C\left(t\right)\right\}=\mathrm{0}.\end{array}$$

We first use the definition given in Eq. (B1):

$$\begin{array}{}\text{(B5)}& \begin{array}{rl}{\text{Tr}}_{{\mathrm{H}}_{\mathrm{2}}}& \left\{{\widehat{\widehat{\mathbf{F}}}}_{{\mathrm{H}}_{\mathrm{2}}z}{\widehat{\mathit{\rho }}}_C\left(t\right)\right\}\\ & ={\text{Tr}}_{{\mathrm{H}}_{\mathrm{2}}}\mathit{\left\{}\right[{\widehat{F}}_{{\mathrm{H}}_{\mathrm{2}}z},{\widehat{\mathit{\rho }}}_C\left(t\right)\left]\mathit{\right\}}\\ & ={\text{Tr}}_{{\mathrm{H}}_{\mathrm{2}}}\mathit{\left\{}{\widehat{F}}_{{\mathrm{H}}_{\mathrm{2}}z}{\widehat{\mathit{\rho }}}_C\right(t\left)\mathit{\right\}}-{\text{Tr}}_{{\mathrm{H}}_{\mathrm{2}}}\mathit{\left\{}{\widehat{\mathit{\rho }}}_C\right(t\left){\widehat{F}}_{{\mathrm{H}}_{\mathrm{2}}z}\mathit{\right\}}.\end{array}\end{array}$$

The operator ${\widehat{F}}_{{\mathrm{H}}_{\mathrm{2}}z}$ is separable in the basis set of the complex since it acts only on the states of the hydrides:

$$\begin{array}{}\text{(B6)}& \begin{array}{rl}{\widehat{F}}_{{\mathrm{H}}_{\mathrm{2}}z}& =\sum _{i,j}\mathit{\left\{}{\widehat{F}}_{{\mathrm{H}}_{\mathrm{2}}z}{\mathit{\}}}_{i,j}\right|i\rangle \langle j|\otimes \widehat{\mathrm{1}}\\ & =\sum _{i,j}\mathit{\left\{}{\widehat{F}}_{{\mathrm{H}}_{\mathrm{2}}z}{\mathit{\}}}_{i,j}\right|i\rangle \langle j|\otimes \sum _{\mathit{\xi }}|\mathit{\xi }\rangle \langle \mathit{\xi }|\\ & =\sum _{i,j,\mathit{\xi }}\mathit{\left\{}{\widehat{F}}_{{\mathrm{H}}_{\mathrm{2}}z}{\mathit{\}}}_{i,j}\right|i\rangle |\mathit{\xi }\rangle \langle j|\langle \mathit{\xi }|,\end{array}\end{array}$$

where the Latin letters denote the states of the hydrides, and the Greek letters denote the state of the substrate in the complex. Now we consider the density matrix of the complex in this product basis set:

$$\begin{array}{}\text{(B7)}& {\widehat{\mathit{\rho }}}_C\left(t\right)=\sum _{i,j,\mathit{\xi },\mathit{\mu }}\mathit{\left\{}{\widehat{\mathit{\rho }}}_C\right(t\left){\mathit{\}}}_{i\mathit{\xi },j\mathit{\mu }}\right|i\rangle |\mathit{\xi }\rangle \langle j|\langle \mathit{\mu }|.\end{array}$$

Now we consider the matrix elements of ${\widehat{F}}_{{\mathrm{H}}_{\mathrm{2}}z}{\widehat{\mathit{\rho }}}_C$ using Eqs. (B6) and (B7):

$$\begin{array}{}\text{(B8)}& \mathit{\{}{\widehat{F}}_{{\mathrm{H}}_{\mathrm{2}}z}{\widehat{\mathit{\rho }}}_C{\mathit{\}}}_{k\mathit{\alpha },l\mathit{\beta }}=\sum _j\mathit{\{}{\widehat{F}}_{{\mathrm{H}}_{\mathrm{2}}z}{\mathit{\}}}_{k,j}\mathit{\left\{}{\widehat{\mathit{\rho }}}_C\right(t){\mathit{\}}}_{j\mathit{\alpha },l\mathit{\beta }}.\end{array}$$

We then take the partial trace over the states of the hydrides in Eq. (B8):

$$\begin{array}{}\text{(B9)}& \begin{array}{rl}{\left\{{\text{Tr}}_{{\mathrm{H}}_{\mathrm{2}}}\mathit{\left\{}{\widehat{F}}_{{\mathrm{H}}_{\mathrm{2}}z}{\widehat{\mathit{\rho }}}_C\right(t\left)\mathit{\right\}}\right\}}_{\mathit{\alpha },\mathit{\beta }}& =\sum _k\mathit{\left\{}{\widehat{F}}_{{\mathrm{H}}_{\mathrm{2}}z}{\widehat{\mathit{\rho }}}_C\right(t){\mathit{\}}}_{k\mathit{\alpha },k\mathit{\beta }}\\ & =\sum _{k,j}\mathit{\left\{}{\widehat{F}}_{{\mathrm{H}}_{\mathrm{2}}z}{\mathit{\}}}_{k,j}\mathit{\right\{}{\widehat{\mathit{\rho }}}_C\left(t\right){\mathit{\}}}_{j\mathit{\alpha },k\mathit{\beta }}.\end{array}\end{array}$$

The matrix elements of ${\widehat{\mathit{\rho }}}_C\left(t\right){\widehat{F}}_{{\mathrm{H}}_{\mathrm{2}}z}$ are defined as follows:

$$\begin{array}{}\text{(B10)}& \mathit{\left\{}{\widehat{\mathit{\rho }}}_C\right(t){\widehat{F}}_{{\mathrm{H}}_{\mathrm{2}}z}{\mathit{\}}}_{k\mathit{\alpha },l\mathit{\beta }}=\sum _j\mathit{\{}{\widehat{\mathit{\rho }}}_C\left(t\right){\mathit{\}}}_{k\mathit{\alpha },j\mathit{\beta }}\mathit{\{}{\widehat{F}}_{{\mathrm{H}}_{\mathrm{2}}z}{\mathit{\}}}_{j,l}.\end{array}$$

The partial trace of Eq. (B10) leads to the following result:

$$\begin{array}{}\text{(B11)}& \begin{array}{rl}\mathit{\{}{\text{Tr}}_{{\mathrm{H}}_{\mathrm{2}}}& \mathit{\left\{}{\widehat{\mathit{\rho }}}_C\right(t\left){\widehat{F}}_{{\mathrm{H}}_{\mathrm{2}}z}\mathit{\right\}}{\mathit{\}}}_{\mathit{\alpha },\mathit{\beta }}\\ & =\sum _k\mathit{\left\{}{\widehat{\mathit{\rho }}}_C\right(t){\widehat{F}}_{{\mathrm{H}}_{\mathrm{2}}z}{\mathit{\}}}_{k\mathit{\alpha },k\mathit{\beta }}\\ & =\sum _{k,j}\mathit{\left\{}{\widehat{\mathit{\rho }}}_C\right(t\left){\mathit{\}}}_{k\mathit{\alpha },j\mathit{\beta }}\mathit{\right\{}{\widehat{F}}_{{\mathrm{H}}_{\mathrm{2}}z}{\mathit{\}}}_{j,k}\\ & {=}_{(j\leftrightarrow k)}\sum _{k,j}\mathit{\left\{}{\widehat{F}}_{{\mathrm{H}}_{\mathrm{2}}z}{\mathit{\}}}_{k,j}\mathit{\right\{}{\widehat{\mathit{\rho }}}_C\left(t\right){\mathit{\}}}_{j\mathit{\alpha },k\mathit{\beta }}.\end{array}\end{array}$$

Now we subtract Eq. (B11) from Eq. (B9), which gives Eq. (B4):

$$\begin{array}{}\text{(B12)}& \begin{array}{rl}{\text{Tr}}_{{\mathrm{H}}_{\mathrm{2}}}& \mathit{\left\{}{\widehat{\widehat{\mathbf{F}}}}_{{\mathrm{H}}_{\mathrm{2}}z}{\widehat{\mathit{\rho }}}_C\right(t\left)\mathit{\right\}}\\ & ={\text{Tr}}_{{\mathrm{H}}_{\mathrm{2}}}\mathit{\left\{}\right[{\widehat{F}}_{{\mathrm{H}}_{\mathrm{2}}z},{\widehat{\mathit{\rho }}}_C\left(t\right)\left]\mathit{\right\}}\\ & ={\text{Tr}}_{{\mathrm{H}}_{\mathrm{2}}}\mathit{\left\{}{\widehat{F}}_{{\mathrm{H}}_{\mathrm{2}}z}{\widehat{\mathit{\rho }}}_C\right(t\left)\mathit{\right\}}-{\text{Tr}}_{{\mathrm{H}}_{\mathrm{2}}}\mathit{\left\{}{\widehat{\mathit{\rho }}}_C\right(t\left){\widehat{F}}_{{\mathrm{H}}_{\mathrm{2}}z}\mathit{\right\}}\equiv \mathrm{0}.\end{array}\end{array}$$

We then apply Eq. (B4) to Eq. (B3):

$$\begin{array}{}\text{(B13)}& {\text{Tr}}_{{\mathrm{H}}_{\mathrm{2}}}\left\{{\widehat{\widehat{\mathbf{F}}}}_{Cz}{\widehat{\mathit{\rho }}}_C\left(t\right)\right\}={\text{Tr}}_{{\mathrm{H}}_{\mathrm{2}}}\left\{{\widehat{\widehat{\mathbf{F}}}}_{Sz}{\widehat{\mathit{\rho }}}_C\left(t\right)\right\}.\end{array}$$

From the statement of the theorem, we have ${\widehat{\widehat{\mathbf{F}}}}_{Cz}{\widehat{\mathit{\rho }}}_C\left(t\right)=\mathrm{0}$. Therefore, using Eq. (B13), we have

$$\begin{array}{}\text{(B14)}& \mathrm{0}={\text{Tr}}_{{\mathrm{H}}_{\mathrm{2}}}\left\{{\widehat{\widehat{\mathbf{F}}}}_{Sz}{\widehat{\mathit{\rho }}}_C\left(t\right)\right\}={\widehat{\widehat{\mathbf{F}}}}_{Sz}{\text{Tr}}_{{\mathrm{H}}_{\mathrm{2}}}\mathit{\left\{}{\widehat{\mathit{\rho }}}_C\right(t\left)\mathit{\right\}},\end{array}$$

where at the last stage we used the fact that ${\widehat{\widehat{\mathbf{F}}}}_{Sz}$ and ${\text{Tr}}_{{\mathrm{H}}_{\mathrm{2}}}$ commute as they act on different variables. Finally, Eq. (B14) proves that zero coherence order is conserved after dissociation from the complex.

B2 Association with the complex

$$\begin{array}{}\text{(B15)}& \text{If }{\widehat{\widehat{\mathbf{F}}}}_{Sz}{\widehat{\mathit{\rho }}}_S\left(t\right)=\mathrm{0}\text{, then }{\widehat{\widehat{\mathbf{F}}}}_{Cz}\mathit{\left\{}{\widehat{\mathit{\rho }}}_S\right(t)\otimes {\widehat{\mathit{\rho }}}_{{\text{pH}}_{\mathrm{2}}}\mathit{\}}=\mathrm{0}\end{array}$$

Proof

This is trivial since for parahydrogen (the singlet nuclear spin state with zero coherence order), we have ${\widehat{\widehat{\mathbf{F}}}}_{{\mathrm{H}}_{\mathrm{2}}z}{\widehat{\mathit{\rho }}}_{{\text{pH}}_{\mathrm{2}}}=\mathrm{0}$. Therefore,

$$\begin{array}{}\text{(B16)}& \begin{array}{rl}{\widehat{\widehat{\mathbf{F}}}}_{Cz}& \mathit{\left\{}{\widehat{\mathit{\rho }}}_S\right(t)\otimes {\widehat{\mathit{\rho }}}_{{\text{pH}}_{\mathrm{2}}}\mathit{\}}\\ & =\left({\widehat{\widehat{\mathbf{F}}}}_{Sz}+{\widehat{\widehat{\mathbf{F}}}}_{{\mathrm{H}}_{\mathrm{2}}z}\right)\mathit{\left\{}{\widehat{\mathit{\rho }}}_S\right(t)\otimes {\widehat{\mathit{\rho }}}_{{\text{pH}}_{\mathrm{2}}}\mathit{\}}\\ & =\underset{\mathrm{0}}{\underbrace{{\widehat{\widehat{\mathbf{F}}}}_{Sz}{\widehat{\mathit{\rho }}}_S\left(t\right)}}\otimes {\widehat{\mathit{\rho }}}_{{\text{pH}}_{\mathrm{2}}}+{\widehat{\mathit{\rho }}}_S\left(t\right)\otimes \underset{\mathrm{0}}{\underbrace{{\widehat{\widehat{\mathbf{F}}}}_{{\mathrm{H}}_{\mathrm{2}}z}{\widehat{\mathit{\rho }}}_{{\text{pH}}_{\mathrm{2}}}}}=\mathrm{0}.\end{array}\end{array}$$

Appendix C: Theorem 2. Restriction of SABRE dynamics to the zero-quantum coherence subspaces

If the following conditions are fulfilled for the Liouvillians and the initial state of the SABRE system,

$$\begin{array}{}\text{(C1)}& {\displaystyle }& {\displaystyle }\left[{\widehat{\widehat{\mathbf{F}}}}_{Cz},{\widehat{\widehat{\mathbf{L}}}}_C\left(t\right)\right]=\mathrm{0}\text{ for any }t,\text{(C2)}& {\displaystyle }& {\displaystyle }\left[{\widehat{\widehat{\mathbf{F}}}}_{Sz},{\widehat{\widehat{\mathbf{L}}}}_S\left(t\right)\right]=\mathrm{0}\text{ for any }t,\text{(C3)}& {\displaystyle }& {\displaystyle }{\widehat{\widehat{\mathbf{F}}}}_{Cz}{\widehat{\mathit{\rho }}}_C\left(\mathrm{0}\right)=\mathrm{0},\text{(C4)}& {\displaystyle }& {\displaystyle }{\widehat{\widehat{\mathbf{F}}}}_{Sz}{\widehat{\mathit{\rho }}}_S\left(\mathrm{0}\right)=\mathrm{0},\end{array}$$

then

$$\begin{array}{}\text{(C5)}& {\widehat{\widehat{\mathbf{F}}}}_{Cz}{\widehat{\mathit{\rho }}}_C\left(t\right)=\mathrm{0},{\widehat{\widehat{\mathbf{F}}}}_{Sz}{\widehat{\mathit{\rho }}}_S\left(t\right)=\mathrm{0}\text{ for any }t.\end{array}$$

Proof

The proof uses a numerical solution of Eq. (40) discretized on a time grid with a step dt:

$$\begin{array}{}\text{(C6)}& {\displaystyle }{\widehat{\mathit{\rho }}}_C\left(\mathrm{d}t\right)-{\widehat{\mathit{\rho }}}_C\left(\mathrm{0}\right)={\widehat{\widehat{\mathbf{A}}}}_C\left(\mathrm{0}\right){\widehat{\mathit{\rho }}}_C\left(\mathrm{0}\right)\mathrm{d}t+W_{\mathrm{a}}\mathit{\left\{}{\widehat{\mathit{\rho }}}_S\right(\mathrm{0})\otimes {\widehat{\mathit{\rho }}}_{{\text{pH}}_{\mathrm{2}}}\mathit{\}},\text{(C7)}& {\displaystyle }{\widehat{\mathit{\rho }}}_S\left(\mathrm{d}t\right)-{\widehat{\mathit{\rho }}}_S\left(\mathrm{0}\right)={\widehat{\widehat{\mathbf{A}}}}_S\left(\mathrm{0}\right){\widehat{\mathit{\rho }}}_S\left(\mathrm{0}\right)\mathrm{d}t+k_{\mathrm{d}}{\text{Tr}}_{{\mathrm{H}}_{\mathrm{2}}}\mathit{\left\{}{\widehat{\mathit{\rho }}}_C\right(\mathrm{0}\left)\mathit{\right\}}.\end{array}$$

We now apply ${\widehat{\widehat{\mathbf{F}}}}_{Sz}$ for both sides of Eq. (C7):

$$\begin{array}{}\text{(C8)}& \begin{array}{rl}{\widehat{\widehat{\mathbf{F}}}}_{Sz}& {\widehat{\mathit{\rho }}}_S\left(\mathrm{d}t\right)\\ & =\underset{\mathrm{0}}{\underbrace{{\widehat{\widehat{\mathbf{F}}}}_{Sz}{\widehat{\mathit{\rho }}}_S\left(\mathrm{0}\right)}}+\underset{\mathrm{0}}{\underbrace{{\widehat{\widehat{\mathbf{F}}}}_{Sz}{\widehat{\widehat{\mathbf{A}}}}_S\left(\mathrm{0}\right){\widehat{\mathit{\rho }}}_S\left(\mathrm{0}\right)}}\mathrm{d}t\\ & +k_{\mathrm{d}}\underset{\mathrm{0}}{\underbrace{{\widehat{\widehat{\mathbf{F}}}}_{Sz}{\text{Tr}}_{{\mathrm{H}}_{\mathrm{2}}}\mathit{\left\{}{\widehat{\mathit{\rho }}}_C\right(\mathrm{0}\left)\mathit{\right\}}}}\mathrm{d}t=\mathrm{0},\end{array}\end{array}$$

where each term is zero:

1. ${\widehat{\widehat{\mathbf{F}}}}_{Sz}{\widehat{\mathit{\rho }}}_S\left(\mathrm{0}\right)=[{\widehat{F}}_{Sz},{\widehat{\mathit{\rho }}}_S(\mathrm{0}\left)\right]=\mathrm{0}$, as it follows from the initial condition of Eq. (43). According to this, ${\widehat{\mathit{\rho }}}_S\left(\mathrm{0}\right)$ is proportional to the identity matrix that commutes with any operator. Consequently, Eq. (C4) is indeed valid.

2. ${\widehat{\widehat{\mathbf{F}}}}_{Sz}{\widehat{\widehat{\mathbf{A}}}}_S\left(\mathrm{0}\right){\widehat{\mathit{\rho }}}_S\left(\mathrm{0}\right)={\widehat{\widehat{\mathbf{F}}}}_{Sz}\left\{{\widehat{\widehat{\mathbf{L}}}}_S\left(\mathrm{0}\right)-W_{\mathrm{a}}{\widehat{\widehat{\mathbf{1}}}}_S\right\}{\widehat{\mathit{\rho }}}_S\left(\mathrm{0}\right)=\left\{{\widehat{\widehat{\mathbf{L}}}}_S\left(\mathrm{0}\right)-W_{\mathrm{a}}{\widehat{\widehat{\mathbf{1}}}}_S\right\}\underset{\mathrm{0}}{\underbrace{{\widehat{\widehat{\mathbf{F}}}}_{Sz}{\widehat{\mathit{\rho }}}_S\left(\mathrm{0}\right)}}=\mathrm{0}$ Here we first used Eq. (C2) and then Eq. (C4).

3. ${\widehat{\widehat{\mathbf{F}}}}_{Sz}{\text{Tr}}_{{\mathrm{H}}_{\mathrm{2}}}\mathit{\left\{}{\widehat{\mathit{\rho }}}_C\right(\mathrm{0}\left)\mathit{\right\}}=\mathrm{0}$ by Eq. (B1) since for the initial state of the complex Eq. (C3), ${\widehat{\widehat{\mathbf{F}}}}_{Cz}{\widehat{\mathit{\rho }}}_C\left(\mathrm{0}\right)=\mathrm{0}$.

Similarly, applying ${\widehat{\widehat{\mathbf{F}}}}_{Cz}$ for both sides of Eq. (C6) and using Eqs. (C1), (C3), (C4), and (B15), we obtain

$$\begin{array}{}\text{(C9)}& {\widehat{\widehat{\mathbf{F}}}}_{Cz}{\widehat{\mathit{\rho }}}_C\left(\mathrm{d}t\right)=\mathrm{0}.\end{array}$$

The following proof is based on induction by dt.

Induction hypothesis

$$\begin{array}{}\text{(C10)}& {\widehat{\widehat{\mathbf{F}}}}_{Sz}{\widehat{\mathit{\rho }}}_S\left(t_{n-\mathrm{1}}\right)=\mathrm{0},{\widehat{\widehat{\mathbf{F}}}}_{Cz}{\widehat{\mathit{\rho }}}_C\left(t_{n-\mathrm{1}}\right)=\mathrm{0}\end{array}$$

Here $t_k=k\cdot \mathrm{d}t$ is the kth node of the time grid. From this, we prove that

$$\begin{array}{}\text{(C11)}& {\widehat{\widehat{\mathbf{F}}}}_{Sz}{\widehat{\mathit{\rho }}}_S\left(t_n\right)=\mathrm{0},{\widehat{\widehat{\mathbf{F}}}}_{Cz}{\widehat{\mathit{\rho }}}_C\left(t_n\right)=\mathrm{0}.\end{array}$$

Inductive step

From Eqs. (40), we have

$$\begin{array}{}\text{(C12)}& \begin{array}{rl}{\widehat{\mathit{\rho }}}_C\left(t_n\right)-{\widehat{\mathit{\rho }}}_C\left(t_{n-\mathrm{1}}\right)& ={\widehat{\widehat{\mathbf{A}}}}_C\left(t_{n-\mathrm{1}}\right){\widehat{\mathit{\rho }}}_C\left(t_{n-\mathrm{1}}\right)\mathrm{d}t\\ & +W_{\mathrm{a}}\mathit{\left\{}{\widehat{\mathit{\rho }}}_S\right(t_{n-\mathrm{1}})\otimes {\widehat{\mathit{\rho }}}_{{\text{pH}}_{\mathrm{2}}}\mathit{\}}\mathrm{d}t.\end{array}\end{array}$$

Applying ${\widehat{\widehat{\mathbf{F}}}}_{Cz}$ for both sides of Eq. (C12), we obtain the following:

$$\begin{array}{}\text{(C13)}& \begin{array}{rl}{\widehat{\widehat{\mathbf{F}}}}_{Cz}& {\widehat{\mathit{\rho }}}_C\left(t_n\right)\\ & ={\widehat{\widehat{\mathbf{F}}}}_{Cz}{\widehat{\mathit{\rho }}}_C\left(t_{n-\mathrm{1}}\right)+{\widehat{\widehat{\mathbf{A}}}}_C\left(t_{n-\mathrm{1}}\right){\widehat{\widehat{\mathbf{F}}}}_{Cz}{\widehat{\mathit{\rho }}}_C\left(t_{n-\mathrm{1}}\right)\mathrm{d}t\\ & +W_{\mathrm{a}}\underset{\mathrm{0}}{\underbrace{{\widehat{\widehat{\mathbf{F}}}}_{Cz}\mathit{\left\{}{\widehat{\mathit{\rho }}}_S\right(t_{n-\mathrm{1}})\otimes {\widehat{\mathit{\rho }}}_{{\text{pH}}_{\mathrm{2}}}\mathit{\}}}}\mathrm{d}t.\end{array}\end{array}$$

The last term is zero because of the induction hypothesis; see Eqs. (C10) and (B15):

$$\begin{array}{}\text{(C14)}& {\widehat{\widehat{\mathbf{F}}}}_{Cz}\mathit{\left\{}{\widehat{\mathit{\rho }}}_S\right(t_{n-\mathrm{1}})\otimes {\widehat{\mathit{\rho }}}_{{\text{pH}}_{\mathrm{2}}}\mathit{\}}=\mathrm{0}.\end{array}$$

Therefore, using Eq. (C13), we prove the theorem

$$\begin{array}{}\text{(C15)}& \begin{array}{rl}{\widehat{\widehat{\mathbf{F}}}}_{Cz}& {\widehat{\mathit{\rho }}}_C\left(t_n\right)\\ & =\underset{\mathrm{0}}{\underbrace{{\widehat{\widehat{\mathbf{F}}}}_{Cz}{\widehat{\mathit{\rho }}}_C\left(t_{n-\mathrm{1}}\right)}}+{\widehat{\widehat{\mathbf{A}}}}_C\left(t_{n-\mathrm{1}}\right)\underset{\mathrm{0}}{\underbrace{{\widehat{\widehat{\mathbf{F}}}}_{Cz}{\widehat{\mathit{\rho }}}_C\left(t_{n-\mathrm{1}}\right)}}\mathrm{d}t=\mathrm{0},\end{array}\end{array}$$

since ${\widehat{\widehat{\mathbf{F}}}}_{Cz}{\widehat{\mathit{\rho }}}_C\left(t_{n-\mathrm{1}}\right)=\mathrm{0}$ from the induction hypothesis; see Eq. (C10). Therefore, we proved that ${\widehat{\widehat{\mathbf{F}}}}_{Cz}{\widehat{\mathit{\rho }}}_C\left(t\right)=\mathrm{0}$ for any t. Similarly, one can show that ${\widehat{\widehat{\mathbf{F}}}}_{Sz}{\widehat{\mathit{\rho }}}_S\left(t\right)=\mathrm{0}$ for any t.

Appendix D: Asymptotics for the size of the ZQC-reduced SABRE matrix

From Eq. (28), for the dimension of a ZQC subspace we have

$$\begin{array}{}\text{(D1)}& \begin{array}{rl}& \text{dim}\left(\text{ZQC}\right)=\\ & \sum _{F_z=-(N-P)/\mathrm{2}-K}^{(N-P)/\mathrm{2}+K}\mathit{\{}\sum _{m_k=-K}^K\left(\genfrac{}{}{0ex}{}{N-P}{F_z+(N-P)/\mathrm{2}-m_k}\right){\mathit{\}}}^{\mathrm{2}}.\end{array}\end{array}$$

Let us rewrite the inner squared sum:

$$\begin{array}{}\text{(D2)}& \begin{array}{rl}& \text{dim}\left(\text{ZQC}\right)\\ & =\sum _{F_z=-(N-P)/\mathrm{2}-K}^{(N-P)/\mathrm{2}+K}\sum _{m_k=-K}^K\sum _{m_k^{\prime }=-K}^K\left(\genfrac{}{}{0ex}{}{N-P}{F_z+(N-P)/\mathrm{2}-m_k}\right)\\ & \times \left(\genfrac{}{}{0ex}{}{N-P}{F_z+(N-P)/\mathrm{2}-m_k^{\prime }}\right)\\ & =\sum _{m_k=-K}^K\sum _{m_k^{\prime }=-K}^K\mathit{\{}\sum _{F_z=-(N-P)/\mathrm{2}-K}^{(N-P)/\mathrm{2}+K}\left(\genfrac{}{}{0ex}{}{N-P}{F_z+(N-P)/\mathrm{2}-m_k}\right)\\ & \times \left(\genfrac{}{}{0ex}{}{N-P}{F_z+(N-P)/\mathrm{2}-m_k^{\prime }}\right)\mathit{\}},\end{array}\end{array}$$

where at the last stage we changed the order of the sums. Let us consider the inner sum. We introduce new variables $t=F_z+(N-P)/\mathrm{2}-m_k$, $\mathrm{\Delta }=m_k-m_k^{\prime }$:

$$\begin{array}{}\text{(D3)}& \begin{array}{rl}\sum _{F_z=-(N-P)/\mathrm{2}-K}^{(N-P)/\mathrm{2}+K}& \left(\genfrac{}{}{0ex}{}{N-P}{F_z+(N-P)/\mathrm{2}-m_k}\right)\\ & \times \left(\genfrac{}{}{0ex}{}{N-P}{F_z+(N-P)/\mathrm{2}-m_k^{\prime }}\right)\\ & =\sum _{t=-K-m_k}^{N-P+K-m_k}\left(\genfrac{}{}{0ex}{}{N-P}{t}\right)\left(\genfrac{}{}{0ex}{}{N-P}{t+\mathrm{\Delta }}\right).\end{array}\end{array}$$

The sum is defined only if both binomial coefficients are defined. Therefore, we require $\mathrm{0}\le t_{\text{binom}}\le N-P$ and $\mathrm{0}\le t_{\text{binom}}+\mathrm{\Delta }\le N-P$. When combined, we require $t_{\text{binom}}\in [max(\mathrm{0},-\mathrm{\Delta }),min(N-P,N-P-\mathrm{\Delta }\left)\right]$. For example, at Δ>0 both binomial coefficients are defined if $t_{\text{binom}}\in [\mathrm{0},N-P-\mathrm{\Delta }]$. In turn, the sum itself is taken over $t_{\text{sum}}\in [-K-m_k,N-P+K-m_k]$. Fortunately, for all Δ we have t_binom⊂t_sum, and the summation interval always covers t_binom. In this case, Vandermonde's identity is valid:

$$\begin{array}{}\text{(D4)}& \begin{array}{rl}\sum _{t=-K-m_k}^{N-P+K-m_k}& \left(\genfrac{}{}{0ex}{}{N-P}{t}\right)\left(\genfrac{}{}{0ex}{}{N-P}{t+\mathrm{\Delta }}\right)\\ & =\left(\genfrac{}{}{0ex}{}{\mathrm{2}N-\mathrm{2}P}{N-P+\mathrm{\Delta }}\right)\\ & =\left(\genfrac{}{}{0ex}{}{\mathrm{2}N-\mathrm{2}P}{N-P+m_k-m_k^{\prime }}\right).\end{array}\end{array}$$

Thus, we have

$$\begin{array}{}\text{(D5)}& \text{dim}\left(\text{ZQC}\right)=\sum _{m_k=-K}^K\sum _{m_k^{\prime }=-K}^K\left(\genfrac{}{}{0ex}{}{\mathrm{2}N-\mathrm{2}P}{N-P+m_k-m_k^{\prime }}\right).\end{array}$$

If $\sqrt{N-P}\gg K$, the Gauss asymptotics is valid for the binomial coefficient

$$\begin{array}{}\text{(D6)}& \begin{array}{rl}\left(\genfrac{}{}{0ex}{}{\mathrm{2}N-\mathrm{2}P}{N-P+m_k-m_k^{\prime }}\right)& \sim {\displaystyle \frac{{\mathrm{4}}^{N-P}}{\sqrt{\mathit{\pi }(N-P)}}}{e}^{-(m_k-m_k^{\prime }{)}^{\mathrm{2}}/(N-P)}\\ & \approx {\displaystyle \frac{{\mathrm{4}}^{N-P}}{\sqrt{\mathit{\pi }(N-P)}}}.\end{array}\end{array}$$

Therefore, the total sum is

$$\begin{array}{}\text{(D7)}& \begin{array}{rl}\text{dim}\left(\text{ZQC}\right)& \approx {\displaystyle \frac{{\mathrm{4}}^{N-P}}{\sqrt{\mathit{\pi }(N-P)}}}\sum _{m_k=-K}^K\sum _{m_k^{\prime }=-K}^K\mathrm{1}\\ & ={\displaystyle \frac{(\mathrm{2}K+\mathrm{1}{)}^{\mathrm{2}}}{\sqrt{\mathit{\pi }(N-P)}}}{\mathrm{4}}^{N-P}.\end{array}\end{array}$$

The total dimension of the ZQC SABRE matrix is

$$\begin{array}{}\text{(D8)}& \begin{array}{rl}\text{dim}& \left(\text{ZQC SABRE}\right)\\ & =\text{dim}\left({\text{ZQC}}_S\right)+\text{dim}\left({\text{ZQC}}_C\right)\\ & \sim {\displaystyle \frac{(\mathrm{2}K+\mathrm{1}{)}^{\mathrm{2}}}{\sqrt{\mathit{\pi }(N-P)}}}{\mathrm{4}}^{N-P}+{\displaystyle \frac{(\mathrm{2}K+\mathrm{1}{)}^{\mathrm{2}}}{\sqrt{\mathit{\pi }(N-P+\mathrm{2})}}}{\mathrm{4}}^{N-P+\mathrm{2}}\\ & \sim {\displaystyle \frac{(\mathrm{2}K+\mathrm{1}{)}^{\mathrm{2}}}{\sqrt{\mathit{\pi }(N-P)}}}{\mathrm{4}}^{N-P+\mathrm{2}}.\end{array}\end{array}$$

Appendix E: Simulation parameters

The J-coupling constants for the SABRE complex with acetonitrile were taken from Mewis et al. (2015); see Table E1. The J-coupling constants for [¹⁵N,¹³C₄]butyronitrile (Table E2) are based on our experimental measurements. These data will be discussed comprehensively in a subsequent paper (Kiryutin et al., 2026). For simplicity, in all cases we assumed that the J-coupling network between the nuclei within the substrates remains unchanged upon binding to the complex.

Table E1 Simulation NMR parameters used for [¹⁵N,¹³C₂]acetonitrile.

Ha, Hb denote the hydride ligands.
The values of T₁ in brackets correspond to the free substrate.

Table E2 Simulation NMR parameters used for [¹⁵N,¹³C₄]butyronitrile.

Ha, Hb denote the hydride ligands.
The values of T₁ in brackets correspond to the free substrate.
3-CH₂/4-CH₂ are the methylene protons at the 3-¹³C/4-¹³C positions.