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

The supplement related to this article is available onlineat:https://doi.org/10.5194/mr-4-87-2023-supplement.

See, for example,http://spindynamics.org/wiki/index.php?title=Zerofield.m(last access: 21March2023).

2.1. Define the experimental sequence

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

During the preparation, some nuclear polarization is acquired byletting the sample rest in a strong magnetic field (in most conventionalexperiments). Mixing consists of bringing the system to a non-stationary statewhose oscillations are recorded during detection. In common high-field NMRexperiments, all the steps are performed in a strong magnet with a nearlyconstant magnetic field. Nuclear polarization is spontaneously acquired due tothe high magnetic field, and both the mixing and detection are performed throughthe same RF coil using Faraday induction. At ZULF, there is no nuclearpolarization, so the preparation has to be performed in different conditions. Acommon method is to shuttle the sample between a region of highfield and aregion of ZULF.

Figure 1 shows a typical experimental setup. A permanent magnet isused to prepolarize the sample, which is connected with the magnetic shields bya guiding solenoid coil. This coil ensures that the sample experiences amagnetic field with constant direction and sufficient strength during thetransfer from the region of highfield to inside the magnetic shields (i.e., thecoil ensures an adiabatic transfer). Once the sample arrives in the magneticshields at the location of detection, the Helmholtz coil continues to produce amagnetic field in the same direction as the solenoid, and the spin system isstill distributed into Zeeman populations (Blanchard and Budker, 2016; Tayler etal., 2017). All the steps detailed until here are part of the preparation. Inpractice, the guiding solenoid and the Helmholtz coil produce a magnetic fieldwhich is much weaker than the prepolarizing magnet. However, this will not betaken into account in the simulation: we consider that the sample spends enoughtime in the prepolarizing magnet to reach Boltzmann equilibrium and that thetransfer is sufficiently fast for us to neglect the change in polarizationduring the transfer.

A further step can optionally be added to the preparation, whichconsists of ramping down the magnetic field produced by the Helmholtz coil tobring the spins adiabatically to ZF. We will refer to experiments which includeor do not include this step as the adiabatic field drop and sudden field dropexperiments, respectively. In the case of sudden field drop experiments, themixing step simply consists of switching off the magnetic fieldnon-adiabatically (that is, fast enough to be considered instantaneous withrespect to the evolution of the spin system). In the case of adiabatic fielddrop experiments, the sample is already at ZF at the end of preparation, sopopulations have to be mixed by applying a magnetic field pulse before anysignal can be detected. This is analogous to high-field pulses except that ituses constant magnetic field rather than RF pulses. After the mixing, theoscillating magnetic field generated by the sample is detected by an opticalmagnetometer. In Fig.1, the magnetometer is represented below the sample, thatis, aligned along the zaxis with respect to the sample. We assume that the OPM isconfigured so as to be sensitive to magnetic field along the zaxis and that the spins are initially prepolarized along thesame axis. Defining this axis as z is a natural choice for high-field NMR spectroscopists, butnote that other conventions exist (see, for example, Ledbetter et al., 2011).During detection, a weak magnetic field may be applied, either along the zaxis or along an orthogonal axis. In the latter case, theexperiment is said to be performed under the ULF regime. In the absence ofapplied magnetic field (and provided the residual magnetic field is properlyzeroed at the location of the sample), the experiment is said to be performedunder the ZF regime.

In summary, there are several possible combinations of experimentalschemes. All of them start with prepolarizing the spins at high magnetic field.After the sample is transported into the magnetic shields, the field is droppedeither suddenly or adiabatically, in which case a magnetic field pulse isapplied. Finally, the oscillating magnetic field produced by the sample isdetected along the zaxis, with or without a weak magnetic field applied along the xaxis. In the remaining of the paper, these sequences presentedin 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 themolecule whose ZULF spectrum is to be simulated and the interactions which thespins are subject to. This paper is concerned with small molecules in the liquidstate. As is the case for high-field NMR, dipolar interactions are averaged outby rapid molecular tumbling and need not be taken into account (except asstochastic perturbation if one intends to include relaxation effects).Therefore, only the Jcoupling and the Zeeman interactions are considered here.

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

2.3. Compute the spin Hamiltonian

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

Let us first assume that the system contains a single spin1/2 interacting with the magnetic field B along the zaxis. The state of any spin system can be represented as alinear combination of basis vectors, which are called “kets” in Dirac's notationand are represented by the symbols |〉. For a single spin1/2, two basis kets are necessary to represent the state of thesystem. We chose to represent the spin system in the Zeeman basis:

The α and β states correspond to the spin being parallel and antiparallelwith the magnetic field, respectively. The general state in which the spin maybe found is a linear combination of these two basis states. Because these statesand their associated kets form an orthonormal basis, their vector representationhave the canonical form with only 0 and 1 coefficients. The space spanned bythese two vectors is called a “Hilbert space” and has dimension 2, as indicatedby the superscript in BZ2. Note that the choice of the Zeeman basis is convenient fornumerical simulation but is not necessary. For example, one may use the coupledbasis, which will be presented and used in Sect.4. The same basis may be usedto perform simulations at highfield or ZULF, although one particular basismight turn out to be more convenient.

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

