Simulation of NMR spectra at zero and ultralow fields from A to Z – a tribute to Prof. Konstantin L'vovich Ivanov

Simulating NMR experiments may appear mysterious and even daunting for those who are new to the field. Yet, broken down into pieces, the process may turn out to be easier than expected. Quite the opposite, it is in fact a powerful and playful means to get insights into the spin dynamics of NMR experiments. In this tutorial paper, we show step by step how some NMR experiments can be simulated, assuming as little prior knowledge from the reader as possible. We focus on the case of NMR at zero and ultralow fields, an emerging modality of NMR in which the spin dynamics are dominated by spin–spin interactions rather than spin–field interactions, as is usually the case with conventional high-field NMR. We first show how to simulate spectra numerically. In a second step, we detail an approach to construct an eigenbasis for systems of spin-$\mathrm{1}/\mathrm{2}$ nuclei at zero field. We then use it to interpret the numerical simulations.

1 Introduction

NMR spectroscopists know well the advantages of performing experiments at the highest possible magnetic field. Increasing magnetic field strength boosts the sensitivity thanks to higher Boltzmann nuclear polarization and higher Larmor frequency (provided the signal linewidth is maintained constant). In addition to this already convincing advantage, higher magnetic fields also imply a larger frequency shift dispersion and therefore easier resolution of individual resonances in crowded spectra. This has motivated the use of ever-increasing magnetic fields (Thayer and Pines, 1987; Schwalbe, 2017; Wikus et al., 2022). The past year has witnessed the implementation of the first spectrometers operating at no less than 28 T, corresponding to a ¹H Larmor frequency of 1.2 GHz. (Wikus et al., 2022) There is no doubt that these new instruments will allow for unprecedented applications.

On the fringe of these great achievements, there is growing interest in the opposite strategy, namely, zero- to ultralow-field (ZULF) NMR, a modality of NMR experiments where the dominant interactions are spin–spin interactions rather than spin–field interactions (Thayer and Pines, 1987; Weitekamp et al., 1983; Blanchard and Budker, 2016; Blanchard et al., 2021; Tayler et al., 2017; Jiang et al., 2021). To realize such conditions, ZULF experiments are not performed in magnets but rather in mumetal magnetic shields that screen magnetic fields originating from the Earth and other surrounding sources, bringing the residual field down to nanotesla (nT) values. In this paper, “zero field” (ZF) designates the regime where heteronuclear spin–spin interactions dominate over spin–field interactions (Zeeman interactions), and the residual spin–field interactions are small enough for the Larmor period to be much longer than the coherence time (Blanchard and Budker, 2016). When this condition is met, decreasing the residual field to even lower values leaves the NMR spectrum unchanged. “Ultralow field” (ULF) designates the regime where the spin–field interactions can be treated as a perturbation to the heteronuclear spin–spin interactions. This typically corresponds to fields on the order of tens to hundreds of nanotesla (Ledbetter et al., 2011). Liquid-state ZULF experiments result in J spectra which do not feature any chemical shift information (Ledbetter et al., 2009). The regime where the intensity of heteronuclear spin–spin interactions is on the order of that of the spin–field interactions occurs typically in the range of microtesla (µT) to tens of microtesla and is referred to as Earth-field NMR (EF-NMR) (Callaghan and Le Gros, 1982; Appelt et al., 2006).

In the simplest form of ZULF experiments, the sample is thermally prepolarized in a permanent magnet (typically 2 T) (Tayler et al., 2017) and subsequently shuttled into the magnetic shields for detection at ZF or ULF. Alternatively, ZULF experiments may be coupled with hyperpolarization techniques (Theis et al., 2012; Butler et al., 2013b; Barskiy et al., 2019; Picazo-Frutos et al., 2023). In particular, parahydrogen-induced polarization (PHIP) has become common as a method for enhancing ZULF signals (Theis et al., 2011, 2012; Butler et al., 2013b). Once the sample is prepolarized (or hyperpolarized), coherences are excited using constant magnetic field pulses rather than radiofrequency (RF) pulses and are usually detected using optically pumped magnetometers (OPMs) rather than inductive coils (Ledbetter et al., 2009). Contrary to high-field instruments, ZULF spectrometers have the advantage of being cheap and relatively easy to assemble (Tayler et al., 2017). They are small enough to sit on a bench and do not require the use of cryogenics (at least if OPMs are used for detection).

Most people who have been introduced to the theory of high-field NMR have first encountered the vector model. The representation of a single-spin system as a vector in 3D space is a powerful tool to build intuition on what happens during an NMR experiment. Then, in a second step, the product operator formalism is necessary to understand the outcome of experiments involving interacting spins. At ZULF, couplings between spins need to be taken into account even to describe the simplest experiment, which consists of detecting the coherence between the singlet S₀ and triplet T₀ states of a pair of J-coupled heteronuclei, e.g., ¹H and ¹³C (Blanchard and Budker, 2016). Polarization oscillates back and forth from one heteronucleus to the other, producing an observable oscillating signal whose frequency is given by the J coupling between the two spins. The outcome of the experiment is simple – a single line at the J-coupling frequency – although it cannot be predicted by the vector model of high-field NMR and Bloch equations. Nonetheless, it is possible to build intuition regarding ZULF experiments in several ways. First, when dealing with two-spin systems, one can define spin operators at ZF in analogy to that at high field so as to translate some of the intuitions from high field to ZULF (Blanchard and Budker, 2016; Butler et al., 2013b). Second, there is a strong analogy between the energy states of electronic spins in atoms and coupled nuclei at ZF (Butler et al., 2013a; Theis et al., 2013). The formalism of addition of angular momenta (widely used in atomic physics and rotational spectroscopy but less frequently in liquid-state NMR) can therefore be used to describe ZULF experiments. Finally, ZULF experiments can be numerically simulated easily, and – as is the case for high-field NMR – simulation provides a playful means to understand NMR experiments (Blanchard et al., 2020; Put et al., 2021). This tutorial paper is focused on the last two approaches.

We present a step-by-step procedure to numerically simulate ZULF spectra in some simple cases. The process is broken down into the following steps:

1. define the experimental sequence,

2. define the spin system,

3. compute the spin Hamiltonian,

4. define the initial state – compute the initial density matrix,

5. propagate the density matrix under the Hamiltonians,

6. extract expectation values from the propagation,

7. Fourier transform the expectation values to obtain a spectrum.

We assume that the reader is familiar with general concepts of NMR and that they are not necessarily used to performing spin dynamics simulations. We take particular care to detail the technical “tricks” which are generally omitted in research papers but are nonetheless essential to performing successful simulations. We present simulated spectra for XA_n spin systems with n between 1 and 5 with several excitation schemes. The spectra are simulated using MATLAB live scripts, which are available in the Supplement (see also Stern and Sheberstov, 2023). The code is abundantly commented and is constructed so as to follow precisely the recipe presented in this paper. Each object and operation presented in this paper can thus be related to lines in the MATLAB code, and vice versa. PDF versions of the live scripts are available. We strongly advise the reader to read the code in parallel to the paper. In a second step, we interpret the simulated results by performing an analytical analysis of the XA_n system using a theoretical framework coming from atomic physics. We show how to construct an eigenbasis and find the selection rules for the allowed transitions. This section is also supported with a code written in Wolfram Mathematica and with a step-by-step link between the text and lines in the code supporting the derived equations (Stern and Sheberstov, 2023).

The reader might wonder whether it makes sense to go through all the details of simulating NMR experiments from scratch while there are powerful simulation packages which are freely available. SpinDynamica (Bengs and Levitt, 2018) and Spinach (Hogben et al., 2011), which run on Mathematica and MATLAB, respectively, are probably the most appropriate tools for simulations at ZULF. The people who have programmed these have already gone through the hurdles of making them efficient and versatile for us, and they even provide code examples for the simulation of NMR spectra at ZULF.¹ However, it is the authors' opinion that performing simple simulations from scratch is the best way to get familiar with the quantum mechanical objects of NMR theory. Once one is confident with these objects and their language, one will make the best use of powerful packages such as SpinDynamica and Spinach. We note that several PhD theses from Pines' group at University of California, Berkeley, present simulations of NMR spectra at ZULF (Theis, 2012; Blanchard, 2014; Sjolander, 2017). These theses contain code examples and are a useful resource for the beginner.

In writing this paper, the authors wish to pay tribute to their regretted lecturer and mentor Prof. Konstantin L'vovich Ivanov, known as Kostya by many, who was taken by COVID-19 on 5 March 2021 (Yurkovskaya and Bodenhausen, 2021). Kirill Sheberstov had Konstantin Ivanov as a PhD co-supervisor, performing research on long-lived states, parahydrogen-induced polarization, and chemically induced dynamic nuclear polarization (CIDNP). Konstantin Ivanov's deep understanding of the underlying physics allowed his group to work in very different directions, e.g., to combine CIDNP and ZULF NMR. During his PhD in Sami Jannin's team in Lyon, France, Quentin Stern collaborated with Konstantin Ivanov on a research project. In the course of their collaboration, Konstantin Ivanov gave Quentin Stern guidance on how to simulate experiments at ZF. A few pieces of advice turned into precious teachings for Quentin Stern. Sadly, these teachings were brutally interrupted by Konstantin Ivanov's death. Konstantin Ivanov's kindness and availability to give help and advice will forever remain an example for Quentin Stern and Kirill Sheberstov.

2 Theory – numerical simulation of spin dynamics

2.1 Define the experimental sequence

Most 1D NMR experiments can be broken down into three steps:

$$\mathrm{preparation}-\mathrm{mixing}-\mathrm{detection}.$$

During the preparation, some nuclear polarization is acquired by letting the sample rest in a strong magnetic field (in most conventional experiments). Mixing consists of bringing the system to a non-stationary state whose oscillations are recorded during detection. In common high-field NMR experiments, all the steps are performed in a strong magnet with a nearly constant magnetic field. Nuclear polarization is spontaneously acquired due to the high magnetic field, and both the mixing and detection are performed through the same RF coil using Faraday induction. At ZULF, there is no nuclear polarization, so the preparation has to be performed in different conditions. A common method is to shuttle the sample between a region of high field and a region of ZULF.

Figure 1 shows a typical experimental setup. A permanent magnet is used to prepolarize the sample, which is connected with the magnetic shields by a guiding solenoid coil. This coil ensures that the sample experiences a magnetic field with constant direction and sufficient strength during the transfer from the region of high field to inside the magnetic shields (i.e., the coil ensures an adiabatic transfer). Once the sample arrives in the magnetic shields at the location of detection, the Helmholtz coil continues to produce a magnetic field in the same direction as the solenoid, and the spin system is still distributed into Zeeman populations (Blanchard and Budker, 2016; Tayler et al., 2017). All the steps detailed until here are part of the preparation. In practice, the guiding solenoid and the Helmholtz coil produce a magnetic field which is much weaker than the prepolarizing magnet. However, this will not be taken into account in the simulation: we consider that the sample spends enough time in the prepolarizing magnet to reach Boltzmann equilibrium and that the transfer is sufficiently fast for us to neglect the change in polarization during the transfer.

A further step can optionally be added to the preparation, which consists of ramping down the magnetic field produced by the Helmholtz coil to bring the spins adiabatically to ZF. We will refer to experiments which include or do not include this step as the adiabatic field drop and sudden field drop experiments, respectively. In the case of sudden field drop experiments, the mixing step simply consists of switching off the magnetic field non-adiabatically (that is, fast enough to be considered instantaneous with respect to the evolution of the spin system). In the case of adiabatic field drop experiments, the sample is already at ZF at the end of preparation, so populations have to be mixed by applying a magnetic field pulse before any signal can be detected. This is analogous to high-field pulses except that it uses constant magnetic field rather than RF pulses. After the mixing, the oscillating magnetic field generated by the sample is detected by an optical magnetometer. In Fig. 1, the magnetometer is represented below the sample, that is, aligned along the z axis with respect to the sample. We assume that the OPM is configured so as to be sensitive to magnetic field along the z axis and that the spins are initially prepolarized along the same axis. Defining this axis as z is a natural choice for high-field NMR spectroscopists, but note that other conventions exist (see, for example, Ledbetter et al., 2011). During detection, a weak magnetic field may be applied, either along the z axis or along an orthogonal axis. In the latter case, the experiment is said to be performed under the ULF regime. In the absence of applied magnetic field (and provided the residual magnetic field is properly zeroed at the location of the sample), the experiment is said to be performed under the ZF regime.

Figure 1 (a) Typical experimental setup for ZULF experiments. Note that the sample is represented in two places in the same drawing even if there is a single sample. (b) Schemes of the experimental sequences for measurements at ZF using a sudden field drop or an adiabatic field drop followed by a pulse of static magnetic field.

In summary, there are several possible combinations of experimental schemes. All of them start with prepolarizing the spins at high magnetic field. After the sample is transported into the magnetic shields, the field is dropped either suddenly or adiabatically, in which case a magnetic field pulse is applied. Finally, the oscillating magnetic field produced by the sample is detected along the z axis, with or without a weak magnetic field applied along the x axis. In the remaining of the paper, these sequences presented in Fig. 3b will be broken down into the following steps:

1. prepolarization,

2. transfer and coherence excitation, and

3. detection.

2.2 Define the spin system

This step consists of listing the different magnetic sites of the molecule whose ZULF spectrum is to be simulated and the interactions which the spins are subject to. This paper is concerned with small molecules in the liquid state. As is the case for high-field NMR, dipolar interactions are averaged out by rapid molecular tumbling and need not be taken into account (except as stochastic perturbation if one intends to include relaxation effects). Therefore, only the J coupling and the Zeeman interactions are considered here.

In this paper, we consider spin systems of the form XA_n, where X is a ¹³C spin coupled to n equivalent ¹H spins A through a coupling J_AX. The A spins are coupled together through J_AA.

2.3 Compute the spin Hamiltonian

The Hamiltonian is the operator which represents the total energy of the system. Information about the spin system is mathematically encoded in the spin Hamiltonian. We will first present how the Hamiltonian for the Zeeman interaction of a single spin is computed based on Pauli matrices. Then, we will present the construction of a two-spin system using the Kronecker product of individual spin spaces to compute the Zeeman and the J-coupling Hamiltonians. Finally, we will show how the procedure is extended to an arbitrary number of spins.

Let us first assume that the system contains a single spin $\mathrm{1}/\mathrm{2}$ interacting with the magnetic field B along the z axis. The state of any spin system can be represented as a linear combination of basis vectors, which are called “kets” in Dirac's notation and are represented by the symbols |〉. For a single spin $\mathrm{1}/\mathrm{2}$, two basis kets are necessary to represent the state of the system. We chose to represent the spin system in the Zeeman basis:

$$\begin{array}{}\text{(1)}& B_{\mathrm{Z}}^{\mathrm{2}}=\left\{\left|\mathit{\alpha }〉\right.,\left|\mathit{\beta }〉\right.\right\}=\left\{\left(\begin{array}{c}\mathrm{1}\\ \mathrm{0}\end{array}\right),\left(\begin{array}{c}\mathrm{0}\\ \mathrm{1}\end{array}\right)\right\}.\end{array}$$