These operators act on the Zeeman states in certain ways, e.g., I^xα=12β. To summarize the set of rules, it is convenient to use thematrix representation of the operators, with the matrix elements determined bythe action of the operator on the α and β states: Iμrs=rI^ks, where r,s∈α,β;μ∈{x,y,z}. This definition makes use of 〈|, i.e., the “bra”, an object which is complementary to the ketand corresponds to the “Hermitean conjugate” of the ket. In the matrixrepresentation of quantum mechanics, the Hermitean conjugate of a ketcorresponds to the complex transpose of the vector representing the ket. Thematrix representations of operators in quantum mechanics is very important forperforming simulations, as they are constructed in such a way that any state ofoperation on the quantum system can be represented using linear algebra. Thematrix representations of the three Cartesian components of the spin angularmomentum operators are proportional to Pauli matrices σ^x, σ^y, and σ^z:

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

where γ is the gyromagnetic ratio of the spin (inrad s ¹T^-1). The dot product of the vectors ofthe magnetic field and of the spin angular momentum (vectors and vectoroperators are denoted in bold throughout the text) is expanded on the rightmember of Eq.(4). Note that we have omitted the reduced Planck constant ℏ in Eq.(4), which implies that the energy is expressed inradians per second (rad s^-1) rather than in joules. This is thecase throughout this paper. In many cases, the magnetic field is aligned withone of the axes. If it points along the zaxis, i.e., B=00B0, Eq.(4) simplifies to

where ω0=-γB0 is the Larmor frequency of the spin. This expression is validregardless of the intensity of the magnetic field, i.e., at high field as wellas at ZULF. The Zeeman states, α and β, which correspond to the spin being parallel and antiparallelwith the magnetic field, respectively, are eigenstates of the ZeemanHamiltonian; that is, they satisfy the eigenequations H^Zα=+1/2α and H^Zβ=-1/2β. The eigenstates of a Hamiltonian are of particularimportance; they are states which do not evolve under the effect of thatHamiltonian (ignoring the accumulation of the phase factor, which turns out tobe irrelevant in most experiments), i.e., stationary states.

The single spin whose Hamiltonian is given by Eq.(5) lives in aHilbert space of dimension 2. To represent a pair of spins1/2, we need to use a Hilbert space with a dimension of 4. To doso, we redefine the angular momentum operators in this higher-dimension space.The matrix representations of the angular momentum operators I^1x, I^1y, and I^1z and I^2x, I^2y, and I^2z of spin 1 and spin 2, respectively, are given by the Kroneckerproduct of matrices of single-spin angular momentum operator and the identityoperator, in the appropriate order. For the z-axis angular momentum operators, we have

and

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

The two operators defined by Eqs.(6) and (7) are the same asthe one given by Eq.(5), except that the world of spin1 now contains spin2,and vice versa. This representation corresponds to a basis that is the Kroneckerproduct of the basis of the individual spins.

For the case where the magnetic field points along the zaxis, the total Zeeman Hamiltonian for the two spins can nowbe computed using Eq.(5) in the basis of Eq.(9) as the sum of the two ZeemanHamiltonians:

where ω10 and ω20 are the Larmor frequencies of spin 1 and 2, respectively. Notethat in a Hilbert space of several spins, it is useful to define projections oftotal angular momentum operators:

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

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

where J is the Jcoupling between the two spins (in Hz). Compared with theZeeman Hamiltonian (see Eq.10), the J-coupling Hamiltonian has the particularity to haveoff-diagonal elements in the αβ,βα subspace, which implies that the Jinteraction mixes the αβ and βα states. In other words, due to the Jinteraction, these two states are no longer eigenstates of thespin system.

In the case of a system of n spins 1/2, the same procedure can be applied to define the angularmomentum operators and the Hamiltonians. These operators can be represented as 2n×2n 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

where and I^kz is the zangular momentum operator of spin k in an n-spin Hilbert space, and 1^2×2 and I^z2×2 are the identity operator and the z angular momentum operator in a single-spin Hilbert space,respectively. The zprojection of total angular momentum operators is given by

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

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

where Jlk is the Jcoupling between spins l and k (in Hz). Because a spin is not Jcoupled to itself, the sum in Eq. (16) does not include termswith 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 JHamiltonian in Eqs.(15) and (16), respectively, are validboth at high field and at ZULF. What makes the difference between the tworegimes 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 bya density operator. If ψ is a ket representing the state of the system as a linearcombination of basis states (like those defined in Eqs.1 and 9), the densityoperator is given by

where the upper bar represents the ensemble average over allidentical spin systems in the sample – the operation performed by the densityoperator. This averaging makes the density operator formalism well-suited forNMR, where the experiment consists of observing a large number of identical spinsystems at the same time rather than a single spin system. The matrixrepresentation of the density operator (and of any other spin operator) isachieved by calculating all the matrix elements ρrs=rρ^s, where r and s are basis states. For example, the matrix representation ofthe density operator for the α and β states of a single spin yields

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

where H^, kB, and T are the Hamiltonian operator of the spin system, Boltzmann'sconstant, and the temperature, respectively. Operation exp⁡() denotes the matrix exponentiation. Note that this operationdoes not consist of applying fx=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 densitymatrix has unit trace. It is given by

The prepolarizing step of the experiments that we intend tosimulate occurs in a strong magnetic field (in the sense that the Zeemaninteraction is largely dominating all other interactions), as in a standardhigh-field experiment. In this case, we can compute the thermal equilibrium bytaking only the Zeeman terms into account. For a single spin with Larmorfrequency ω0 and gyromagnetic ratio γ in prepolarizing field Bp, the thermal equilibrium density matrix yields

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

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

For an n-spin system, we take the Kronecker product of density matricesof individual spins ρ^eq,l2×2.

The first approximation in this expression comes from the factthat it neglects all spin–spin interactions. The second approximation comes fromthe fact that we neglect multiple-spin order terms, which are proportional tothe products PlPk, where Pl and Pk are the polarizations of spin l and k, respectively. This approximation is valid unless the systemis highly polarized, which is the case even at very high field (withouthyperpolarization). To avoid confusion, we specified that the operators ρ^eq,l2×2, 1^2×2, and I^z2×2 act on a single-spin Hilbert space (2×2 matrix). On the contrary, the operators 1^ and I^lz act on spin states of nspins, and accordingly their matrix representations havedimension of 2n×2n (for spins1/2). As shown by Eq.(23), one may compute the density matrixeither using the Kronecker product of operators in a single-spin Hilbert spaceor by summing the operators in a Hilbert space of nspins.

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

where n is the number of spins in the system. Because all operatorscommute with the identity, this does not affect the result of propagation. Theresulting expression is simpler (ρ^eq=PI^z2×2 for a single spin), which is convenient for calculations byhand. It may also make the numerical propagation faster and more precise.Another common simplification is to drop the polarization factor. For a singlespin, the two combined simplifications yield

Simplifications are useful, but they should be handled withcare. The polarization factor P is different for spins with different gyromagnetic ratio. Ifit is dropped without introducing further corrections, the relative sizes of thepopulation of spins with different gyromagnetic ratios will not be respected. Inthe simulations presented here, we will compute the initial density matrix usingthe 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 thematrix representation of the Hamiltonian. We now describe how the evolution ofthe system (represented by the density matrix) evolves with time under a givenHamiltonian. This will be used at several steps of the simulation: when thesample is brought adiabatically to ZF, during the pulse, and during the signalmeasurement.

The evolution of a quantum system with time is given by thetime-dependent Schrödinger equation. Its equivalent for the evolution of densitymatrix is the Liouville–vonNeumann equation ddtρ^t=-iH^t,ρ^t, which has the solution

where ρ^0 is the density matrix at t=0, and U^ is the propagator during time t, which is defined as

where H^ is the total Hamiltonian. The operation of Eq.(26) “takes”the spin system from ρ^0 to ρ^t. Again, note that exp⁡() denotes the matrix exponentiation and not theelement-by-element exponentiation. An important case of propagator is therotation operator. For an angular momentum operator I^μ, with μ∈{x,y,z}, the propagator exp⁡-iI^μθ is called a rotation operator; it represents a rotation of thespins of angle θ around axis μ, when applied to the density matrix using Eq.(26). For asingle spin subject to a static magnetic field along the zaxis, the total Hamiltonian is the Zeeman Hamiltonian (seeEq.5) which causes the spin to rotate around the zaxis; this rotation can be expressed using the rotationoperator exp⁡-iH^t=exp⁡-iω0I^zt with angle ω0t.

The brute force calculation of the exponential operator in anarbitrary basis is computationally challenging as it requires calculating theTaylor expansion of the H^ operator, which can be calculated analytically only in a fewcases, like the propagators of the Iz, Ix, and Iy operators. To avoid this, the calculation of the propagator(Eq.27) is usually performed by diagonalizing the Hamiltonian and then takingthe complex exponent for each of its eigenvalues, exp⁡-iωkt, where ωk denotes the kth eigenvalue. Therefore, the transformation to the eigenbasisof 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 thelinear algebra packages of the software. One may note that the basis does notaffect the result of the calculation, but the choice of a more appropriate onemay help rationalize the problem. In many cases, the initial choice is theZeeman basis, in which spin operators are readily introduced based on Kroneckerproducts of the Pauli matrices. Depending on the symmetry of the problem, itmight be more convenient to change the basis to another one. As we will see inSect.4.1, a choice of coupled basis is preferable for understanding thezero-field Jspectroscopy of coupled spins.

It is important to remark that Eq.(27) is only valid if theHamiltonian is constant during the evolution period. The case where theHamiltonian is time dependent is treated below. Note that the propagator is aunitary operator and therefore has the convenient property that its inverse isequal to its complex transpose (i.e., U^-1=U^†), which is much faster to compute than the matrix inverse U^-1.

Equations(26) and (27) allow us to know the state of the system atany time t from the initial time t=0. To simulate the signal produced by the spin system during thecourse of the experiment, we must calculate the time domain signal at differenttime points. Note that in this case the Hamiltonian remains constant during freeevolution. To calculate the signal at fixed time steps, it is convenient tofirst calculate the propagator U^(dt) over period dt. We then apply Eq.(26) recursively to get the new densitymatrix ρ^tk+1 from the previous one ρ^tk,