The |α and |β states correspond to the spin being parallel and antiparallel with the magnetic field, respectively. The general state in which the spin may be found is a linear combination of these two basis states. Because these states and their associated kets form an orthonormal basis, their vector representation have the canonical form with only 0 and 1 coefficients. The space spanned by these two vectors is called a “Hilbert space” and has dimension 2, as indicated by the superscript in $B_{\mathrm{Z}}^{\mathrm{2}}$. Note that the choice of the Zeeman basis is convenient for numerical simulation but is not necessary. For example, one may use the coupled basis, which will be presented and used in Sect. 4. The same basis may be used to perform simulations at high field or ZULF, although one particular basis might turn out to be more convenient.

The angular momentum of a single spin is associated with the spin angular momentum operators, which can be represented as a vector with three Cartesian components:

$$\begin{array}{}\text{(2)}& \widehat{\mathit{I}}=\left({\widehat{I}}_x,{\widehat{I}}_y,{\widehat{I}}_z\right).\end{array}$$

These operators act on the Zeeman states in certain ways, e.g., ${\widehat{I}}_x\left|\mathit{\alpha }〉\right.=\frac{\mathrm{1}}{\mathrm{2}}\left|\mathit{\beta }〉\right.$. To summarize the set of rules, it is convenient to use the matrix representation of the operators, with the matrix elements determined by the action of the operator on the |α and |β states: $I_{\mathit{\mu }}^{rs}=〈r\left|{\widehat{I}}_k\right|s〉$, where $r,s\in \left\{\mathit{\alpha },\mathit{\beta }\right\};\mathit{\mu }\in \mathit{\{}x,y,z\mathit{\}}$. This definition makes use of 〈|, i.e., the “bra”, an object which is complementary to the ket and corresponds to the “Hermitean conjugate” of the ket. In the matrix representation of quantum mechanics, the Hermitean conjugate of a ket corresponds to the complex transpose of the vector representing the ket. The matrix representations of operators in quantum mechanics is very important for performing simulations, as they are constructed in such a way that any state of operation on the quantum system can be represented using linear algebra. The matrix representations of the three Cartesian components of the spin angular momentum operators are proportional to Pauli matrices ${\widehat{\mathit{\sigma }}}_x$, ${\widehat{\mathit{\sigma }}}_y$, and ${\widehat{\mathit{\sigma }}}_z$:

$$\begin{array}{}\text{(3)}& \begin{array}{rl}& {\widehat{I}}_x={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\widehat{\mathit{\sigma }}}_x={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left(\begin{array}{cc}\mathrm{0}& \mathrm{1}\\ \mathrm{1}& \mathrm{0}\end{array}\right),\\ & {\widehat{I}}_y={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\widehat{\mathit{\sigma }}}_y={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left(\begin{array}{cc}\mathrm{0}& -i\\ i& \mathrm{0}\end{array}\right),\\ & {\widehat{I}}_z={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\widehat{\mathit{\sigma }}}_z={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left(\begin{array}{cc}\mathrm{1}& \mathrm{0}\\ \mathrm{0}& -\mathrm{1}\end{array}\right).\end{array}\end{array}$$

The interaction of a single spin with a magnetic field B is given by the Zeeman Hamiltonian:

$$\begin{array}{}\text{(4)}& {\widehat{H}}_{\mathrm{Z}}=-\mathit{\gamma }\mathit{B}\cdot \widehat{\mathit{I}}=-\mathit{\gamma }\left(B_x{\widehat{I}}_x+B_y{\widehat{I}}_y+B_z{\widehat{I}}_z\right),\end{array}$$

where γ is the gyromagnetic ratio of the spin (in rad s⁻¹ T⁻¹). The dot product of the vectors of the magnetic field and of the spin angular momentum (vectors and vector operators are denoted in bold throughout the text) is expanded on the right member of Eq. (4). Note that we have omitted the reduced Planck constant ℏ in Eq. (4), which implies that the energy is expressed in radians per second (rad s⁻¹) rather than in joules. This is the case throughout this paper. In many cases, the magnetic field is aligned with one of the axes. If it points along the z axis, i.e., $\mathit{B}=\left(\begin{array}{ccc}\mathrm{0}& \mathrm{0}& B_{\mathrm{0}}\end{array}\right)$, Eq. (4) simplifies to

$$\begin{array}{}\text{(5)}& {\widehat{H}}_{\mathrm{Z}}=-\mathit{\gamma }{B}_{\mathrm{0}}{\widehat{I}}_z={\mathit{\omega }}^{\mathrm{0}}{\widehat{I}}_z={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left(\begin{array}{cc}+{\mathit{\omega }}^{\mathrm{0}}& \mathrm{0}\\ \mathrm{0}& -{\mathit{\omega }}^{\mathrm{0}}\end{array}\right),\end{array}$$

where ${\mathit{\omega }}^{\mathrm{0}}=-\mathit{\gamma }{B}_{\mathrm{0}}$ is the Larmor frequency of the spin. This expression is valid regardless of the intensity of the magnetic field, i.e., at high field as well as at ZULF. The Zeeman states, |α and |β, which correspond to the spin being parallel and antiparallel with the magnetic field, respectively, are eigenstates of the Zeeman Hamiltonian; that is, they satisfy the eigenequations ${\widehat{H}}_{\mathrm{Z}}\left|\mathit{\alpha }〉\right.=+\mathrm{1}/\mathrm{2}\left|\mathit{\alpha }〉\right.$ and ${\widehat{H}}_{\mathrm{Z}}\left|\mathit{\beta }〉\right.=-\mathrm{1}/\mathrm{2}\left|\mathit{\beta }〉\right.$. The eigenstates of a Hamiltonian are of particular importance; they are states which do not evolve under the effect of that Hamiltonian (ignoring the accumulation of the phase factor, which turns out to be irrelevant in most experiments), i.e., stationary states.

The single spin whose Hamiltonian is given by Eq. (5) lives in a Hilbert space of dimension 2. To represent a pair of spins $\mathrm{1}/\mathrm{2}$, we need to use a Hilbert space with a dimension of 4. To do so, we redefine the angular momentum operators in this higher-dimension space. The matrix representations of the angular momentum operators ${\widehat{I}}_{\mathrm{1}x}$, ${\widehat{I}}_{\mathrm{1}y}$, and ${\widehat{I}}_{\mathrm{1}z}$ and ${\widehat{I}}_{\mathrm{2}x}$, ${\widehat{I}}_{\mathrm{2}y}$, and ${\widehat{I}}_{\mathrm{2}z}$ of spin 1 and spin 2, respectively, are given by the Kronecker product of matrices of single-spin angular momentum operator and the identity operator, in the appropriate order. For the z-axis angular momentum operators, we have

$$\begin{array}{}\text{(6)}& \begin{array}{rl}{\widehat{I}}_{\mathrm{1}z}& ={\widehat{I}}_z\otimes \widehat{\mathrm{1}}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left(\begin{array}{cc}\mathrm{1}& \mathrm{0}\\ \mathrm{0}& -\mathrm{1}\end{array}\right)\otimes \left(\begin{array}{cc}\mathrm{1}& \mathrm{0}\\ \mathrm{0}& \mathrm{1}\end{array}\right)\\ & ={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left(\begin{array}{cccc}+\mathrm{1}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& +\mathrm{1}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& -\mathrm{1}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& -\mathrm{1}\end{array}\right)\end{array}\end{array}$$

and

$$\begin{array}{}\text{(7)}& \begin{array}{rl}{\widehat{I}}_{\mathrm{2}z}& =\widehat{\mathrm{1}}\otimes {\widehat{I}}_z=\left(\begin{array}{cc}\mathrm{1}& \mathrm{0}\\ \mathrm{0}& \mathrm{1}\end{array}\right)\otimes {\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left(\begin{array}{cc}\mathrm{1}& \mathrm{0}\\ \mathrm{0}& -\mathrm{1}\end{array}\right)\\ & ={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left(\begin{array}{cccc}+\mathrm{1}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& -\mathrm{1}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& +\mathrm{1}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& -\mathrm{1}\end{array}\right).\end{array}\end{array}$$

Similar expressions are obtained for the matrices of x and y operators. They are not shown here but are available in many textbooks (Hore et al., 2015; Levitt, 2013). Here, we have used the following convention for the Kronecker product

$$\begin{array}{}\text{(8)}& \begin{array}{rl}\left(\begin{array}{cc}a& b\\ c& d\end{array}\right)& \otimes \left(\begin{array}{cc}\mathit{\alpha }& \mathit{\beta }\\ \mathit{\gamma }& \mathit{\delta }\end{array}\right)\\ & =\left(\begin{array}{cc}a\left(\begin{array}{cc}\mathit{\alpha }& \mathit{\beta }\\ \mathit{\gamma }& \mathit{\delta }\end{array}\right)& b\left(\begin{array}{cc}\mathit{\alpha }& \mathit{\beta }\\ \mathit{\gamma }& \mathit{\delta }\end{array}\right)\\ c\left(\begin{array}{cc}\mathit{\alpha }& \mathit{\beta }\\ \mathit{\gamma }& \mathit{\delta }\end{array}\right)& d\left(\begin{array}{cc}\mathit{\alpha }& \mathit{\beta }\\ \mathit{\gamma }& \mathit{\delta }\end{array}\right)\end{array}\right)\\ & =\left(\begin{array}{cccc}a\mathit{\alpha }& a\mathit{\beta }& b\mathit{\alpha }& b\mathit{\beta }\\ a\mathit{\gamma }& a\mathit{\delta }& b\mathit{\gamma }& b\mathit{\delta }\\ c\mathit{\alpha }& c\mathit{\beta }& d\mathit{\alpha }& d\mathit{\beta }\\ c\mathit{\gamma }& c\mathit{\delta }& d\mathit{\gamma }& d\mathit{\delta }\end{array}\right).\end{array}\end{array}$$

The two operators defined by Eqs. (6) and (7) are the same as the one given by Eq. (5), except that the world of spin 1 now contains spin 2, and vice versa. This representation corresponds to a basis that is the Kronecker product of the basis of the individual spins.

$$\begin{array}{}\text{(9)}& B_{\mathrm{Z}}^{\mathrm{4}}=B_{\mathrm{Z}}^{\mathrm{2}}\otimes B_{\mathrm{Z}}^{\mathrm{2}}=\left\{\left|\mathit{\alpha }\mathit{\alpha }〉\right.,\left|\mathit{\alpha }\mathit{\beta }〉\right.,\left|\mathit{\beta }\mathit{\alpha }〉\right.,\left|\mathit{\beta }\mathit{\beta }〉\right.\right\}\end{array}$$

For the case where the magnetic field points along the z axis, the total Zeeman Hamiltonian for the two spins can now be computed using Eq. (5) in the basis of Eq. (9) as the sum of the two Zeeman Hamiltonians:

$$\begin{array}{}\text{(10)}& \begin{array}{rl}{\widehat{H}}_{\mathrm{Z}}& ={\widehat{H}}_{\mathrm{Z},\mathrm{1}}+{\widehat{H}}_{\mathrm{Z},\mathrm{2}}={\mathit{\omega }}_{\mathrm{1}}^{\mathrm{0}}{\widehat{I}}_{\mathrm{1}z}+{\mathit{\omega }}_{\mathrm{2}}^{\mathrm{0}}{\widehat{I}}_{\mathrm{2}z}\\ & ={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left(\begin{array}{cccc}+{\mathit{\omega }}_{\mathrm{1}}^{\mathrm{0}}+{\mathit{\omega }}_{\mathrm{2}}^{\mathrm{0}}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& +{\mathit{\omega }}_{\mathrm{1}}^{\mathrm{0}}-{\mathit{\omega }}_{\mathrm{2}}^{\mathrm{0}}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& -{\mathit{\omega }}_{\mathrm{1}}^{\mathrm{0}}+{\mathit{\omega }}_{\mathrm{2}}^{\mathrm{0}}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& -{\mathit{\omega }}_{\mathrm{1}}^{\mathrm{0}}-{\mathit{\omega }}_{\mathrm{2}}^{\mathrm{0}}\end{array}\right)\end{array},\end{array}$$

where ${\mathit{\omega }}_{\mathrm{1}}^{\mathrm{0}}$ and ${\mathit{\omega }}_{\mathrm{2}}^{\mathrm{0}}$ are the Larmor frequencies of spin 1 and 2, respectively. Note that in a Hilbert space of several spins, it is useful to define projections of total angular momentum operators:

$$\begin{array}{}\text{(11)}& \begin{array}{rl}& {\widehat{I}}_x={\widehat{I}}_{\mathrm{1}x}+{\widehat{I}}_{\mathrm{2}x},\\ & {\widehat{I}}_y={\widehat{I}}_{\mathrm{1}y}+{\widehat{I}}_{\mathrm{2}y},\\ & {\widehat{I}}_z={\widehat{I}}_{\mathrm{1}z}+{\widehat{I}}_{\mathrm{2}z}.\end{array}\end{array}$$

Note that these operators are represented by the same symbol as their equivalent in the single-spin Hilbert space (see Eq. 3). It should be clear from the context whether the operator corresponds to a single-spin or multiple-spin Hilbert space. Where confusion may remain, we will indicate the dimension of the space on which the operator acts.

At this point, the two spins are represented in a common space, but they do not interact. The J-coupling Hamiltonian for the pair of spins is given by

$$\begin{array}{}\text{(12)}& \begin{array}{rl}{\widehat{H}}_J& =\mathrm{2}\mathit{\pi }J{\widehat{\mathit{I}}}_{\mathrm{1}}\cdot {\widehat{\mathit{I}}}_{\mathrm{2}}=\mathrm{2}\mathit{\pi }J\left({\widehat{I}}_{\mathrm{1}x}{\widehat{I}}_{\mathrm{2}x}+{\widehat{I}}_{\mathrm{1}y}{\widehat{I}}_{\mathrm{2}y}+{\widehat{I}}_{\mathrm{1}z}{\widehat{I}}_{\mathrm{2}z}\right)\\ & =\mathit{\pi }J\left(\begin{array}{cccc}\mathrm{1}/\mathrm{2}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& -\mathrm{1}/\mathrm{2}& \mathrm{1}& \mathrm{0}\\ \mathrm{0}& \mathrm{1}& -\mathrm{1}/\mathrm{2}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{1}/\mathrm{2}\end{array}\right),\end{array}\end{array}$$

where J is the J coupling between the two spins (in Hz). Compared with the Zeeman Hamiltonian (see Eq. 10), the J-coupling Hamiltonian has the particularity to have off-diagonal elements in the {|αβ,|βα} subspace, which implies that the J interaction mixes the |αβ and |βα states. In other words, due to the J interaction, these two states are no longer eigenstates of the spin system.

In the case of a system of n spins $\mathrm{1}/\mathrm{2}$, the same procedure can be applied to define the angular momentum operators and the Hamiltonians. These operators can be represented as 2ⁿ×2ⁿ matrices. Their Zeeman basis can be constructed as in Eq. (9), taking all possible combinations of |α and |β states of the individuals spins. Equations (6) and (7) generalize to

$$\begin{array}{}\text{(13)}& {\widehat{I}}_{kz}={\otimes }_{l=\mathrm{1}}^n{\widehat{u}}_{lz}^{\mathrm{2}\times \mathrm{2}}\text{where}{\widehat{u}}_{lz}^{\mathrm{2}\times \mathrm{2}}=\left\{\begin{array}{l}{\widehat{\mathrm{1}}}^{\mathrm{2}\times \mathrm{2}}\mathrm{if}l\ne k,\\ {\widehat{I}}_z^{\mathrm{2}\times \mathrm{2}}\mathrm{if}l=k,\end{array}\right.\end{array}$$

where and ${\widehat{I}}_{kz}$ is the z angular momentum operator of spin k in an n-spin Hilbert space, and ${\widehat{\mathrm{1}}}^{\mathrm{2}\times \mathrm{2}}$ and ${\widehat{I}}_z^{\mathrm{2}\times \mathrm{2}}$ are the identity operator and the z angular momentum operator in a single-spin Hilbert space, respectively. The z projection of total angular momentum operators is given by

$$\begin{array}{}\text{(14)}& {\widehat{I}}_z=\sum _{l=\mathrm{1}}^n{\widehat{I}}_{lz}.\end{array}$$

Equations (13) and (14) are shown for z operators but apply similarly for x and y operators. The reader is encouraged to compute the matrix representations of these operators using the MATLAB codes provided in the Supplement. The Zeeman Hamiltonian for a system of n spins is given by

$$\begin{array}{}\text{(15)}& \begin{array}{rl}{\widehat{H}}_{\mathrm{Z}}& =-\sum _{l=\mathrm{1}}^n{\mathit{\gamma }}_l\mathit{B}\cdot {\widehat{\mathit{I}}}_l\\ & =-\sum _{l=\mathrm{1}}^n{\mathit{\gamma }}_l\left(B_x{\widehat{I}}_{lx}+B_y{\widehat{I}}_{ly}+B_z{\widehat{I}}_{lz}\right),\end{array}\end{array}$$

where γ_l is the gyromagnetic ratio of spin l. The J Hamiltonian in the same space is given by

$$\begin{array}{}\text{(16)}& \begin{array}{rl}{\widehat{H}}_J& =\mathrm{2}\mathit{\pi }\sum _{l>k}^n{J}_{lk}{\widehat{\mathit{I}}}_l\cdot {\widehat{\mathit{I}}}_k\\ & =\mathrm{2}\mathit{\pi }\sum _{l>k}^n{J}_{lk}\left({\widehat{I}}_{lx}{\widehat{I}}_{kx}+{\widehat{I}}_{ly}{\widehat{I}}_{ky}+{\widehat{I}}_{lz}{\widehat{I}}_{kz}\right),\end{array}\end{array}$$

where J_lk is the J coupling between spins l and k (in Hz). Because a spin is not J coupled to itself, the sum in Eq. (16) does not include terms with l=k. Furthermore, to avoid counting terms twice, terms with l<k are not included either, leaving only l>k terms. The expression of the Zeeman Hamiltonian and J Hamiltonian in Eqs. (15) and (16), respectively, are valid both at high field and at ZULF. What makes the difference between the two regimes is the relative intensity of the two contributions.

2.4 Define the initial state – compute the initial density matrix

The state of a spin system during an NMR experiment is described by a density operator. If |ψ is a ket representing the state of the system as a linear combination of basis states (like those defined in Eqs. 1 and 9), the density operator is given by

$$\begin{array}{}\text{(17)}& \widehat{\mathit{\rho }}=\stackrel{\mathrm{‾}}{\left|\mathit{\psi }〉\right.〈\left.\mathit{\psi }\right|},\end{array}$$

where the upper bar represents the ensemble average over all identical spin systems in the sample – the operation performed by the density operator. This averaging makes the density operator formalism well-suited for NMR, where the experiment consists of observing a large number of identical spin systems at the same time rather than a single spin system. The matrix representation of the density operator (and of any other spin operator) is achieved by calculating all the matrix elements ${\mathit{\rho }}_{rs}=〈\left.r\right|\widehat{\mathit{\rho }}\left|s〉\right.$, where |r and |s are basis states. For example, the matrix representation of the density operator for the |α and |β states of a single spin yields

$$\begin{array}{}\text{(18)}& \begin{array}{rl}& {\widehat{\mathit{\rho }}}_{\mathit{\alpha }}=\stackrel{\mathrm{‾}}{\left|\mathit{\alpha }〉\right.〈\left.\mathit{\alpha }\right|}=\left(\begin{array}{c}\mathrm{1}\\ \mathrm{0}\end{array}\right)\left(\begin{array}{cc}\mathrm{1}& \mathrm{0}\end{array}\right)=\left(\begin{array}{cc}\mathrm{1}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}\end{array}\right),\\ & {\widehat{\mathit{\rho }}}_{\mathit{\beta }}=\stackrel{\mathrm{‾}}{\left|\mathit{\beta }〉\right.〈\left.\mathit{\beta }\right|}=\left(\begin{array}{c}\mathrm{0}\\ \mathrm{1}\end{array}\right)\left(\begin{array}{cc}\mathrm{0}& \mathrm{1}\end{array}\right)=\left(\begin{array}{cc}\mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{1}\end{array}\right).\end{array}\end{array}$$

To start a simulation, we need to determine the density matrix of the system at the initial point of the experiment. We assume that the sample has spent enough time in the prepolarizing magnet to reach thermal equilibrium; that is, the spin system follows Boltzmann's distribution of states. In this case, the density matrix is given by

$$\begin{array}{}\text{(19)}& {\widehat{\mathit{\rho }}}_{\mathrm{eq}}={\displaystyle \frac{\exp \left(-\frac{\widehat{H}}{k_{\mathrm{B}}T}\right)}{Z}},\end{array}$$

where $\widehat{H}$, k_B, and T are the Hamiltonian operator of the spin system, Boltzmann's constant, and the temperature, respectively. Operation exp ( ) denotes the matrix exponentiation. Note that this operation does not consist of applying f(x)=exp (x) to each element of the matrix. It is a more complex operation, which is realized in MATLAB by the built-in function expm (rather than exp). Z is a normalization constant, which ensures that the density matrix has unit trace. It is given by

$$\begin{array}{}\text{(20)}& Z=\mathrm{Tr}\left\{\exp \left(-{\displaystyle \frac{\widehat{H}}{k_{\mathrm{B}}T}}\right)\right\}.\end{array}$$

The prepolarizing step of the experiments that we intend to simulate occurs in a strong magnetic field (in the sense that the Zeeman interaction is largely dominating all other interactions), as in a standard high-field experiment. In this case, we can compute the thermal equilibrium by taking only the Zeeman terms into account. For a single spin with Larmor frequency ω⁰ and gyromagnetic ratio γ in prepolarizing field B_p, the thermal equilibrium density matrix yields

$$\begin{array}{}\text{(21)}& \begin{array}{rl}{\widehat{\mathit{\rho }}}_{\mathrm{eq}}& ={\displaystyle \frac{\exp \left(-\frac{\mathrm{\hslash }{\mathit{\omega }}^{\mathrm{0}}{\widehat{I}}_z}{k_{\mathrm{B}}T}\right)}{Z}}={\displaystyle \frac{\exp \left(+\frac{\mathrm{\hslash }\mathit{\gamma }{B}_{\mathrm{p}}{\widehat{I}}_z}{k_{\mathrm{B}}T}\right)}{Z}}\\ & =\left(\begin{array}{cc}\frac{\mathrm{1}+P}{\mathrm{2}}& \mathrm{0}\\ \mathrm{0}& \frac{\mathrm{1}-P}{\mathrm{2}}\end{array}\right)={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\widehat{\mathrm{1}}+P{\widehat{I}}_z,\end{array}\end{array}$$

where P is the polarization of the nucleus along the z axis (for positive γ, it corresponds to the population excess of the |α state with respect to the |β state), defined by

$$\begin{array}{}\text{(22)}& P=\tanh \left({\displaystyle \frac{\mathrm{\hslash }\mathit{\gamma }{B}_{\mathrm{p}}}{\mathrm{2}{k}_{\mathrm{B}}T}}\right).\end{array}$$

Note that the use of ℏ in the expression of the Hamiltonian (i.e., expressing the energy in joules) cannot be avoided here, to ensure consistency of units. To obtain the expression on the right-hand side of Eq. (21), we have jumped several steps of calculation which are all based on the definition of polarization. This expression for the density matrix is exact for a spin whose only interaction is the Zeeman interaction, which we have assumed here.

For an n-spin system, we take the Kronecker product of density matrices of individual spins ${\widehat{\mathit{\rho }}}_{\mathrm{eq},l}^{\mathrm{2}\times \mathrm{2}}$.

$$\begin{array}{}\text{(23)}& \begin{array}{rl}{\widehat{\mathit{\rho }}}_{\mathrm{eq}}& \approx {\otimes }_{l=\mathrm{1}}^n{\widehat{\mathit{\rho }}}_{\mathrm{eq},l}^{\mathrm{2}\times \mathrm{2}}={\otimes }_{l=\mathrm{1}}^n\left({\displaystyle \frac{{\widehat{\mathrm{1}}}^{\mathrm{2}\times \mathrm{2}}}{\mathrm{2}}}+P_l{\widehat{I}}_z^{\mathrm{2}\times \mathrm{2}}\right)\\ & \approx {\displaystyle \frac{\widehat{\mathrm{1}}}{{\mathrm{2}}^n}}+{\displaystyle \frac{\mathrm{1}}{{\mathrm{2}}^{n-\mathrm{1}}}}\sum _{l=\mathrm{1}}^n{P}_l{\widehat{I}}_{lz}\end{array}\end{array}$$

The first approximation in this expression comes from the fact that it neglects all spin–spin interactions. The second approximation comes from the fact that we neglect multiple-spin order terms, which are proportional to the products P_lP_k, where P_l and P_k are the polarizations of spin l and k, respectively. This approximation is valid unless the system is highly polarized, which is the case even at very high field (without hyperpolarization). To avoid confusion, we specified that the operators ${\widehat{\mathit{\rho }}}_{\mathrm{eq},l}^{\mathrm{2}\times \mathrm{2}}$, ${\widehat{\mathrm{1}}}^{\mathrm{2}\times \mathrm{2}}$, and ${\widehat{I}}_z^{\mathrm{2}\times \mathrm{2}}$ act on a single-spin Hilbert space (2×2 matrix). On the contrary, the operators $\widehat{\mathrm{1}}$ and ${\widehat{I}}_{lz}$ act on spin states of n spins, and accordingly their matrix representations have dimension of 2ⁿ×2ⁿ (for spins $\mathrm{1}/\mathrm{2}$). As shown by Eq. (23), one may compute the density matrix either using the Kronecker product of operators in a single-spin Hilbert space or by summing the operators in a Hilbert space of n spins.

In many textbooks (Hore et al., 2015; Levitt, 2013), one encounters simplified expressions of the density operator. First, it is common to remove the identity component:

$$\begin{array}{}\text{(24)}& {\widehat{\mathit{\rho }}}_{\mathrm{eq}}\to {\widehat{\mathit{\rho }}}_{\mathrm{eq}}-{\displaystyle \frac{\widehat{\mathrm{1}}}{{\mathrm{2}}^n}},\end{array}$$

where n is the number of spins in the system. Because all operators commute with the identity, this does not affect the result of propagation. The resulting expression is simpler (${\widehat{\mathit{\rho }}}_{\mathrm{eq}}=P{\widehat{I}}_z^{\mathrm{2}\times \mathrm{2}}$ for a single spin), which is convenient for calculations by hand. It may also make the numerical propagation faster and more precise. Another common simplification is to drop the polarization factor. For a single spin, the two combined simplifications yield

$$\begin{array}{}\text{(25)}& {\widehat{\mathit{\rho }}}_{\mathrm{eq}}\to {\widehat{I}}_z.\end{array}$$

Simplifications are useful, but they should be handled with care. The polarization factor P is different for spins with different gyromagnetic ratio. If it is dropped without introducing further corrections, the relative sizes of the population of spins with different gyromagnetic ratios will not be respected. In the simulations presented here, we will compute the initial density matrix using the transformation of Eq. (24) but not that of Eq. (25).

2.5 Propagate the density matrix under the Hamiltonians

We have seen how to compute the initial density matrix and the matrix representation of the Hamiltonian. We now describe how the evolution of the system (represented by the density matrix) evolves with time under a given Hamiltonian. This will be used at several steps of the simulation: when the sample is brought adiabatically to ZF, during the pulse, and during the signal measurement.

The evolution of a quantum system with time is given by the time-dependent Schrödinger equation. Its equivalent for the evolution of density matrix is the Liouville–von Neumann equation $\frac{\mathrm{d}}{\mathrm{d}t}\widehat{\mathit{\rho }}\left(t\right)=-i\left[\widehat{H}\left(t\right),\widehat{\mathit{\rho }}\left(t\right)\right]$, which has the solution

$$\begin{array}{}\text{(26)}& \widehat{\mathit{\rho }}\left(t\right)=\widehat{U}\left(t\right){\widehat{\mathit{\rho }}}_{\mathrm{0}}{\widehat{U}}^{-\mathrm{1}}\left(t\right),\end{array}$$

where ${\widehat{\mathit{\rho }}}_{\mathrm{0}}$ is the density matrix at t=0, and $\widehat{U}$ is the propagator during time t, which is defined as

$$\begin{array}{}\text{(27)}& \widehat{U}\left(t\right)=\exp \left(-i\widehat{H}t\right),\end{array}$$

where $\widehat{H}$ is the total Hamiltonian. The operation of Eq. (26) “takes” the spin system from ${\widehat{\mathit{\rho }}}_{\mathrm{0}}$ to $\widehat{\mathit{\rho }}\left(t\right)$. Again, note that exp ( ) denotes the matrix exponentiation and not the element-by-element exponentiation. An important case of propagator is the rotation operator. For an angular momentum operator ${\widehat{I}}_{\mathit{\mu }}$, with $\mathit{\mu }\in \mathit{\{}x,y,z\mathit{\}}$, the propagator $\exp \left(-i{\widehat{I}}_{\mathit{\mu }}\mathit{\theta }\right)$ is called a rotation operator; it represents a rotation of the spins of angle θ around axis μ, when applied to the density matrix using Eq. (26). For a single spin subject to a static magnetic field along the z axis, the total Hamiltonian is the Zeeman Hamiltonian (see Eq. 5) which causes the spin to rotate around the z axis; this rotation can be expressed using the rotation operator $\exp \left(-i\widehat{H}t\right)=\exp \left(-i{\mathit{\omega }}_{\mathrm{0}}{\widehat{I}}_{z}t\right)$ with angle ω₀t.

The brute force calculation of the exponential operator in an arbitrary basis is computationally challenging as it requires calculating the Taylor expansion of the $\widehat{H}$ operator, which can be calculated analytically only in a few cases, like the propagators of the I_z, I_x, and I_y operators. To avoid this, the calculation of the propagator (Eq. 27) is usually performed by diagonalizing the Hamiltonian and then taking the complex exponent for each of its eigenvalues, exp (−iω_kt), where ω_k denotes the kth eigenvalue. Therefore, the transformation to the eigenbasis of the Hamiltonian implicitly happens during most spin dynamics simulations, meaning that, even if it was not set by the user, this is likely done by the linear algebra packages of the software. One may note that the basis does not affect the result of the calculation, but the choice of a more appropriate one may help rationalize the problem. In many cases, the initial choice is the Zeeman basis, in which spin operators are readily introduced based on Kronecker products of the Pauli matrices. Depending on the symmetry of the problem, it might be more convenient to change the basis to another one. As we will see in Sect. 4.1, a choice of coupled basis is preferable for understanding the zero-field J spectroscopy of coupled spins.

It is important to remark that Eq. (27) is only valid if the Hamiltonian is constant during the evolution period. The case where the Hamiltonian is time dependent is treated below. Note that the propagator is a unitary operator and therefore has the convenient property that its inverse is equal to its complex transpose (i.e., ${\widehat{U}}^{-\mathrm{1}}={\widehat{U}}^{\mathit{†}}$), which is much faster to compute than the matrix inverse ${\widehat{U}}^{-\mathrm{1}}$.

Equations (26) and (27) allow us to know the state of the system at any time t from the initial time t=0. To simulate the signal produced by the spin system during the course of the experiment, we must calculate the time domain signal at different time points. Note that in this case the Hamiltonian remains constant during free evolution. To calculate the signal at fixed time steps, it is convenient to first calculate the propagator $\widehat{U}\left(\mathrm{d}t\right)$ over period dt. We then apply Eq. (26) recursively to get the new density matrix $\widehat{\mathit{\rho }}\left(t_{k+\mathrm{1}}\right)$ from the previous one $\widehat{\mathit{\rho }}\left(t_k\right)$,

$$\begin{array}{}\text{(28)}& \widehat{\mathit{\rho }}\left(t_{k+\mathrm{1}}\right)=\widehat{U}\left(\mathrm{d}t\right)\widehat{\mathit{\rho }}\left(t_k\right){\widehat{U}}^{-\mathrm{1}}\left(\mathrm{d}t\right),\end{array}$$

where $t_{k+\mathrm{1}}-t_k=\mathrm{d}t$. To simulate ZULF spectra, we will also encounter situations where the Hamiltonian is time dependent. First, the Hamiltonian can vary with time but be “constant by block”. This is, for example, the case for the sudden field drop; the system is under a certain Zeeman Hamiltonian in the beginning of the experiment and suddenly under the ZULF Hamiltonian during detection. This situation does not present particular difficulties; the evolution of the system can be described step by step by both Eqs. (26) and (28).

Second, the Hamiltonian can vary continuously, as in the case of the adiabatic field drop, where the intensity of the magnetic field is ramped down to zero. This event can be simulated by propagating the evolution of the system during time intervals which are sufficiently short for the Hamiltonian to be considered constant during this time interval. The propagator must then be computed for each time increment. The form of the equation for propagation is similar to Eq. (28),

$$\begin{array}{}\text{(29)}& \widehat{\mathit{\rho }}\left(t_{k+\mathrm{1}}\right)=\widehat{U}\left(t_k\to t_{k+\mathrm{1}}\right)\widehat{\mathit{\rho }}\left(t_k\right){\widehat{U}}^{-\mathrm{1}}\left(t_k\to t_{k+\mathrm{1}}\right),\end{array}$$

where the propagator is given by

$$\begin{array}{}\text{(30)}& \widehat{U}\left(t_k\to t_{k+\mathrm{1}}\right)=\exp \left(-i\widehat{H}\left(t_k\right)\mathrm{d}t\right),\end{array}$$

where $\widehat{H}\left(t_k\right)$ is the Hamiltonian at time t_k. Note that the choice of $\widehat{H}\left(t_k\right)$ rather than $\widehat{H}\left(t_{k+\mathrm{1}}\right)$ in Eq. (30) is arbitrary, but in the limit of small intervals, the choice has no consequence.

2.6 Extract expectation values from the propagation

The propagation procedure described above gives access to the density matrix along time. To simulate the time domain signal, we need to extract a physical quantity from the density matrix as it evolves with time. The measured physical quantity of a ZULF experiment is the magnetic field produced by the nuclear spins of the sample at the location of an OPM. In a first approximation, we can consider that the whole sample is a point dipole interacting with the OPM and that this total dipole is the sum of the dipoles of the individual spin systems (Fig. 2 gives a visual representation of the approximation). Whether this approximation is appropriate or not depends on the geometry of the experimental setup. We have chosen the z axis as the quantization axis (defined by the detector, i.e., the OPM). Therefore, the physical quantity that we need to compute is the total magnetic field produced by the spins along the z axis at the location of the vapor cell:

$$\begin{array}{}\text{(31)}& \langle B_z\rangle ={\displaystyle \frac{{\mathit{\mu }}_{\mathrm{0}}}{\mathrm{2}\mathit{\pi }}}{\displaystyle \frac{\langle {\widehat{\mathit{\mu }}}_z^{\mathrm{tot}}\rangle }{r^{\mathrm{3}}}}={\displaystyle \frac{{\mathit{\mu }}_{\mathrm{0}}}{\mathrm{2}\mathit{\pi }}}{\displaystyle \frac{N\langle {\widehat{\mathit{\mu }}}_z\rangle }{r^{\mathrm{3}}}},\end{array}$$

where $\langle {\widehat{\mathit{\mu }}}_z^{\mathrm{tot}}\rangle$, μ₀, N, $\langle {\widehat{\mathit{\mu }}}_z\rangle$, and r are the magnetic moment of the sample along the z axis, the permeability of free space, the number of identical spin systems in the sample, their individual magnetic moments along the z axis, and the distance between the center of the sample and the center of the vapor cell, respectively.

Figure 2 Comparison of the real geometry of the sample of the OPM with the approximated one. The arrows represent local magnetization vectors parallel to the total magnetization vector.

For each identical spin system, we then compute the magnetic moment as the sum of the contributions of each spin l.

$$\begin{array}{}\text{(32)}& \langle B_z\rangle ={\displaystyle \frac{{\mathit{\mu }}_{\mathrm{0}}}{\mathrm{2}\mathit{\pi }}}{\displaystyle \frac{N}{r^{\mathrm{3}}}}\sum _{l=\mathrm{1}}^n\langle {\widehat{\mathit{\mu }}}_{lz}\rangle ={\displaystyle \frac{{\mathit{\mu }}_{\mathrm{0}}}{\mathrm{2}\mathit{\pi }}}{\displaystyle \frac{N\mathrm{\hslash }}{r^{\mathrm{3}}}}\sum _{l=\mathrm{1}}^n{\mathit{\gamma }}_l\langle {\widehat{I}}_{lz}\rangle ,\end{array}$$

where ${\widehat{\mathit{\mu }}}_{lz}$, γ_l, and ${\widehat{I}}_{lz}$ are the magnetic moment, the gyromagnetic ratio, and the angular momentum along the z axis of spin l, respectively. Note that n and N represent the number of spins in the molecule and the number of molecules in the sample, respectively. The notation 〈 〉 denotes the expectation value of a quantity. Particularly important ones are those that can be physically measured in the experiment. In the density matrix formalism that we are using, the expectation value of a physical quantity related to an operator $\widehat{A}$ is given by

$$\begin{array}{}\text{(33)}& \langle A\rangle =\mathrm{Tr}\left\{\widehat{A}\widehat{\mathit{\rho }}\right\},\end{array}$$

where Tr{ } denotes the matrix trace, i.e., the sum of all diagonal elements of the matrix representation of the operator. Note that the expectation value of ${\widehat{\mathit{\mu }}}_{lz}$ (or ${\widehat{I}}_{lz})$ is proportional to the polarization level of spin l which was accounted for in Eqs. (21) and (22). Therefore, the total magnetic moment calculated with Eq. (32) depends on the polarization of the different spin species.

If $\widehat{\mathit{\rho }}\left(t\right)$ is the density matrix at time t, we obtain the signal S(t) measured by the OPM by plugging Eq. (32) into Eq. (33)

$$\begin{array}{}\text{(34)}& S\left(t\right)=\langle B_z\rangle \left(t\right)={\displaystyle \frac{{\mathit{\mu }}_{\mathrm{0}}}{\mathrm{2}\mathit{\pi }}}{\displaystyle \frac{N\mathrm{\hslash }}{r^{\mathrm{3}}}}\mathrm{Tr}\left\{\widehat{O}\widehat{\mathit{\rho }}\left(t\right)\right\},\end{array}$$

where we have defined a “detection operator”,

$$\begin{array}{}\text{(35)}& \widehat{O}=\sum _{l=\mathrm{1}}^n{\mathit{\gamma }}_l{\widehat{I}}_{lz}.\end{array}$$

To obtain Eq. (34), we have used the fact that taking the trace of a matrix is a linear operation, and so the trace of a sum is the sum of the traces.

In the case of a sample with volume V=100 µL of ¹³C-formic acid prepolarized at 2 T at 298 K, with molar mass of 46 g mol⁻¹ and density of 1.22 g mL⁻¹, one finds that the amplitude of the magnetic field generated by the sample at a distance of r =1 cm is on the order of 10 pT using the above equations. This estimation does not take into account demagnetization effects caused by distribution of spins in space, giving the upper limit for the expected field. Experimentally measured magnetic fields are about 10 times smaller (Tayler et al., 2017).

2.7 Fourier transform the expectation values to obtain a spectrum

The time domain signal is what is measured by the ZULF NMR spectrometer. The final step of the simulation is to transform the measured signal from the time domain to the frequency domain using a discrete Fourier transform. Programming environments such as MATLAB or Mathematica are equipped with built-in functions for fast Fourier transformation. We will not discuss the mathematics behind this process, but we will give a few practical hints. Contrary to high-field NMR, ZULF spectra can be obtained with real magnetic field units (rather than arbitrary units). We will show how such units can be obtained.

Let us call t and S the arrays of numbers containing the time and corresponding time domain signal values, respectively, which resulted from the previous steps (note that, in MATLAB's programming environment, such arrays are usually called vectors). Let us call K the number of elements of both arrays (which corresponds to the number of points in the time domain signal). For now, S consists of a sum of oscillating signals which do not decay with time as our simulation did not include relaxation effects. If we perform a Fourier transform on S, we will obtain non-Lorentzian line shapes (with distinctive sinc patterns). We must therefore artificially include relaxation by multiplying the signal with an apodization function, to force the signal to decay to 0. For liquid-state signals, the most common choice is a monoexponential decay which can be expressed as

$$\begin{array}{}\text{(36)}& S_k^{\prime }=S_k\exp \left(-\mathit{\pi }{l}_{\mathrm{B}}{t}_k\right)=S_k\exp \left(-{\displaystyle \frac{t_k}{T_{\mathrm{2}}}}\right),\end{array}$$

where S_k, t_k, l_B, and T₂ are the kth elements of t and S, the line broadening (in Hz), and the coherence time constant (in s), respectively. Note that the coherence time constant is often referred to as the spin–spin relaxation constants or transverse relaxation time constant. The signal intensities $S_k^{\prime }$ define the apodized signal array S^′. As shown in Eq. (36), we may choose to express the apodization function using either the coherence time constant T₂ or the line broadening l_B (in Hz), which are related by $\mathit{\pi }{l}_{\mathrm{B}}=\mathrm{1}/T_{\mathrm{2}}$. The former is the time constant at which the time domain signal decays, while the latter is the full width at half height (FWHH) of the frequency domain signals. In order to avoid “truncating” the decay of the time domain signal and the related spectral artifacts, we must fulfill the condition T₂≪t_aq, where $t_{\mathrm{aq}}=max\left\{t_k\right\}$ is the acquisition time (or the length of the signal in the time domain). Typically, we may choose T₂ and t_aq so that t_aq=5T₂. Table 1 summarizes the parameters which were used in this paper.

The apodization function of Eq. (36) yields Lorentzian signals as one would expect. However, without further apodization, the baseline of the spectra will have some distortions (Zhu et al., 1993), with the main distortion being a small offset of the baseline. This problem arises because the time domain signal has its first point at time t=0, so the Fourier transform gives the integral of the first segment of twice larger amplitude than it should be. As proposed by Otting, this baseline offset can be removed by weighting the first and last point of the time domain signal by factor $\mathrm{1}/\mathrm{2}$ (Otting et al., 1986). However, because the integral of the Fourier transform is proportional to the first point of the time domain signal, this apodization does not preserve the integral. To obtain spectra without baseline offset and preserving the integral, we propose to use an apodization function which weights all points by 2 expect for the first and last ones:

$$\begin{array}{}\text{(37)}& S_k^{\prime \prime }=\left\{\begin{array}{c}{S}_k^{\prime }\\ \mathrm{2}{S}_k^{\prime }\end{array}\begin{array}{c}\mathrm{if}k=\mathrm{1}\mathrm{or}L,\\ \mathrm{otherwise},\end{array}\right.\end{array}$$

where L is the number of points in the Fourier transform (see below). We show in the Supplement that this apodization function preserves the integral (see Supplement S2).

In MATLAB programming language, the function for fast Fourier transformation fft() takes array S^′′ as input and returns the frequency domain array which corresponds to the simulated spectrum. Optionally, one may add a second argument L to fft() to include zero-filling in the Fourier transform. Including zero-filling has the advantage of increasing the number of points per FWHH on the spectrum without increasing the computation time of the propagation. Due to MATLAB's Fourier transform convention, it is convenient to retransform the signal with fftshift() in order to obtain a Fourier transformed signal with 0 as the middle frequency. We then divide the output of MATLAB's Fourier transform by the number of points L:

$$\begin{array}{}\text{(38)}& S\left(\mathit{\nu }\right)={\displaystyle \frac{\mathrm{1}}{L}}FS\left(t\right),\end{array}$$

where F designates the Fourier transform. The frequency domain signal obtained after this whole procedure has units of magnetic field (e.g., pT). Changing the zero-filling L changes the intensity of the frequency domain signal but preserves the integrals.

MATLAB's fft() function does not generate the frequency array associated with the Fourier transformed signal. The frequency array ν (in Hz) can be generated based on the following expression:

$$\begin{array}{}\text{(39)}& {\mathit{\nu }}_k={\displaystyle \frac{k}{L}}f,k\in \left[\left[-{\displaystyle \frac{L}{\mathrm{2}}};{\displaystyle \frac{L}{\mathrm{2}}}-\mathrm{1}\right]\right],\end{array}$$

where the sampling frequency (in Hz) is given by

$$\begin{array}{}\text{(40)}& f={\displaystyle \frac{K-\mathrm{1}}{t_{\mathrm{aq}}}}.\end{array}$$

Figure 3 Illustration of signal sampling and the effect of undersampling. Panels (a), (c), and (e) represent a cosine oscillating at 1 Hz in gray sampled with various frequencies f (1.3, 3.7, and 7.3 Hz, respectively). The blue dots represent the samples. In each case, the Fourier transform is shown in panels (b), (d), and (f), respectively. When the sampling frequency is lower than 1 Hz, the peak cannot appear at 1 Hz and is therefore found at a fictitious position.

The sampling frequency of the time domain signal gives the maximum frequency that can be appropriately sampled. Figure 3 illustrates the consequence of choosing a sampling frequency which is lower than the maximum frequency. If the sampling frequency is lower than the signal to be sampled, the Fourier transformed signal lies outside the spectral width (between $-f/\mathrm{2}$ and $+f/\mathrm{2}$). However, due to the “refolding effect” of the Fourier transform, the signal still appears in the spectrum but at irrelevant positions. To avoid this, one may repeat the simulation by increasing the sampling frequency and keeping other parameters constant. If the sampling frequency is sufficient, the spectrum should not be affected.

The choice of the parameters discussed in this section and above influences the outcome of the simulation in the same way as it does for the experiment. Once an NMR simulation is running, one might want to play with combinations of f, K, t_aq, T₂, and L until the simulated spectra display convenient features. If one intends to simulate spectra to match experimental data, one might simply perform the simulation with the same f, K, L, and t_aq values. Table 1 summarizes the parameters which were used in this paper.

Table 1 List of parameters that were used to simulate the time domain signals and spectra in Fig. 4.

The procedure described here yields an NMR signal which is symmetric around 0. As a consequence, each signal is found both in the positive and negative frequencies, and the integral is split into the two duplicates. Because the experimental procedure that we are simulating does not differentiate negative and positive frequencies, we discard the frequency domain signal corresponding to negative frequencies and multiply the abscissa of the frequency domain signal corresponding to positive frequencies by a factor of 2. This operation corresponds to “folding” the spectrum around ν=0. Note that in high-field NMR, the measured signal is complex and is therefore not split into positive and negative halves. The central frequency of the spectrum at high field is given by the carrier frequency of the spectrometer (e.g., typically 400 MHz for ¹H at 9.4 T). Section 2.8 describes this difference between high-field and ZULF NMR in more detail.

Whether the time domain signal which results from the simulation is real or complex, the Fourier transform yields a complex frequency domain signal. To get a spectrum consisting of a signal intensity as a function of the frequency, we must use the real part of the frequency domain signal. Depending on the experiment that we are simulating, we might find that some or all spectral components of the frequency domain signal are not in phase. To compensate for this, one might apply a phase correction by multiplying each point of the frequency domain signal by a complex constant exp iϕ, where φ is the phase correction before taking its real part.

$$\begin{array}{}\text{(41)}& I_{\mathrm{r}}\left(\mathit{\nu }\right)=\mathrm{Re}\left\{I_{\mathrm{c}}\left(\mathit{\nu }\right)\exp i\mathit{\varphi }\right\},\end{array}$$

where I_r(ν) and I_c(ν) are the real and complex frequency domain signals, respectively.

In summary, the Fourier transform procedure that we have described has the following steps:

1. Apply a monoexponential apodization window to the time domain signal so that it decays to 0 (see Eq. 36).

2. Apply the apodization described by Eq. (37) to avoid baseline artifacts in the frequency domain signal.

3. Obtain the complex frequency domain signal by Fourier transforming the time domain signal using a fast Fourier transform algorithm.

4. Generate the corresponding frequency axis using Eqs. (39) and (40).

5. Remove the negative frequencies from both the frequency axis and the frequency domain signal and multiply the abscissa of the frequency domain signal by 2 to account for the partition of the signal integral between positive and negative frequencies.

6. Take the real part of the signal after applying an optional phase correction (see Eq. 41).

2.8 Comparison with high-field NMR

We conclude this theory section by listing the main differences between high-field and ZULF NMR, which are summarized in Table 2. As is the case for the rest of the paper, our description is limited to small molecules containing spin $\mathrm{1}/\mathrm{2}$ in the liquid state.

At high magnetic field, the Zeeman interaction is much larger than the J coupling and therefore dominates the dynamics. Furthermore, the Larmor frequency of the spins (which results from the Zeeman interaction) is slightly shifted by the presence of the electron cloud around the nuclei. This phenomenon, called the chemical shift, gives a slightly different Larmor frequency for nuclei in different positions in a molecule, which spreads over typically 10 and 200 ppm around the Larmor frequency for ¹H and ¹³C spins, respectively. At ZULF, the J coupling dominates while the Zeeman interaction is a perturbation, and the chemical shift plays no role (in that it is a small perturbation of a small perturbation).

In Fig. 1 and in the simulations presented in this paper, we have assumed that the detector was positioned below the sample (along the z axis in our axis convention) and that it was sensitive to magnetic field along the z axis. Although this choice is typical, it is not the only possibility. In common high-field experiments, the oscillating signal emitted by the spins is recorded perpendicular to the static magnetic field. Detection at ZULF is performed with magnetometers that are sensitive to the total magnetic field produced by the sample. The operator corresponding to this observable is the sum of the magnetic moment of the spins along the sensitive axis of the OPM (see Eq. 34). In typical experiments, a single detector is used, which results in a real signal. Note that an imaginary ZULF signal could be obtained if the OPM were to have several sensitive axes or more than one detector were used. High-field NMR uses Faraday induction in pickup coils. Signals originating from different nuclei are usually not observed in the same experiment as their Larmor frequencies are too far apart, and the NMR coils are only sensitive over a limited bandwidth. The operator corresponding to inductive detection in pulsed NMR is non-Hermitean and therefore yields complex signals. An extra step of the acquisition process at high field that is not required at ZULF is modulating the signal recorded by the coil with a carrier frequency. Indeed, the NMR coil picks up a signal at the Larmor frequency of the spins, which is too high to be digitized (e.g., 400 MHz for ¹H at 9.4 T). Instead, the signal is mixed with a carrier frequency, and only the difference is digitized, over a small bandwidth (e.g., over 10 ppm, corresponding to 4 kHz for ¹H spins at 9.4 T). The signals at ZULF can be detected without mixing the frequency as the they are sufficiently low to be digitalized directly. For more details on the signal modulation at high field, the reader is referred to chapter 4 of James Keeler's book Understanding NMR spectroscopy (Keeler, 2010).

The code in Supplement S3 presents in great detail the simulation of the spectra for a pair of J-coupled ¹H and ¹³C spin pairs at ZF and ULF and at high field for both ¹H and ¹³C (9.4 T; Stern and Sheberstov, 2023). The code is decomposed in sections corresponding to Sect. 2.1 to 2.7 of the text above, and, whenever possible, the equations presented in this paper are referenced in the code. The reader is encouraged to open this code to understand the difference between simulating a spectrum at high field and at ZULF. The code can be opened in PDF, including the figures, for those who do not have a MATLAB license.

Table 2 Comparison between high-field and ZULF NMR for typical experiments. Note that quadrature detection (and thus imaginary signals) is possible at ZULF, although uncommon. SQUID stands for superconducting quantum interference device.

3 Results of numerical simulations

3.1 Excitation schemes on an XA spin system

The ZF and ULF spectra of an XA spin system with a J coupling of 140 Hz were simulated for different experimental sequences, assuming that the sample consists of 100 µL of solution where the spin system has a concentration of 27 mol L⁻¹. The code and its PDF version are presented in Supplement S4. Figure 4 shows the experimental sequences, as well as the simulated time domain and frequency domain signals. For all simulations, the sample was assumed to have spent sufficient time in a prepolarizing field of 2 T at 298 K to be at Boltzmann's equilibrium. The polarizations of the ¹³C and ¹H spins were calculated using Boltzmann's distribution (see Eq. 19) and used to compute the single-spin density matrices of the ¹³C and ¹H spins, ${\widehat{\mathit{\rho }}}_{\mathrm{eq}}{(}^{\mathrm{13}}\mathrm{C})$ and ${\widehat{\mathit{\rho }}}_{\mathrm{eq}}{(}^{\mathrm{1}}\mathrm{H})$ (see Eq. 21). The density matrix of the two-spin system was computed by taking the Kronecker product of the single-spin density matrices ${\widehat{\mathit{\rho }}}_{\mathrm{0}}={\widehat{\mathit{\rho }}}_{\mathrm{eq}}{(}^{\mathrm{13}}\mathrm{C})\otimes {\widehat{\mathit{\rho }}}_{\mathrm{eq}}{(}^{\mathrm{1}}\mathrm{H})$ (see Eq. 23). The identity was removed from the two-spin density matrix using Eq. (24). The resulting density matrix was assumed to represent the initial state of the simulation (as explained above, only the Zeeman terms are considered to contribute to the initial state). For each experimental sequence, the spectrum was simulated both at 0 nT (including only the J Hamiltonian ${\widehat{H}}_J$; see Eq. 12) and with a basis field of 0.5 µT along the x axis, that is, orthogonal to both the direction of the prepolarizing field and the sensitive axis (including both the J Hamiltonian ${\widehat{H}}_J$ and the Zeeman Hamiltonian ${\widehat{H}}_{\mathrm{Z}}$; see Eqs. 12 and 10). The time domain signal was computed by propagating the density matrix under the effect of the Hamiltonian for a total time of 5 s (parameter t_aq), discretized into 4096 points (parameter K), and corresponding to time intervals dt of 1.2207 ms (parameter τ_D). Prior to the propagation loop, the ZF and ULF propagators $\widehat{U}$ for this particular time step dt (see Eq. 30) and the observable operator $\widehat{O}$ (see Eq. 35) were computed only once.

Figure 4 Excitation schemes for an XA spin system corresponding to ¹³C and ¹H spins with a J coupling of 140 Hz and corresponding simulated time domain signals and spectra. The vertical dashed line indicates the J coupling. The time domain signal was computed by propagating the density matrix under the effect of the Hamiltonian for a total time of 5 s (parameter t_aq), discretized into 4096 points (parameter K), and corresponding to time intervals dt of 1.2207 ms (parameter τ_D). A monoexponential apodization function was applied to the time domain signal, with a coherence time constant T₂ of 1 s. The apodized time domain signal was Fourier transformed with zero-filling to 65 536 points.

The density matrix was propagated from time t_k to time $t_{k+\mathrm{1}}=t_k+\mathrm{d}t$ under the Hamiltonian (ZF or ULF) using the formula ${\widehat{\mathit{\rho }}}_{k+\mathrm{1}}=\widehat{U}{\widehat{\mathit{\rho }}}_k{\widehat{U}}^{-\mathrm{1}}=\widehat{U}{\widehat{\mathit{\rho }}}_k{\widehat{U}}^{\mathit{†}}$ (see Eq. 29). At each time point k of the propagation (realized by a “for” loop), the signal intensity of the time domain signal was extracted from the density matrix using the trace $\mathrm{Tr}\left\{\widehat{O}{\widehat{\mathit{\rho }}}_k\right\}$ (see Eq. 33) (in pT). In theory, the trace of a Hermitean operator should be real. However, due to the finite machine precision of the numeric algorithm, the trace can sometimes contain a nonzero imaginary part. This residual imaginary part is discarded by taking the real part of the trace $\mathrm{Re}\left(\mathrm{Tr}\left\{\widehat{O}{\widehat{\mathit{\rho }}}_k\right\}\right)$. This point might appear secondary, but dealing with complex numbers while thinking they are real can lead to mistakes. After propagation, a monoexponential apodization function was applied to the time domain signal (see Eq. 36), with a coherence time constant T₂ of 1 s. A second apodization function was applied to avoid baseline artifacts (see Eq. 37). The apodized time domain signal was Fourier transformed with zero-filling to 65 536 points, using MATLAB's built-in functions. The real part of the Fourier transform is shown in Fig. 4. The frequency axis of the spectra was computed using Eqs. (39) and (40). The spectra are symmetric around zero, and so it is common to work only with the positive frequencies as shown in Fig. 4.

Simulating the sudden field drop experiment is the simplest case presented here. Because the coherence excitation scheme (or mixing) only consists of bringing the spin from high magnetic field to ZF or ULF, the simulation only consists of propagating the high-field thermal equilibrium density matrix under the ZF or ULF Hamiltonian. The ZF spectrum consists of one line at the J coupling and one at zero frequency (see Fig. 4a). Including a bias field of 0.5 µT along the x axis (ULF case) splits the J peak as well as the line at zero frequency.

The simulations presented in Fig. 4b–d feature an adiabatic field drop. We used a monoexponential field drop from B_start=200 µT to 0 occurring over t_decay=0.5 s with a decay time constant of τ=0.05 s, described by

$$\begin{array}{}\text{(42)}& B\left(t\right)=B_{\mathrm{start}}{\displaystyle \frac{\exp \left(-\frac{t}{\mathit{\tau }}\right)-\exp \left(-\frac{t_{\mathrm{decay}}}{\mathit{\tau }}\right)}{\mathrm{1}-\exp \left(-\frac{t_{\mathrm{decay}}}{\mathit{\tau }}\right)}},\end{array}$$

which fulfills the conditions B(0)=B_start and B(t_decay)=0. During the field drop, the Hamiltonian $\widehat{H}\left(t\right)={\widehat{H}}_J+{\widehat{H}}_{\mathrm{Z}}\left(t\right)$ is time dependent. This step thus cannot be simulated in a single propagation step. Instead, it must be discretized into substeps dt that are sufficiently short for the Hamiltonian to be considered time independent. Here, the 0.5 s time interval was discretized into 5000 steps of 0.1 ms. At time t=0, the density matrix is the thermal equilibrium density matrix ${\widehat{\mathit{\rho }}}_{\mathrm{eq}}$ obtained above. At each time step t_k, the propagator $\widehat{U}\left(t_k\to t_{k+\mathrm{1}}\right)=\exp \left(-i\widehat{H}\left(t_k\right)\mathrm{d}t\right)$ is computed (see Eq. 30), and the density matrix is propagated from time t_k to time $t_{k+\mathrm{1}}=t_k+\mathrm{d}t$ (see Eq. 29) under the Hamiltonian $\widehat{H}\left(t_k\right)={\widehat{H}}_J+{\widehat{H}}_{\mathrm{Z}}\left(t_k\right)$. We name ${\widehat{\mathit{\rho }}}_{\mathrm{adia}}$ the density matrix obtained after this process. A question arises here: is this magnetic field drop that we have chosen sufficiently slow to be considered adiabatic? In other words, is ${\widehat{\mathit{\rho }}}_{\mathrm{adia}}$ stationary? A simple way to ensure that it is the case is to simulate the spectrum at ZF after the magnetic field drop without any excitation pulse, that is, taking ${\widehat{\mathit{\rho }}}_{\mathrm{adia}}$ as the density matrix at time t=0, ${\widehat{\mathit{\rho }}}_{\mathrm{0}}$. If the transition is adiabatic, then the system should remain stationary; that is, the time domain signal should feature no oscillation and the spectrum no peak. Figure 4b shows the result of this procedure, which confirms that the transition is adiabatic. The only feature of the ZF spectrum in Fig. 4b is the line at zero frequency. This line originates from the non-oscillating magnetization decaying with T₂, which is the result of the apodization function that we have applied. Verifying that the ZF spectrum is flat also ensures that the field drop was discretized into sufficiently short time intervals dt.

The density matrix after the adiabatic field drop ${\widehat{\mathit{\rho }}}_{\mathrm{adia}}$ obtained above was used for the simulations presented in Fig. 4c–d. In the experimental sequences of Fig. 4c–d, the adiabatic field drop is followed by a magnetic field pulse either along the z or x axis. This was simulated by propagating ${\widehat{\mathit{\rho }}}_{\mathrm{adia}}$ under the pulse Hamiltonian to obtain ${\widehat{\mathit{\rho }}}_{\mathrm{0}}={\widehat{U}}_{\mathrm{p}}\left({\mathit{\tau }}_{\mathrm{p}}\right){\widehat{\mathit{\rho }}}_{\mathrm{adia}}{\widehat{U}}_{\mathrm{p}}^{\mathit{†}}\left({\mathit{\tau }}_{\mathrm{p}}\right)$, where ${\widehat{U}}_{\mathrm{p}}\left({\mathit{\tau }}_{\mathrm{p}}\right)$ is the propagator of the pulse Hamiltonian ${\widehat{H}}_{\mathrm{p}}={\widehat{H}}_J+{\widehat{H}}_{\mathrm{Z}}$, which acts on the density matrix during pulse length τ_p. The Zeeman Hamiltonian depends on the magnetic field intensity of the pulse B_p and its direction (see Eq. 15). For the z-axis pulse, we used a pulse intensity and length of 50 µT and 150 µs, respectively. For the x-axis pulse, we used a pulse intensity and length of 50 µT and 910 µs, respectively. These choices are justified in the next section. The resulting density matrices ${\widehat{\mathit{\rho }}}_{\mathrm{0}}$ were used as the density matrix at time t=0 of the time domain signal, which was computed and Fourier transformed as described above. In the case of the z-axis pulse experiment, the peaks of interest (J peak at 140 Hz) were found to be out of phase; a phase correction e^iϕ with $\mathit{\varphi }=\mathit{\pi }/\mathrm{2}$ was thus applied to the Fourier transform. Adjusting the phase for the J peak caused the lower-frequency peaks to be out of phase. Interestingly, in Fig. 4d, the intensity of the J peak is higher than for the other excitation schemes while the lower-frequency peaks are suppressed, indicating that all the available polarization has been transferred to the J peak.

3.2 Rabi oscillation curves

The pairs of magnetic field intensity and length of the pulses used for the simulation in Fig. 4d were chosen by simulating Rabi curves for both the z- and x-axis pulses. The high-field NMR equivalent to the Rabi curve is the “nutation experiment”, which consists of recording a series of NMR detections while keeping the RF pulse power constant and varying the pulse length (or the pulse length is kept constant and the amplitude is varied; Tayler et al., 2017). The nutation or Rabi curve is the plot of the signal intensity as a function of the varied parameter. It allows one to determine the pair of RF power and pulse length which maximizes the signal intensity. Except in the presence of rapid relaxation effects or RF field inhomogeneities, the observed curve is sinusoidal. At ZULF, the Rabi curve is more complex and depends on the spin system under scrutiny. To simulate the Rabi curve at ZF, we repeated the simulation of the ZF spectra for an experiment with an adiabatic field drop (using the same parameters as above) followed by a pulse of 50 µT along the z and x axes, varying the pulse length from 0 to 3000 µs (the code and its PDF version are presented in Supplement S4). The time domain signal was Fourier transformed as described above, and the frequency domain signal was integrated from 138 to 142 Hz. The signal integral of the J peak is plotted as a function of the pulse length in Fig. 5. The signal integral of the sudden drop experiment is shown as a horizontal dashed line for comparison. When a pulse along the z axis is used, a simple sinusoidal curve is obtained, and its maximum matches that of the sudden drop experiment (see Fig. 5a). The first maximum is reached for a pulse length of 150 µs. When a pulse along the x axis is used, a more complex pattern is obtained, and the maximum is found to be 1.64 times higher than the sudden drop experiment (see Fig. 5b). The first global maximum is reached for a pulse length of 910 µs.

Figure 5 Rabi curves at ZF with excitation pulses along z (a) and x (b) axes applied to an XA spin system. The horizontal dashed line represents the signal integral of the sudden drop ZF experiment. The time domain signal was computed by propagating the density matrix under the effect of the Hamiltonian for a total time of 5 s (parameter t_aq), discretized into 4096 points (parameter K), and corresponding to time intervals dt of 1.2207 ms (parameter τ_D). A monoexponential apodization function was applied to the time domain signal, with a coherence time constant T₂ of 1 s. The apodized time domain signal was Fourier transformed with zero-filling to 65 536 points. The frequency domain signal was then integrated from 138 to 142 Hz. The Rabi curve represents the integral compared with the excitation pulse length.

3.3 XA_n spin system

The simulations shown up to this point only deal with an XA spin system, which typically corresponds to ¹³C-formate (or ¹³C-formic acid), where the ¹³C spin interacts with a single ¹H through a J coupling of 195–222 Hz (Blanchard and Budker, 2016; Tayler et al., 2017) (depending on experimental conditions). ¹³C,¹⁵N-cyanide groups are also interesting two-spin systems which were used in ZULF experiments (Blanchard et al., 2020, 2015). We now extend the simulation to incorporate multiple A spins. An XA₂ spin system is, for example, met in ¹³C-glycine (Put et al., 2021). XA₃ spins are met in a number of molecules containing methyl groups such as ¹³C-pyruvate (Barskiy et al., 2019). XA₄ (for example ¹⁵N-ammonium cation; Barskiy et al., 2019) and XA₅ are less common, but they are presented here to show the pattern that arises when adding spins.

Figure 6 shows the simulated spectra for sudden drop experiments with detection at ZF and ULF of XA_n spin systems with $n=\mathrm{1},\mathrm{2},\mathrm{\dots },\mathrm{5}$, where X represents a ¹³C spin, and A_n represents ¹H spins with a J coupling of 140 Hz between X and A spins and 10 Hz among A spins (the code and its PDF version are presented in Supplement S5). All the relevant mathematics to construct the operators of an $m=n+\mathrm{1}$ spin system are given in the Theory section. For an XA₅ spin system, the Hilbert space has 2⁶=64 dimensions (and related operators). To avoid constructing each operator manually, recursive formulae were used (see Eqs. 13 and 23). The time domain signal was computed by propagating the density matrix under the effect of the Hamiltonian for a total time of 5 s (parameter t_aq), discretized into 8192 points (parameter K), and corresponding to time intervals dt of 0.6104 ms (parameter τ_D). A monoexponential apodization function was applied to the time domain signal, with a coherence time constant T₂ of 1 s. The apodized time domain signal was Fourier transformed with zero-filling to 32 768 points.

Figure 6 Simulation of ZF and ULF spectra after sudden field drop for XA, XA₂, XA₃, XA₄, and XA₅ spin systems with a J coupling of 140 Hz between X and A spins and 10 Hz among A spins. The time domain signal was computed by propagating the density matrix under the effect of the Hamiltonian for a total time of 5 s (parameter t_aq), discretized into 4096 points (parameter K), and corresponding to time intervals dt of 1.2207 ms (parameter τ_D). A monoexponential apodization function was applied to the time domain signal, with a coherence time constant T₂ of 1 s. The apodized time domain signal was Fourier transformed with zero-filling to 32 768 points.

Increasing the number of A spins increases the number of spectral components in the spectrum. A known result of ZULF NMR appears in this simulation: for odd numbers of n, the ZF spectrum features lines at integer multiples of the J coupling k⋅J_AX with $k\in \left[\right[\mathrm{1};(n+\mathrm{1})/\mathrm{2}\left]\right]$, while for even numbers of n, it features lines at half-integer multiples of the J coupling $k\cdot J_{\mathrm{AX}}/\mathrm{2}$ with $k\in \left[\right[\mathrm{1};n/\mathrm{2}\left]\right]$. Adding a 0.5 µT bias field along the x axis during detection (that is, performing ULF detection) splits the J lines. The higher the multiple of the J line, the greater the number of splittings. Note that the intensity of NMR signals at high field increases upon adding more equivalent spins to the spin system. The analysis of Fig. 6 shows that this logic does not apply to the J lines for the ZULF case, where the spectrum completely changes upon changes in the spin topology. For example, note that the amplitude of the J line for the XA system has the same intensity as the J line for the XA₂ system (appearing at $\mathrm{3}/\mathrm{2}\cdot J_{\mathrm{AX}}$ frequency). Likewise, the two J lines for the XA₃ system have the same total intensity as the two J lines for the XA₄ system. An empirical law of conservation of the total spectral intensity for the J lines can be deduced by looking at Fig. 6: indeed, the total intensity of all J lines is the same for any XA_n system, assuming equal sample volume, prepolarization, etc. On the other hand, the intensity of low-frequency peaks shown in Fig. 6 is proportional to the total number of spins in the spin system, like in high-field NMR. This is of course expected as these signals are associated with the precession of total magnetization around residual ULF field, and total magnetization is proportional to the number of spins.

4 Interpretation

We are now going to show how to calculate ZULF NMR spectra considering energy levels and transition probabilities rather than through the numerical propagation of the density matrix. We will derive analytical solutions for the XA_n system, but the same approach can be used for more complex spin systems. This approach was investigated in the following references: Butler et al., 2013a; Theis et al., 2013; and Emondts et al., 2014. Here we aim to present it with more explanations and explicit derivations, but we limit ourselves to only the simplest spin systems.

The relative contributions of ${\widehat{H}}_{\mathrm{Z}}$ (see Eq. 10) and ${\widehat{H}}_J$ (see Eq. 12) terms depends on the magnetic field strength. In the high-field extreme, for a heteronuclear spin system, ${\widehat{H}}_{\mathrm{Z}}$ is the dominant term, and ${\widehat{H}}_J$ is considered as a first-order perturbation. In this case, heteronuclei are said to be weakly coupled, and their eigenstates coincide with the Zeeman states (e.g., those in Eq. 9). At zero field, the weak coupling approximation is not valid; the Zeeman states do not correspond to the eigenstates of system. However, it is still possible to calculate analytically the eigenstates for some spin systems, and the simplest case is when all the spins are identical (A_n system). In this case, the Hamiltonian is represented by only the ${\widehat{H}}_J$ term, and it commutes with the square of the total angular momentum operator.

$$\begin{array}{}\text{(43)}& \begin{array}{rl}& {\widehat{\mathit{F}}}^{\mathrm{2}}={\widehat{F}}_x^{\mathrm{2}}+{\widehat{F}}_y^{\mathrm{2}}+{\widehat{F}}_z^{\mathrm{2}},\\ & {\widehat{F}}_{\mathit{\mu }}=\sum _{l=\mathrm{1}}^n{\widehat{I}}_{l\mathit{\mu }};\mathit{\mu }\in \left\{x,y,z\right\},\\ & \left[{\widehat{H}}_J,{\widehat{\mathit{F}}}^{\mathrm{2}}\right]={\widehat{H}}_J{\widehat{\mathit{F}}}^{\mathrm{2}}-{\widehat{\mathit{F}}}^{\mathrm{2}}{\widehat{H}}_J=\mathrm{0},\end{array}\end{array}$$

where n is the number of spins in the system. It is well known that any pair of commuting Hermitean operators share their eigenspaces (Levitt, 2013). The set of eigenstates which forms an eigenbasis for both operators simultaneously is unique in cases where there are no degeneracies (i.e., all the eigenvalues for both operators are different). When there are degeneracies, the common eigenbasis is not unique. It turns out that ${\widehat{H}}_J$ and ${\widehat{\mathit{F}}}^{\mathrm{2}}$ operators do have degeneracies, and this results in the existence of an infinite number of different shared eigenbases. Let us describe how to find such a set of eigenstates.

4.1 Eigenstates at zero field

The eigenstates of a ${\widehat{\mathit{F}}}^{\mathrm{2}}$ operator can be expressed in terms of the total spin and its projection quantum numbers. The conventional way to express them is to use the |F,m_F〉 notation, where F denotes the total spin, and m_F denotes the projection onto a quantization axis ($m_F\in \left\{-F,-F+\mathrm{1},\mathrm{\dots },F-\mathrm{1},F\right\})$. For example, by definition, for a single spin $\mathrm{1}/\mathrm{2}$, we have the sates $\left|\mathit{\alpha }〉\right.\equiv \left|\mathrm{1}/\mathrm{2},\mathrm{1}/\mathrm{2}〉\mathrm{and}\right.\left|\mathit{\beta }〉\right.\equiv \left|\mathrm{1}/\mathrm{2},-\mathrm{1}/\mathrm{2}〉\right.$. For a pair of spins, we have the three triplet states $\left|T_{+\mathrm{1}}〉\right.\equiv \left|\mathrm{1},\mathrm{1}〉;\right.\left|T_{\mathrm{0}}〉\right.\equiv \left|\mathrm{1},\mathrm{0}〉;\mathrm{and}\right.\left|T_{-\mathrm{1}}〉\right.\equiv \left|\mathrm{1},-\mathrm{1}〉\right.$; and the singlet state $\left|S_{\mathrm{0}}〉\right.\equiv \left|\mathrm{0},\mathrm{0}〉\right.$. Any |F,m_F〉 state is an eigenstate of the ${\widehat{\mathit{F}}}^{\mathrm{2}}$ and ${\widehat{F}}_z$ operators with the following eigenvalues:

$$\begin{array}{}\text{(44)}& \begin{array}{rl}& {\widehat{\mathit{F}}}^{\mathrm{2}}\left|\right.F,m_F〉=F(F+\mathrm{1})\left|\right.F,m_F〉,\\ & {\widehat{F}}_z\left|\right.F,m_F〉=m_F\left|\right.F,m_F〉.\end{array}\end{array}$$

To find the total spin of a system constituted by n spins, one must sum up the angular momenta of the individual spins, which is a common procedure in the field of atomic physics but not so much in NMR. All possible values of the angular momentum of the interacting spins are added up to constitute a set of uncoupled quasiparticles with different total spin. The total spin F of a system constituted by two spins I and S can take the values with steps of 1 between the sum I+S and the absolute value of their difference:

$$\begin{array}{}\text{(45)}& \left|I-S\right|\le F\le I+S.\end{array}$$

For a pair of spins $\mathrm{1}/\mathrm{2}$, the possible values are $F=\mathrm{0},\mathrm{1}$. For n spins, the summation should proceed until all the possible pairs of the angular momenta of the individual spins are summed up. As an illustration, consider a coupled system of three spins $\mathrm{1}/\mathrm{2}$ (see Fig. 7). First, any two spins are added up together to give F=1 (a triplet) and F=0 (a singlet). Then, the remaining spin $\mathrm{1}/\mathrm{2}$ is added up to the quasiparticles formed in the previous step (spins 1 and 0 in this case). As a result, the initial A₃ system is decomposed into three subsystems with total spins of $F=\mathrm{3}/\mathrm{2}$, $\mathrm{1}/\mathrm{2}$ (addition of 1 and $\mathrm{1}/\mathrm{2}$), and $\mathrm{1}/\mathrm{2}$ (addition of 0 and $\mathrm{1}/\mathrm{2}$).

Figure 7 Procedure for adding up the angular momenta for the A₃ spin system.

A useful property of such a decomposition can be illustrated at this point: the total spin operator commutes with all rotation operators (e.g., $\exp (-i\mathit{\theta }{\widehat{I}}_z)$); therefore, 3D rotations will never mix terms of the wave function belonging to different total spin, e.g., spin $\mathrm{3}/\mathrm{2}$ with $\mathrm{1}/\mathrm{2}$. At ZF, there is no distinction between directions; therefore, the eigenstates must be invariant with respect to 3D rotations. This also partially explains the existence of an infinite number of eigenbases for ${\widehat{\mathit{F}}}^{\mathrm{2}}$, as all different orientations of the $\mathit{\{}x,y,z\mathit{\}}$ system correspond to different bases.

One can check that the total number of the spin states remains the same after the procedure of adding up the spins. On the one hand, the number of states formed by n coupled spins I equals to (2I+1)ⁿ, which is 8 in the considered case. On the other hand, a manifold with a total spin F has 2F+1 different states associated with different possible projections of the spin on the quantization axis. Therefore, there are $\mathrm{4}+\mathrm{2}+\mathrm{2}$ states in the considered case.

The explicit form of the resulting eigenstates can be obtained in terms of “uncoupled” spin states, which are constructed as a Kronecker product of the individual Zeeman states (see Eq. 9). The resulting state |F,m_F〉 of the addition of two angular momenta (I and S) can be represented as the following linear combination:

$$\begin{array}{}\text{(46)}& \left|F,m_F\right.〉=\sum _{m_I,m_S}{C}_{I,m_I,S,m_S}^{F,m_F}\left|I,m_I,S,m_S\right.〉,\end{array}$$

where $C_{I,m_I,J,m_J}^{F,m_F}$ represents Clebsch–Gordan coefficients, which are defined by

$$\begin{array}{}\text{(47)}& C_{I,m_I,S,m_S}^{F,m_F}=〈I,m_I;S,m_S\left|\right.F,m_F〉.\end{array}$$

Each Clebsch–Gordan coefficient is specified by six numbers: the total spin of the coupled state F, its projection m_F, and the total spins of the uncoupled states and their projections (I, S, m_I, m_S). Coefficient $C_{I,m_I,S,m_S}^{F,m_F}$ represents “how much” of uncoupled state $\left|I,m_I,S,m_S\right.〉$ there is in a coupled state |F,m_F〉. The analytical values of the Clebsch–Gordan coefficients can be calculated using recursive expressions, which are available in many software packages and textbooks. Table S1.1 in the Supplement provides the relation between the coupled and uncoupled states for the considered A₃ system and shows explicitly how to calculate them. The full set of all possible |F,m_F〉 states forms the new basis that is better suited than the Zeeman basis for ZULF NMR. In fact, this basis coincides with the eigenstates at ZULF for A_n and for XA_n systems, but this basis is also a good starting point for more complicated cases. We will refer to this new basis as the “coupled” basis, because it is appropriate for the description of strongly coupled spins.

4.2 Eigenenergies at zero field

Having the eigenstates, we can now proceed with finding the eigenvalues of the Hamiltonian; these values correspond to the energy of the states and therefore determine the frequencies of ZULF NMR transitions. It turns out that A_n systems are not detectable at ZULF; it is shown in the next section (where intensities of transitions are calculated) that they give rise to no observable transition. At least two types of nuclei with different gyromagnetic ratios are necessary for an observable transition to exist. Therefore, we consider an XA_n system from now on. We will denote the operators associated with the X spin as $\widehat{\mathit{S}}$ and with A spins as $\widehat{\mathit{I}}$. It is also convenient to introduce total spin operators for A spins: ${\widehat{F}}_{\mathrm{A}\mathit{\mu }}=\sum _{l=\mathrm{1}}^n{\widehat{I}}_{l\mathit{\mu }}$, $\mathit{\mu }\in \mathit{\{}x,y,z\mathit{\}}$. The Hamiltonian at ZF for this spin system is given by

$$\begin{array}{}\text{(48)}& {\widehat{H}}_{\mathrm{AX}}=\mathrm{2}\mathit{\pi }{J}_{\mathrm{AX}}\sum _{l=\mathrm{1}}^n\widehat{\mathit{S}}\cdot {\widehat{\mathit{I}}}_l+\mathrm{2}\mathit{\pi }{J}_{\mathrm{AA}}\sum _{l=\mathrm{1}}^{n-\mathrm{1}}\sum _{k>l}^n{\widehat{\mathit{I}}}_l\cdot {\widehat{\mathit{I}}}_k.\end{array}$$

The ${\widehat{H}}_{\mathrm{AX}}$ Hamiltonian can be expressed in terms of the total spin operators using algebraic tricks. We find an expression for the first term of Eq. (48) in terms of ${\widehat{\mathit{F}}}^{\mathrm{2}}$, ${\widehat{\mathit{F}}}_{\mathrm{A}}^{\mathrm{2}}$, and ${\widehat{\mathit{S}}}^{\mathrm{2}}$ by developing ${\widehat{F}}^{\mathrm{2}}$:

$$\begin{array}{}\text{(49)}& \begin{array}{rl}{\widehat{\mathit{F}}}^{\mathrm{2}}& ={\left(\widehat{\mathit{S}}+\sum _{l=\mathrm{1}}^n{\widehat{\mathit{I}}}_l\right)}^{\mathrm{2}}={\widehat{\mathit{S}}}^{\mathrm{2}}+\mathrm{2}\widehat{\mathit{S}}\cdot \sum _{l=\mathrm{1}}^n{\widehat{\mathit{I}}}_l+{\left(\sum _{l=\mathrm{1}}^n{\widehat{\mathit{I}}}_l\right)}^{\mathrm{2}}\\ & ={\widehat{\mathit{S}}}^{\mathrm{2}}+\mathrm{2}\sum _{l=\mathrm{1}}^n\widehat{\mathit{S}}\cdot {\widehat{\mathit{I}}}_l+{\widehat{\mathit{F}}}_{\mathrm{A}}^{\mathrm{2}}\\ & \iff \sum _{l=\mathrm{1}}^n\widehat{\mathit{S}}\cdot {\widehat{\mathit{I}}}_l={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left({\widehat{\mathit{F}}}^{\mathrm{2}}-{\widehat{\mathit{S}}}^{\mathrm{2}}-{\widehat{\mathit{F}}}_{\mathrm{A}}^{\mathrm{2}}\right).\end{array}\end{array}$$

Similarly, we find an expression for the second term of Eq. (48) in terms of ${\widehat{\mathit{F}}}_{\mathrm{A}}^{\mathrm{2}}$ and ${\widehat{\mathit{I}}}_l^{\mathrm{2}}$ by developing ${\widehat{\mathit{F}}}_{\mathrm{A}}^{\mathrm{2}}$:

$$\begin{array}{}\text{(50)}& \begin{array}{rl}{\widehat{\mathit{F}}}_{\mathrm{A}}^{\mathrm{2}}& ={\left(\sum _{l=\mathrm{1}}^n{\widehat{\mathit{I}}}_l\right)}^{\mathrm{2}}=\sum _{l=\mathrm{1}}^n{\widehat{\mathit{I}}}_l^{\mathrm{2}}+\mathrm{2}\sum _{l=\mathrm{1}}^{n-\mathrm{1}}\sum _{k>l}^n{\widehat{\mathit{I}}}_l\cdot {\widehat{\mathit{I}}}_k\\ & \iff \sum _{l=\mathrm{1}}^{n-\mathrm{1}}\sum _{k>l}^n{\widehat{\mathit{I}}}_l\cdot {\widehat{\mathit{I}}}_k={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left({\widehat{\mathit{F}}}_{\mathrm{A}}^{\mathrm{2}}-\sum _{l=\mathrm{1}}^n{\widehat{\mathit{I}}}_l^{\mathrm{2}}\right).\end{array}\end{array}$$

By substituting the results of Eqs. (49) and (50) into Eq. (48), we obtain a form of the Hamiltonian for which the energies will be more easily calculated:

$$\begin{array}{}\text{(51)}& \begin{array}{rl}{\widehat{H}}_{\mathrm{AX}}& =\mathrm{2}\mathit{\pi }{J}_{\mathrm{AX}}{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left({\widehat{\mathit{F}}}^{\mathrm{2}}-{\widehat{\mathit{S}}}^{\mathrm{2}}-{\widehat{\mathit{F}}}_{\mathrm{A}}^{\mathrm{2}}\right)\\ & +\mathrm{2}\mathit{\pi }{J}_{\mathrm{AA}}{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left({\widehat{\mathit{F}}}_{\mathrm{A}}^{\mathrm{2}}-\sum _{l=\mathrm{1}}^n{\widehat{\mathit{I}}}_l^{\mathrm{2}}\right).\end{array}\end{array}$$

The ${\widehat{H}}_{\mathrm{AX}}$ Hamiltonian commutes with the ${\widehat{\mathit{F}}}^{\mathrm{2}}$ operator; therefore, they share eigenstates |F,m_F〉. So, the eigenenergies can be written as the expectation values of |F,m_F〉 with respect to ${\widehat{H}}_{\mathrm{AX}}$:

$$\begin{array}{}\text{(52)}& E_{F,mF}=〈\left.F,m_F\right|{\widehat{H}}_{\mathrm{AX}}\left|F,m_F\right.〉.\end{array}$$

To calculate explicitly the eigenvalues, we substitute the Hamiltonian of Eq. (51) into Eq. (52) and use the following properties:

$$\begin{array}{}\text{(53)}& \begin{array}{rl}& {\widehat{\mathit{F}}}^{\mathrm{2}}\left|\right.F,m_F〉=F(F+\mathrm{1})\left|\right.F,m_F〉,\\ & {\widehat{\mathit{S}}}^{\mathrm{2}}\left|\right.F,m_F〉=S(S+\mathrm{1})\left|\right.F,m_F〉,\\ & {\widehat{\mathit{F}}}_{\mathrm{A}}^{\mathrm{2}}\left|\right.F,m_F〉=F_{\mathrm{A}}(F_{\mathrm{A}}+\mathrm{1})\left|\right.F,m_F〉,\\ & {\widehat{\mathit{I}}}_l^{\mathrm{2}}\left|\right.F,m_F〉=I_l(I_l+\mathrm{1})\left|\right.F,m_F〉,\end{array}\end{array}$$

to obtain the final expression for the energy of level |F,m_F〉.

$$\begin{array}{}\text{(54)}& \begin{array}{rl}{E}_{F,mF}& ={\displaystyle \frac{J_{\mathrm{AX}}}{\mathrm{2}}}\left[F(F+\mathrm{1})-S(S+\mathrm{1})-F_{\mathrm{A}}(F_{\mathrm{A}}+\mathrm{1})\right]\\ & +{\displaystyle \frac{J_{\mathrm{AA}}}{\mathrm{2}}}\left[F_{\mathrm{A}}(F_{\mathrm{A}}+\mathrm{1})-n{I}_l(I_l+\mathrm{1})\right],\end{array}\end{array}$$

expressed in hertz. Here, quantum number F corresponds to the total spin of the full XA_n system, S corresponds to the spin of the nucleus X, F_A is the total spin of the A_n spins, and I_l is the spin of individual nuclei A. The energy does not depend on the spin projections, resulting in degeneracy of all 2F+1 levels with equal F.

The spin number S is the same for all eigenstates (e.g., it is $\mathrm{1}/\mathrm{2}$ for ¹³C); similarly, all spin numbers I_l are the same (e.g., for ¹H spins they are equal to $\mathrm{1}/\mathrm{2}$). The remaining two quantum numbers F and F_A can have different values depending on the state, therefore removing degeneracy between some of the levels. Figure 8 presents the energy levels of XA, XA₂, and XA₃ systems at ZF calculated using Eq. (54). Mathematica codes to perform these calculations are available in the Supplement (Supplement S6).

Figure 8 Energy levels for XA, XA₂, and XA₃ spin systems calculated according to Eq. (54). The numbers above the energy levels represent the z projection of the angular momentum of the states m_F. Allowed transitions are shown by green arrows. J_AX was set to 140 Hz, and J_AA was set to −12 Hz; these are typical J coupling values for ¹J_CH and ²J_HH. The energy difference for the allowed transitions equals to J_AX for the XA system, $\mathrm{3}/\mathrm{2}{J}_{\mathrm{AX}}$ for the XA₂ system, and two frequencies of 2J_AX and of J_AX for the XA₃ system. This agrees with the numerical simulations shown in Fig. 6.

4.3 Selection rules

We have now found the eigenstates and their energies, but not all transitions between any pair of states are allowed. The last step is to find the transition intensities and thus get the analytical appearance for the ZF NMR spectrum of an XA_n system. There are certain selection rules specifying which transitions are in principle possible and which are forbidden, like those in high-field NMR, where only single quantum transitions are allowed. A general expression for the transition intensity between any two eigenstates |F,m_F〉 and $\left|F^{\prime },{m^{\prime }}_F〉\right.$ is given by

$$\begin{array}{}\text{(55)}& Y=〈F^{\prime },{m^{\prime }}_F\mathrm{|}{\widehat{\mathit{\rho }}}_{\mathrm{0}}\mathrm{|}F,m_F〉〈F^{\prime },{m^{\prime }}_F\mathrm{|}\widehat{O}\mathrm{|}F,m_F〉.\end{array}$$

We will explicitly calculate the transition intensity for the sudden field drop experiment. In this case, both the initial density operator ${\widehat{\mathit{\rho }}}_{\mathrm{0}}$ and the observation operator $\widehat{O}$ (see Eq. 35) are proportional to ${\mathit{\gamma }}_I{\widehat{F}}_{\mathrm{A},z}+{\mathit{\gamma }}_S{\widehat{S}}_z$ (as a reminder, ${\widehat{F}}_{\mathrm{A},z}=\sum _{l=\mathrm{1}}^n{\widehat{I}}_{l,z}$). Therefore, the transition intensity becomes

$$\begin{array}{}\text{(56)}& Y={〈F^{\prime },{m^{\prime }}_F\mathrm{|}{\mathit{\gamma }}_I{\widehat{F}}_{\mathrm{A},z}+{\mathit{\gamma }}_S{\widehat{S}}_z\mathrm{|}F,m_F〉}^{\mathrm{2}}.\end{array}$$

This expression is an example of Fermi's golden rule, which is used to calculate transition amplitudes in different problems in quantum mechanics. Similar expressions can be found for the high-field NMR. By expressing the coupled states |F,m_F〉 in terms of the uncoupled basis (see Eq. 46), we find that

$$\begin{array}{}\text{(57)}& \begin{array}{rl}& \left({\mathit{\gamma }}_I{\widehat{F}}_{\mathrm{A},z}+{\mathit{\gamma }}_S{\widehat{S}}_z\right)\left|F,m_F\right.〉\\ & =\left({\mathit{\gamma }}_I{\widehat{F}}_{\mathrm{A},z}+{\mathit{\gamma }}_S{\widehat{S}}_z\right)\sum _{m_{\mathrm{A}},m_S}{C}_{F_{\mathrm{A}},m_{\mathrm{A}},S,m_S}^{F,m_F}\left|F_{\mathrm{A}},m_{\mathrm{A}},S,m_S\right.〉\\ & =\sum _{m_{\mathrm{A}},m_S}{C}_{F_{\mathrm{A}},m_{\mathrm{A}},S,m_S}^{F,m_F}\left({\mathit{\gamma }}_I{m}_{\mathrm{A}}+{\mathit{\gamma }}_S{m}_S\right)\left|F_{\mathrm{A}},m_{\mathrm{A}},S,m_S\right.〉,\end{array}\end{array}$$

where m_A and m_S are the z projection of the total spins F_A (for n protons, the maximum value of F_A equals $n/\mathrm{2}$, and for each value of F_A, $m_{\mathrm{A}}\in \mathit{\{}-F_{\mathrm{A}},-F_{\mathrm{A}}+\mathrm{1},\mathrm{\dots },F_{\mathrm{A}}-\mathrm{1},F_{\mathrm{A}}\mathit{\}}$) and the z projection of the spin S (in the case where S is a carbon-13 spin, $m_S\in \mathit{\{}-\frac{\mathrm{1}}{\mathrm{2}},\frac{\mathrm{1}}{\mathrm{2}}\mathit{\}}$), respectively. Now let us express the remaining $〈\left.F^{\prime },{m^{\prime }}_F\right|$ bra in terms of the uncoupled basis as well and combine Eqs. (56) and (57).

$$\begin{array}{}\text{(58)}& \begin{array}{rl}Y& =(\sum _{m_{\mathrm{A}}^{\prime },m_S^{\prime }}{C}_{F_{\mathrm{A}}^{\prime },m_{\mathrm{A}}^{\prime },S^{\prime },m_S^{\prime }}^{F^{\prime },m_F^{\prime }}〈\left.F_{\mathrm{A}}^{\prime },m_{\mathrm{A}}^{\prime },S,m_S^{\prime }\right|\\ & \sum _{m_{\mathrm{A}},m_S}{C}_{F_{\mathrm{A}},m_{\mathrm{A}},S,m_S}^{F,m_F}\left({\mathit{\gamma }}_I{m}_{\mathrm{A}}+{\mathit{\gamma }}_S{m}_S\right)\left|F_{\mathrm{A}},m_{\mathrm{A}},S,m_S\right.〉{)}^{\mathrm{2}}\\ & =(\sum _{{m^{\prime }}_{\mathrm{A}},{m^{\prime }}_S,m_{\mathrm{A}},m_S}{C}_{F_{\mathrm{A}}^{\prime },m_{\mathrm{A}}^{\prime },S^{\prime },{m^{\prime }}_S}^{F^{\prime },{m^{\prime }}_F}{C}_{F_{\mathrm{A}},m_{\mathrm{A}},S,m_S}^{F,m_F}\\ & \left({\mathit{\gamma }}_I{m}_I+{\mathit{\gamma }}_S{m}_S\right)〈F_{\mathrm{A}}^{\prime },m_{\mathrm{A}}^{\prime },S^{\prime },{m^{\prime }}_S\left|\right.F_{\mathrm{A}},m_{\mathrm{A}},S,m_S〉{)}^{\mathrm{2}}\end{array}\end{array}$$

The last term $〈F_{\mathrm{A}}^{\prime },m_{\mathrm{A}}^{\prime },S^{\prime },{m^{\prime }}_S\left|\right.F_{\mathrm{A}},m_{\mathrm{A}},S,m_S〉$ is nonzero only if

$$\begin{array}{}\text{(59)}& \begin{array}{rl}& \mathrm{\Delta }{F}_{\mathrm{A}}=F_{\mathrm{A}}^{\prime }-F_{\mathrm{A}}=\mathrm{0},\\ & \mathrm{\Delta }{m}_{\mathrm{A}}=m_{\mathrm{A}}^{\prime }-m_{\mathrm{A}}=\mathrm{0},\\ & \mathrm{\Delta }S=S^{\prime }-S=\mathrm{0},\\ & \mathrm{\Delta }{m}_S=m_S^{\prime }-m_S=\mathrm{0}.\end{array}\end{array}$$

These selection rules mean that the only allowed transitions are those which conserve the total spins F_A and S (S is conserved automatically as it can be only $\mathrm{1}/\mathrm{2}$, but F_A can have different values), as well as their projections onto the reference axis. Equation (58) therefore simplifies to

$$\begin{array}{}\text{(60)}& Y={\left(\sum _{m_{\mathrm{A}},m_S}{C}_{F_{\mathrm{A}},m_{\mathrm{A}},S,m_S}^{F^{\prime },{m^{\prime }}_F}{C}_{F_{\mathrm{A}},m_{\mathrm{A}},S,m_S}^{F,m_F}\left({\mathit{\gamma }}_I{m}_{\mathrm{A}}+{\mathit{\gamma }}_S{m}_S\right)\right)}^{\mathrm{2}}.\end{array}$$

It is important to notice that, in the case where γ_I=γ_S, this sum is always null. This is shown in the Wolfram Mathematica code for all observable transitions in XA, XA₂, and XA₃ systems and can be rationalized in the general case by the following (see Supplement S6). The operator ${\mathit{\gamma }}_I{\widehat{F}}_{\mathrm{A},z}+{\mathit{\gamma }}_S{\widehat{S}}_z$ (which is proportional to the initial density operator ${\widehat{\mathit{\rho }}}_{\mathrm{0}}$) can be rewritten as ${\mathit{\gamma }}_I\left({\widehat{F}}_{\mathrm{A},z}+{\widehat{S}}_z\right)+\left({\mathit{\gamma }}_S-{\mathit{\gamma }}_I\right){\widehat{S}}_z$. The first term in this expression commutes with the ${\widehat{H}}_{\mathrm{AX}}$ Hamiltonian (see Eq. 48); therefore, it does not produce any observable coherences, whereas the second term does not commute with the ${\widehat{H}}_{\mathrm{AX}}$ and leads to ZULF signals.

Finally, there are two more selection rules that are derived by implementing the Wigner–Eckart theorem. The considered case is equivalent to a “dipole” transition, where the transition is observed between two states connected by operator of rank 1 (e.g., Eq. 56). This is a common situation in atomic physics, and we adapt this result without evaluation: the reduced matrix element coming from Wigner–Eckart theorem is shown to be nonzero if and only if

$$\begin{array}{}\text{(61)}& \begin{array}{c}\mathrm{\Delta }F=\pm \mathrm{1},\\ \mathrm{\Delta }{m}_F=\mathrm{0}.\end{array}\end{array}$$

The whole set of selection rules given by Eqs. (59) and (61) allows one to find which transitions are observable in XA_n systems at ZF. These transitions are shown in Fig. 8 by the green arrows. It can be seen that J_AA couplings shift the energy levels but do not affect the frequencies of the observable transitions. This is a common situation that J couplings between magnetically equivalent spins do not contribute to the observed NMR spectrum. As can be seen from the analysis presented above, this statement holds for each case of the ZF NMR spectra of XA_n systems.

In this section, we analytically found the allowed transitions for XA, XA₂, and XA₃ for the case of the sudden field drop to ZF. The XA spin system has a single transition at J_AX, the XA₂ spin system has a single transition at $\mathrm{3}/\mathrm{2}{J}_{\mathrm{AX}}$, and the XA₃ spin system has one allowed transition at J_AX and another one at 2 J_AX. The allowed transitions analytically found here correspond to the numerical simulation: XA single line at J_AX, XA₂ single line at $\mathrm{3}/\mathrm{2}{J}_{\mathrm{AX}}$, etc. This derivation explained the appearance of the ZF spectra but not that of the ultralow-field spectra. To understand how the degeneracy of the ZF eigenstates are split by the presence of a bias field, one has to use perturbation theory. We refer the interested reader to Ledbetter et al. (2011) and Appelt et al. (2010).

4.4 Rabi oscillation curves

We finish this section on the interpretation of the numerical simulations by giving a short explanation of the Rabi oscillation curves presented in Sect. 3.2 (see Fig. 5). The successful implementation of excitation pulses in ZULF-NMR experiments requires two conditions to be fulfilled (Butler et al., 2013b). First, the constant magnetic field of the pulse should be strong enough so that heteronuclei (here, ¹H and ¹³C spins) can be considered weakly coupled during the pulse. The field of 50 µT satisfies this condition, as the difference in Larmor frequencies between ¹H and ¹³C spins is larger than 1.5 kHz ≫ J_XA=140 Hz. Second, the pulse must be much shorter than the evolution under the J coupling. Here, the longest pulses that were simulated had a duration of 3 ms, while the characteristic time of the evolution under the J coupling is $\mathrm{1}/J_{\mathrm{XA}}\approx \mathrm{7.1}$ ms. Provided these two conditions are met, the product operator formalism can provide a convenient explanation for the results of Fig. 5. Both Rabi oscillation curves in Fig. 5 are rather unusual for high-field NMR, but the reader who is familiar with the product operator formalism at high field will see that there is a strong connection between the algebra describing pulsed experiments at high field and at ZULF. Here, we give a brief summary of how this formalism can be used to understand the Rabi oscillation curves. We recommend the interested reader to look at the following references for a more detailed derivation (Butler et al., 2013b; Blanchard, 2014; Tayler et al., 2017).

After the adiabatic field drop, the magnetization of the sample is proportional to ${\mathit{\gamma }}_{\mathrm{H}}{\widehat{I}}_z+{\mathit{\gamma }}_{\mathrm{C}}{\widehat{S}}_z$ and does not evolve. In addition, part of the population is also on the zero-quantum term ${\widehat{Z}}_z=\mathrm{2}\left({\widehat{I}}_x{\widehat{S}}_x+{\widehat{I}}_y{\widehat{S}}_y\right)$, which produces no observable magnetization. Magnetic field pulses are applied to convert one or both of these terms into the observable zero-quantum term ${\widehat{Z}}_x={\widehat{I}}_z-{\widehat{S}}_z$. In the case of Fig. 5a where the pulse is applied along the z axis, the pulse converts ${\widehat{Z}}_z$ into ${\widehat{Z}}_y=\mathrm{2}\left({\widehat{I}}_x{\widehat{S}}_y+{\widehat{I}}_y{\widehat{S}}_x\right)$. The state of the system after the pulse is $\widehat{\mathit{\rho }}\left({\mathit{\tau }}_{\mathrm{p}}\right)={\widehat{Z}}_z\cos \left[\left({\mathit{\gamma }}_{\mathrm{H}}-{\mathit{\gamma }}_{\mathrm{C}}\right)B_z{\mathit{\tau }}_{\mathrm{p}}\right]+{\widehat{Z}}_y\sin \left[\left({\mathit{\gamma }}_{\mathrm{H}}-{\mathit{\gamma }}_{\mathrm{C}}\right)B_z{\mathit{\tau }}_{\mathrm{p}}\right]$. The excited unobservable ${\widehat{Z}}_y$ term then starts to evolve into the observable term ${\widehat{Z}}_x$ under the action of the ${\widehat{H}}_J$ Hamiltonian, generating an oscillating magnetic field along the z axis. The resulting ZULF signal has a sine rather than a cosine time dependence and requires a $\mathit{\pi }/\mathrm{2}$ phase correction of the spectrum to have an absorption line as was described in Sect. 3.2. The signal is maximized when the pulse has duration ${\mathit{\tau }}_{\mathrm{p}}=\mathit{\pi }/\left[\mathrm{2}\left({\mathit{\gamma }}_{\mathrm{H}}-{\mathit{\gamma }}_{\mathrm{C}}\right)B_z\right]$, which is around 157 µs in the considered case. This is consistent with the simulated Rabi oscillation curve of Fig. 5a. In the case of Fig. 5b where the pulse is applied along the x axis, both initial terms of the density operator, ${\mathit{\gamma }}_{\mathrm{H}}{\widehat{I}}_z+{\mathit{\gamma }}_{\mathrm{C}}{\widehat{S}}_z$ and ${\widehat{Z}}_z$, are converted into the observable term ${\widehat{Z}}_x$. The conversion follows a $\sin \left[\left({\mathit{\gamma }}_{\mathrm{H}}+{\mathit{\gamma }}_{\mathrm{C}}\right)B_x{\mathit{\tau }}_{\mathrm{p}}/\mathrm{2}\right]\sin \left[\left({\mathit{\gamma }}_{\mathrm{H}}-{\mathit{\gamma }}_{\mathrm{C}}\right)B_x{\mathit{\tau }}_{\mathrm{p}}/\mathrm{2}\right]$ function, allowing one to excite slightly stronger signals over a slightly longer pulse duration.

5 Conclusion

We have shown how to numerically simulate spectra at both zero and ultralow fields for sudden drop and pulsed experiments. We have then explained the results of the numerical simulation for sudden field drop experiments at ZF by constructing the eigenbasis of the ZF Hamiltonian and finding the allowed transitions among the eigenstates. The other numerically simulated cases (i.e., pulsed experiments) can be explained using the analytical approach that we have presented here. It requires an additional step which is to describe how a pulse converts the populations of the states. The reader who is acquainted with the product operator formalism commonly used in high-field NMR might be interested in an alternative approach based of commutation rules as presented in Blanchard and Budker (2016) and Butler et al. (2013b). We have chosen to describe the simplest cases, i.e., experiments with thermal prepolarization with AX_n systems. Using this methodology, the reader can proceed with simulating more advanced cases, where analytical solutions do not exist. This includes calculation of ZULF spectra of molecules with multiple spins (Wilzewski et al., 2017) and molecules containing three or more types of nuclei, e.g., ¹H, ¹³C, ¹⁵N, and ²D (Alcicek et al., 2021); the evolution during dynamical decoupling sequences (Bodenstedt et al., 2022a); the ZULF-TOCSY type of spin-locking experiments (Kiryutin et al., 2021); spin evolution at intermediate fields, where perturbation approaches are not valid (Bodenstedt et al., 2021); or the complicated spin dynamics that may occur under the action of composite pulses (Jiang et al., 2018; Bodenstedt et al., 2022b). The formalism we presented here is a good starting point for the description and understanding of hyperpolarized ZULF experiments, e.g., those involving transfer of spin order in parahydrogen experiments at low fields. Simulations are also useful to study different kinds of imperfections such as field inhomogeneities, timing errors, etc. We hope that this tutorial paper has allowed us to share our excitement with the reader.