where tk+1-tk=dt. To simulate ZULF spectra, we will also encounter situationswhere the Hamiltonian is time dependent. First, the Hamiltonian can vary withtime but be “constant by block”. This is, for example, the case for the suddenfield drop; the system is under a certain Zeeman Hamiltonian in the beginning ofthe experiment and suddenly under the ZULF Hamiltonian during detection. Thissituation does not present particular difficulties; the evolution of the systemcan be described step by step by both Eqs.(26) and (28).

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

where the propagator is given by

where H^(tk) is the Hamiltonian at time tk. Note that the choice of H^(tk) rather than H^(tk+1) 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 thedensity matrix along time. To simulate the time domain signal, we need toextract a physical quantity from the density matrix as it evolves with time. Themeasured physical quantity of a ZULF experiment is the magnetic field producedby the nuclear spins of the sample at the location of an OPM. In a firstapproximation, we can consider that the whole sample is a point dipoleinteracting with the OPM and that this total dipole is the sum of the dipoles ofthe individual spin systems (Fig.2 gives a visual representation of theapproximation). Whether this approximation is appropriate or not depends on thegeometry of the experimental setup. We have chosen the zaxis as the quantization axis (defined by the detector, i.e.,the OPM). Therefore, the physical quantity that we need to compute is the totalmagnetic field produced by the spins along the zaxis at the location of the vapor cell:

where 〈μ^ztot〉, μ0, N, 〈μ^z〉, and r are the magnetic moment of the sample along the zaxis, the permeability of free space, the number of identicalspin systems in the sample, their individual magnetic moments along the zaxis, and the distance between the center of the sample andthe center of the vapor cell, respectively.

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

where μ^lz, γl, and I^lz are the magnetic moment, the gyromagnetic ratio, and theangular momentum along the zaxis of spin l, respectively. Note that n and N represent the number of spins in the molecule and the numberof molecules in the sample, respectively. The notation 〈〉 denotes the expectation value of a quantity. Particularlyimportant ones are those that can be physically measured in the experiment. Inthe density matrix formalism that we are using, the expectation value of aphysical quantity related to an operator A^ is given by

where Tr{} denotes the matrix trace, i.e., the sum of all diagonalelements of the matrix representation of the operator. Note that the expectationvalue of μ^lz (or I^lz) is proportional to the polarization level of spin l which was accounted for in Eqs.(21) and (22). Therefore, thetotal magnetic moment calculated with Eq.(32) depends on the polarization ofthe different spin species.

If ρ^t is the density matrix at time t, we obtain the signal S(t) measured by the OPM by plugging Eq.(32) into Eq.(33)

where we have defined a “detection operator”,

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

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

2.7. Fourier transform the expectation values to obtain a spectrum

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

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

where Sk, tk, lB, and T2 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 toas the spin–spin relaxation constants or transverse relaxation time constant.The signal intensities Sk′ define the apodized signal array S′. As shown in Eq.(36), we may choose to express theapodization function using either the coherence time constant T2 or the line broadening lB (in Hz), which are related by πlB=1/T2. The former is the time constant at which the time domainsignal decays, while the latter is the full width at half height (FWHH) of thefrequency domain signals. In order to avoid “truncating” the decay of the timedomain signal and the related spectral artifacts, we must fulfill the condition T2≪taq, where taq=max⁡tk is the acquisition time (or the length of the signal in thetime domain). Typically, we may choose T2 and taq so that taq=5T2. Table1 summarizes the parameters which were used in thispaper.

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

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

In MATLAB programming language, the function for fast Fouriertransformation fft() takes array S′′ as input and returns the frequency domain array whichcorresponds to the simulated spectrum. Optionally, one may add a second argument L to fft() to include zero-filling in theFourier transform. Including zero-filling has the advantage of increasing thenumber of points per FWHH on the spectrum without increasing the computationtime of the propagation. Due to MATLAB's Fourier transform convention, it isconvenient to retransform the signal with fftshift() in orderto obtain a Fourier transformed signal with 0 as the middle frequency. We thendivide the output of MATLAB's Fourier transform by the number of points L:

where F designates the Fourier transform. The frequency domain signalobtained 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 butpreserves the integrals.

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

where the sampling frequency (in Hz) is given by

The sampling frequency of the time domain signal gives the maximumfrequency that can be appropriately sampled. Figure3 illustrates theconsequence of choosing a sampling frequency which is lower than the maximumfrequency. If the sampling frequency is lower than the signal to be sampled, theFourier transformed signal lies outside the spectral width (between -f/2 and +f/2). However, due to the “refolding effect” of the Fouriertransform, the signal still appears in the spectrum but at irrelevant positions.To avoid this, one may repeat the simulation by increasing the samplingfrequency and keeping other parameters constant. If the sampling frequency issufficient, the spectrum should not be affected.

The choice of the parameters discussed in this section and aboveinfluences the outcome of the simulation in the same way as it does for theexperiment. Once an NMR simulation is running, one might want to play withcombinations of f, K, taq, T2, and L until the simulated spectra display convenient features. Ifone intends to simulate spectra to match experimental data, one might simplyperform the simulation with the same f, K, L, and taq values. Table1 summarizes the parameters which were used inthis paper.

The procedure described here yields an NMR signal which is symmetricaround 0. As a consequence, each signal is found both in the positive andnegative frequencies, and the integral is split into the two duplicates. Becausethe experimental procedure that we are simulating does not differentiatenegative and positive frequencies, we discard the frequency domain signalcorresponding to negative frequencies and multiply the abscissa of the frequencydomain signal corresponding to positive frequencies by a factor of 2. Thisoperation corresponds to “folding” the spectrum around ν=0. Note that in high-field NMR, the measured signal is complexand is therefore not split into positive and negative halves. The centralfrequency of the spectrum at high field is given by the carrier frequency of thespectrometer (e.g., typically 400 MHz for ¹H at 9.4 T). Section2.8describes this difference between high-field and ZULF NMR in more detail.

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

where Irν and Icν are the real and complex frequency domain signals,respectively.

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

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

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

3. Obtain the complex frequency domain signal by Fouriertransforming the time domain signal using a fast Fourier transformalgorithm.

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

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

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

2.8. Comparison with high-field NMR

We conclude this theory section by listing the main differencesbetween high-field and ZULF NMR, which are summarized in Table2. As is the casefor the rest of the paper, our description is limited to small moleculescontaining spin1/2 in the liquid state.

At high magnetic field, the Zeeman interaction is much larger thanthe Jcoupling and therefore dominates the dynamics. Furthermore,the Larmor frequency of the spins (which results from the Zeeman interaction) isslightly shifted by the presence of the electron cloud around the nuclei. Thisphenomenon, called the chemical shift, gives a slightly different Larmorfrequency for nuclei in different positions in a molecule, which spreads overtypically 10 and 200 ppm around the Larmor frequency for¹H and¹³C spins, respectively. AtZULF, the Jcoupling dominates while the Zeeman interaction is aperturbation, and the chemical shift plays no role (in that it is a smallperturbation of a small perturbation).

In Fig.1 and in the simulations presented in this paper, we haveassumed that the detector was positioned below the sample (along the zaxis in our axis convention) and that it was sensitive tomagnetic field along the zaxis. Although this choice is typical, it is not the onlypossibility. In common high-field experiments, the oscillating signal emitted bythe spins is recorded perpendicular to the static magnetic field. Detection atZULF is performed with magnetometers that are sensitive to the total magneticfield produced by the sample. The operator corresponding to this observable isthe 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 resultsin a real signal. Note that an imaginary ZULF signal could be obtained if theOPM were to have several sensitive axes or more than one detector were used.High-field NMR uses Faraday induction in pickup coils. Signals originating fromdifferent nuclei are usually not observed in the same experiment as their Larmorfrequencies are too far apart, and the NMR coils are only sensitive over alimited bandwidth. The operator corresponding to inductive detection in pulsedNMR is non-Hermitean and therefore yields complex signals. An extra step of theacquisition process at highfield that is not required at ZULF is modulating thesignal recorded by the coil with a carrier frequency. Indeed, the NMR coil picksup a signal at the Larmor frequency of the spins, which is too high to bedigitized (e.g., 400 MHz for ¹H at 9.4 T). Instead, the signalis mixed with a carrier frequency, and only the difference is digitized, over asmall bandwidth (e.g., over 10 ppm, corresponding to 4 kHz for¹H spins at 9.4 T). Thesignals at ZULF can be detected without mixing the frequency as the they aresufficiently low to be digitalized directly. For more details on the signalmodulation at high field, the reader is referred to chapter4 of JamesKeeler'sbook Understanding NMR spectroscopy (Keeler, 2010).

The code in SupplementS3 presents in great detail the simulation ofthe spectra for a pair of J-coupled ¹H and ¹³C spin pairs at ZF and ULF andat 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 ofthe text above, and, whenever possible, the equations presented in this paperare referenced in the code. The reader is encouraged to open this code tounderstand the difference between simulating a spectrum at high field and atZULF. The code can be opened in PDF, including the figures, for those who do nothave a MATLAB license.

3.1. Excitation schemes on an XA spin system

The ZF and ULF spectra of an XA spin system with a Jcoupling of 140 Hz were simulated for different experimentalsequences, assuming that the sample consists of 100 µL of solution where the spin system has a concentration of27 mol L^-1. The code and its PDF version arepresented in SupplementS4. Figure4 shows the experimental sequences, as wellas the simulated time domain and frequency domain signals. For all simulations,the sample was assumed to have spent sufficient time in a prepolarizing field of2 T at 298 K to be at Boltzmann's equilibrium. The polarizations of the¹³C and¹H spins were calculatedusing Boltzmann's distribution (see Eq.19) and used to compute the single-spindensity matrices of the ¹³C and¹H spins, ρ^eq(13C) and ρ^eq(1H) (see Eq.21). The density matrix of the two-spin system wascomputed by taking the Kronecker product of the single-spin density matrices ρ^0=ρ^eq(13C)⊗ρ^eq(1H) (see Eq.23). The identity was removed from the two-spindensity matrix using Eq.(24). The resulting density matrix was assumed torepresent the initial state of the simulation (as explained above, only theZeeman terms are considered to contribute to the initial state). For eachexperimental sequence, the spectrum was simulated both at 0 nT (including onlythe JHamiltonian H^J; see Eq.12) and with a basis field of 0.5 µT along the xaxis, that is, orthogonal to both the direction of theprepolarizing field and the sensitive axis (including both the JHamiltonian H^J and the Zeeman Hamiltonian H^Z; see Eqs.12 and 10). The time domain signal was computed bypropagating the density matrix under the effect of the Hamiltonian for a totaltime of 5 s (parameter taq), 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 U^ for this particular time step dt (see Eq.30) and the observable operator O^ (see Eq.35) were computed only once.

The density matrix was propagated from time tk to time tk+1=tk+dt under the Hamiltonian (ZF or ULF) using the formula ρ^k+1=U^ρ^kU^-1=U^ρ^kU^† (see Eq.29). At each time point k of the propagation (realized by a “for” loop), the signalintensity of the time domain signal was extracted from the density matrix usingthe trace TrO^ρ^k (see Eq.33) (in pT). In theory, the trace of a Hermiteanoperator should be real. However, due to the finite machine precision of thenumeric algorithm, the trace can sometimes contain a nonzero imaginary part.This residual imaginary part is discarded by taking the real part of the trace ReTrO^ρ^k. This point might appear secondary, but dealing with complexnumbers while thinking they are real can lead to mistakes. After propagation, amonoexponential apodization function was applied to the time domain signal (seeEq.36), with a coherence time constant T2 of 1 s. A second apodization function was applied to avoidbaseline artifacts (see Eq.37). The apodized time domain signal was Fouriertransformed with zero-filling to 65 536 points, using MATLAB's built-infunctions. The real part of the Fourier transform is shown in Fig.4. Thefrequency axis of the spectra was computed using Eqs.(39) and (40). The spectraare symmetric around zero, and so it is common to work only with the positivefrequencies as shown in Fig.4.

Simulating the sudden field drop experiment is the simplest casepresented here. Because the coherence excitation scheme (or mixing) onlyconsists of bringing the spin from high magnetic field to ZF or ULF, thesimulation only consists of propagating the high-field thermal equilibriumdensity matrix under the ZF or ULF Hamiltonian. The ZF spectrum consists of oneline at the Jcoupling and one at zero frequency (see Fig.4a). Including abias field of 0.5 µT along the xaxis (ULF case) splits the J peak as well as the line at zero frequency.

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

which fulfills the conditions B0=Bstart and Btdecay=0. During the field drop, the Hamiltonian H^t=H^J+H^Zt is time dependent. This step thus cannot be simulated in asingle propagation step. Instead, it must be discretized into substeps dt that are sufficiently short for the Hamiltonian to beconsidered time independent. Here, the 0.5 s time interval was discretized into5000steps of 0.1 ms. At time t=0, the density matrix is the thermal equilibrium density matrix ρ^eq obtained above. At each time step tk, the propagator U^tk→tk+1=exp⁡-iH^(tk)dt is computed (see Eq.30), and the density matrix is propagatedfrom time tk to time tk+1=tk+dt (see Eq.29) under the Hamiltonian H^tk=H^J+H^Ztk. We name ρ^adia the density matrix obtained after this process. A questionarises here: is this magnetic field drop that we have chosen sufficiently slowto be considered adiabatic? In other words, is ρ^adia stationary? A simple way to ensure that it is the case is tosimulate the spectrum at ZF after the magnetic field drop without any excitationpulse, that is, taking ρ^adia as the density matrix at time t=0, ρ^0. If the transition is adiabatic, then the system should remainstationary; that is, the time domain signal should feature no oscillation andthe spectrum no peak. Figure4b shows the result of this procedure, whichconfirms that the transition is adiabatic. The only feature of the ZF spectrumin Fig.4b is the line at zero frequency. This line originates from thenon-oscillating magnetization decaying with T2, which is the result of the apodization function that we haveapplied. Verifying that the ZF spectrum is flat also ensures that the field dropwas discretized into sufficiently short time intervals dt.

The density matrix after the adiabatic field drop ρ^adia obtained above was used for the simulations presented inFig.4c–d. In the experimental sequences of Fig.4c–d, the adiabatic field dropis followed by a magnetic field pulse either along the z or xaxis. This was simulated by propagating ρ^adia under the pulse Hamiltonian to obtain ρ^0=U^pτpρ^adiaU^p†τp, where U^pτp is the propagator of the pulse Hamiltonian H^p=H^J+H^Z, which acts on the density matrix during pulse length τp. The Zeeman Hamiltonian depends on the magnetic fieldintensity of the pulse Bp 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 nextsection. The resulting density matrices ρ^0 were used as the density matrix at time t=0 of the time domain signal, which was computed and Fouriertransformed as described above. In the case of the z-axis pulse experiment, the peaks of interest (Jpeak at 140 Hz) were found to be out of phase; a phasecorrection eiϕ with ϕ=π/2 was thus applied to the Fourier transform. Adjusting the phasefor the Jpeak caused the lower-frequency peaks to be out of phase.Interestingly, in Fig.4d, the intensity of the Jpeak is higher than for the other excitation schemes while thelower-frequency peaks are suppressed, indicating that all the availablepolarization has been transferred to the Jpeak.

3.2. Rabi oscillation curves

The pairs of magnetic field intensity and length of the pulses usedfor 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 curveis the “nutation experiment”, which consists of recording a series of NMRdetections while keeping the RF pulse power constant and varying the pulselength (or the pulse length is kept constant and the amplitude is varied; Tayleret al., 2017). The nutation or Rabi curve is the plot of the signal intensity asa function of the varied parameter. It allows one to determine the pair of RFpower and pulse length which maximizes the signal intensity. Except in thepresence of rapid relaxation effects or RF field inhom*ogeneities, the observedcurve is sinusoidal. At ZULF, the Rabi curve is more complex and depends on thespin system under scrutiny. To simulate the Rabi curve at ZF, we repeated thesimulation 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 inSupplementS4). The time domain signal was Fourier transformed as describedabove, and the frequency domain signal was integrated from 138 to 142 Hz. Thesignal integral of the Jpeak is plotted as a function of the pulse length in Fig.5.The signal integral of the sudden drop experiment is shown as a horizontaldashed line for comparison. When a pulse along the zaxis is used, a simple sinusoidal curve is obtained, and itsmaximum matches that of the sudden drop experiment (see Fig.5a). The firstmaximum is reached for a pulse length of 150 µs. When a pulse along the xaxis is used, a more complex pattern is obtained, and themaximum is found to be 1.64 times higher than the sudden drop experiment (seeFig.5b). The first global maximum is reached for a pulse length of 910 µs.

3.3. XA_n spin system

The simulations shown up to this point only deal with an XA spinsystem, which typically corresponds to ¹³C-formate (or¹³C-formic acid), where the¹³C spin interacts with asingle ¹H through a Jcoupling of 195–222 Hz (Blanchard and Budker, 2016; Tayler etal., 2017) (depending on experimental conditions). ¹³C,¹⁵N-cyanide groups are alsointeresting two-spin systems which were used in ZULF experiments (Blanchard etal., 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 moleculescontaining methyl groups such as ¹³C-pyruvate (Barskiy et al.,2019). XA₄ (for example ¹⁵N-ammonium cation; Barskiy etal., 2019) and XA₅ are less common, but they are presentedhere to show the pattern that arises when adding spins.

Figure 6 shows the simulated spectra for sudden drop experimentswith detection at ZF and ULF ofXA_n spin systems with n=1,2,…,5, where X represents a ¹³C spin, andA_n represents¹H spins with a Jcoupling of 140 Hz between X and A spins and 10 Hz among Aspins (the code and its PDF version are presented in SupplementS5). All therelevant mathematics to construct the operators of an m=n+1 spin system are given in the Theory section. For anXA₅ spin system, the Hilbert space has 26=64 dimensions (and related operators). To avoid constructing eachoperator manually, recursive formulae were used (see Eqs.13 and 23). The timedomain signal was computed by propagating the density matrix under the effect ofthe Hamiltonian for a total time of 5 s (parameter taq), 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 thetime domain signal, with a coherence time constant T2 of 1 s. The apodized time domain signal was Fouriertransformed with zero-filling to 32 768 points.

Increasing the number of A spins increases the number of spectralcomponents in the spectrum. A known result of ZULF NMR appears in thissimulation: for odd numbers of n, the ZF spectrum features lines at integer multiples of the Jcoupling k⋅JAX with k∈[[1;(n+1)/2]], while for even numbers of n, it features lines at half-integer multiples of the Jcoupling k⋅JAX/2 with k∈[[1;n/2]]. Adding a 0.5 µT bias field along the xaxis during detection (that is, performing ULF detection)splits the Jlines. The higher the multiple of the Jline, the greater the number of splittings. Note that theintensity of NMR signals at highfield increases upon adding more equivalentspins to the spin system. The analysis of Fig.6 shows that this logic does notapply to the Jlines for the ZULF case, where the spectrum completely changesupon changes in the spin topology. For example, note that the amplitude of the Jline for the XA system has the same intensity as the Jline for the XA₂ system (appearing at 3/2⋅JAX frequency). Likewise, the two Jlines for the XA₃ system have the same total intensity asthe two Jlines for the XA₄ system. An empirical law ofconservation of the total spectral intensity for the Jlines can be deduced by looking at Fig.6: indeed, the totalintensity of all Jlines is the same for anyXA_n system, assuming equal sample volume,prepolarization, etc. On the other hand, the intensity of low-frequency peaksshown 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 areassociated with the precession of total magnetization around residual ULF field,and total magnetization is proportional to the number of spins.

F^2=F^x2+F^y2+F^z2,F^μ=∑l=1nI^lμ;μ∈x,y,z,H^J,F^2=H^JF^2-F^2H^J=0,

43

4.1. Eigenstates at zero field

The eigenstates of a F^2 operator can be expressed in terms of the total spin and itsprojection quantum numbers. The conventional way to express them is to use the F,mF notation, where F denotes the total spin, and mF denotes the projection onto a quantization axis (mF∈-F,-F+1,…,F-1,F). For example, by definition, for a single spin1/2, we have the sates α≡1/2,1/2andβ≡1/2,-1/2. For a pair of spins, we have the three triplet states T+1≡1,1;T0≡1,0;andT-1≡1,-1; and the singlet state S0≡0,0. Any F,mF state is an eigenstate of the F^2 and F^z operators with the following eigenvalues:

To find the total spin of a system constituted by n spins, one must sum up the angular momenta of the individualspins, which is a common procedure in the field of atomic physics but not somuch in NMR. All possible values of the angular momentum of the interactingspins are added up to constitute a set of uncoupled quasiparticles withdifferent 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:

For a pair of spins1/2, the possible values are F=0,1. For n spins, the summation should proceed until all the possiblepairs of the angular momenta of the individual spins are summed up. As anillustration, consider a coupled system of three spins1/2 (see Fig.7). First, any two spins are added up together togive F=1 (a triplet) and F=0 (a singlet). Then, the remaining spin1/2 is added up to the quasiparticles formed in the previous step(spins 1 and 0 in this case). As a result, the initialA₃ system is decomposed into threesubsystems with total spins of F=3/2, 1/2 (addition of 1 and 1/2), and 1/2 (addition of 0 and 1/2).

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

One can check that the total number of the spin states remains thesame after the procedure of adding up the spins. On the one hand, the number ofstates formed by n coupled spins I equals to 2I+1n, which is 8 in the considered case. On the other hand, amanifold with a total spin F has 2F+1 different states associated with different possibleprojections of the spin on the quantization axis. Therefore, there are 4+2+2 states in the considered case.

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

where CI,mI,J,mJF,mF represents Clebsch–Gordan coefficients, which are defined by

Each Clebsch–Gordan coefficient is specified by six numbers: thetotal spin of the coupled state F, its projection mF, and the total spins of the uncoupled states and theirprojections (I, S, mI, mS). Coefficient CI,mI,S,mSF,mF represents “how much” of uncoupled state I,mI,S,mS there is in a coupled state F,mF. The analytical values of the Clebsch–Gordan coefficients canbe calculated using recursive expressions, which are available in many softwarepackages and textbooks. TableS1.1 in the Supplement provides the relationbetween the coupled and uncoupled states for the consideredA₃ system and shows explicitly how tocalculate them. The full set of all possible F,mF states forms the new basis that is better suited than theZeeman basis for ZULF NMR. In fact, this basis coincides with the eigenstates atZULF for A_n and forXA_n systems, but this basis is also a goodstarting point for more complicated cases. We will refer to this new basis asthe “coupled” basis, because it is appropriate for the description of stronglycoupled spins.

4.2. Eigenenergies at zero field

Having the eigenstates, we can now proceed with finding theeigenvalues of the Hamiltonian; these values correspond to the energy of thestates and therefore determine the frequencies of ZULF NMR transitions. It turnsout that A_n systems are not detectable at ZULF; itis shown in the next section (where intensities of transitions are calculated)that they give rise to no observable transition. At least two types of nucleiwith different gyromagnetic ratios are necessary for an observable transition toexist. Therefore, we consider anXA_n system from now on. We will denote theoperators associated with the X spin as S^ and with A spins as I^. It is also convenient to introduce total spin operators for Aspins: F^Aμ=∑l=1nI^lμ, μ∈{x,y,z}. The Hamiltonian at ZF for this spin system is given by

The H^AX Hamiltonian can be expressed in terms of the total spinoperators using algebraic tricks. We find an expression for the first term ofEq.(48) in terms of F^2, F^A2, and S^2 by developing F^2:

Similarly, we find an expression for the second term of Eq.(48)in terms of F^A2 and I^l2 by developing F^A2:

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 easilycalculated:

The H^AX Hamiltonian commutes with the F^2 operator; therefore, they share eigenstates F,mF. So, the eigenenergies can be written as the expectationvalues of F,mF with respect to H^AX:

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

to obtain the final expression for the energy of level F,mF.

expressed in hertz. Here, quantum number F corresponds to the total spin of the fullXA_n system, S corresponds to the spin of the nucleus X, FA is the total spin of theA_n spins, and Il is the spin of individual nuclei A. The energy does not dependon 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 1/2 for ¹³C); similarly, all spinnumbers Il are the same (e.g., for ¹H spins they are equal to 1/2). The remaining two quantum numbers F and FA can have different values depending on the state, thereforeremoving degeneracy between some of the levels. Figure8 presents the energylevels of XA, XA₂, and XA₃ systems at ZF calculated usingEq.(54). Mathematica codes to perform these calculations are available in theSupplement (SupplementS6).

4.3. Selection rules

We have now found the eigenstates and their energies, but not alltransitions between any pair of states are allowed. The last step is to find thetransition intensities and thus get the analytical appearance for the ZF NMRspectrum of an XA_n system. There are certain selectionrules specifying which transitions are in principle possible and which areforbidden, like those in high-field NMR, where only single quantum transitionsare allowed. A general expression for the transition intensity between any twoeigenstates F,mF and F′,m′F is given by

We will explicitly calculate the transition intensity for thesudden field drop experiment. In this case, both the initial density operator ρ^0 and the observation operator O^ (see Eq.35) are proportional to γIF^A,z+γSS^z (as a reminder, F^A,z=∑l=1nI^l,z). Therefore, the transition intensity becomes

This expression is an example of Fermi's golden rule, which isused to calculate transition amplitudes in different problems in quantummechanics. Similar expressions can be found for the high-field NMR. Byexpressing the coupled states F,mF in terms of the uncoupled basis (see Eq.46), we find that

where mA and mS are the zprojection of the total spins FA (for n protons, the maximum value of FA equals n/2, and for each value of FA, mA∈{-FA,-FA+1,…,FA-1,FA}) and the zprojection of the spin S (in the case where S is a carbon-13 spin, mS∈{-12,12}), respectively. Now let us express the remaining F′,m′F bra in terms of the uncoupled basis as well and combineEqs.(56) and (57).