Math

The Hamiltonian $\hat{H}$ we are using for the NMR systems is of the form

$$\begin{array}{rl}\hat{H}& =-\sum _{\ell }{\gamma }_{\ell }{B}^z\left(1+{\delta }_{\ell }\right){\hat{I}}_{\ell }^z+2\pi \sum _{k<l}{J}_{kl}{\hat{\mathbf{I}}}_k\cdot {\hat{\mathbf{I}}}_l\end{array}$$

where ${\gamma }_{\ell }$ are the gyromagnetic factors and ${\delta }_{\ell }$ the chemical shifts of nuclear spin $\ell$. $J_{kl}$ denotes the coupling between spins $k$ and $l$ and ${\hat{I}}^{\alpha }={\hat{S}}^{\alpha }/\hslash$, where ${\hat{S}}^{\alpha }$ are the usual $SU(2)$ spin operators ($\alpha \in \{x,y,z\}$).

Note in the following, whenever we denote a spin operator as ${\hat{I}}^{\alpha }$ without an additional site index, we mean the sum over the individual spin operators for each nuclear spin ${\hat{I}}^{\alpha }={\sum }_i{\hat{I}}_i^{\alpha }$.

Within NMR we have a strong magnetic field $B^z$ in the $z$-direction, and electromagnetic pulses / oscillating fields are applied to flip the spins into the $x/y$-plane. $B^z$ is typically of the order of 500 Mhz, the bandwidth of the pulses is about 10 kHz, and the required resolution is less than 1 Hz.

The spectrum measured in an NMR experiment corresponds to the spectral function calculated in the HQS NMR Tool.

Calculation of the spectrum

Time domain

To discuss how to calculate a NMR spectrum let us first elaborate on how it is obtained experimentally. For this, after a probe has been placed inside a NMR spectrometer a strong magnetic field is applied along the $z$-axis, which leads to a net magnetization along the same axis. Note that at room temperature this magnetization is typically quite small. One then applies a $\pi /2$-pulse along the $x$-axis which flips the magnetization to the minus $y$-direction. The spins then start to precess in the $x/y$-plane, which creates a signal that can be recorded by measuring the magnetic field along the $x$- or $y$-axis.

Therefore, to calculate the spectrum, we simply have to calculate the time-dependent expectation value along the $x$- or $y$-direction. Here we will only consider the field along the $x$-direction. The time-dependent expectation value can be calculated using the density matrix $\rho (t)$ and the full spin operator ${\hat{I}}^x$

$$⟨{\hat{I}}^x⟩(t)=⟨{\hat{I}}^x\rho (t)⟩$$

This implies that we need to find an expression for the time evolution of the density matrix, which is described by the Liouville-von-Neumann equation

$$\frac{d\rho }{dt}=\frac{i}{\hslash }[H,\rho ]$$

In the following we will assume a Hamiltonian in units of radians per second, which allows us to drop the $\hslash$ factor. Therefore the time evolution can be written as

$$\rho (t)=e^{-iHt}\rho (t_0)e^{iHt}$$

Note that $H$ is the Hamiltonian of the process taking place within the time interval $[t_0,t]$. Therefore we have to perform two time evolutions: We start with the density matrix at some time $t_0$ and perform a time evolution using the Hamiltonian associated with the $\pi /2$ pulse. Then we perform a second time evolution using the NMR Hamiltonian, describing the precession of the spins in the $x/y$-plane.

Due to the setup of the experiment we can assume a hard (or instantaneous) pulse, meaning we will not be interested in the explicit time evolution during the pulse and can therefore write the time evolution operator of the pulse as

$$P_{\pi /2}^x=e^{-i\frac{\pi }{2}{\hat{I}}^x}$$

This leads to the following density matrix after the pulse

$$\rho (t_p)=P_{\pi /2}^x\rho (t_0)(P_{\pi /2}^x{)}^{†}$$

where $t_p$ is the time duration of the pulse. Performing the second time evolution under the NMR Hamiltonian we arrive at

$$\rho (t)=e^{-iHt}\rho (t_p)e^{iHt}=e^{-iHt}{P}_{\pi /2}^x\rho (t_0)(P_{\pi /2}^x{)}^{†}{e}^{iHt}$$

We now just have to define the density matrix at time $t_0$. Due to the experimental setup it is simply given via the Boltzmann distribution

$$\rho (t_0)=\frac{e^{-\beta H}}{Z}$$

with the inverse temperature $\beta =\frac{1}{k_{b}T}$, where $k_b$ is the Boltzmann constant and $T$ the temperature, as well as the partition function $Z$.

While HQS NMR tool provides a solver that evaluates a spectrum using the full density matrix, one should note that for typical NMR experiments it is sufficient to approximate it using a series expansion

$$\rho (t_0)=\frac{e^{-\beta H}}{Z}\approx \frac{1}{Z}(I-\beta H)\approx \frac{\beta B^z}{Z}(I-{\hat{I}}^z)$$

In the second approximation applied here, the first term with the identity $I$ can be dropped as it does not contribute to the overall expectation value and the second term is approximated with $\beta B^z{\hat{I}}^z$. This is valid as the contribution from the J-coupling term to the distribution is insignificant compared to that of the magnetic field term.

Note that in the NMR literature this step is often skipped and sometimes even ${\hat{I}}^z$ is called a density matrix. This is technically incorrect, as ${\hat{I}}^z$ is not positive semi-definite, a key requirement for a density matrix. Furthermore, the prefactor $\frac{\beta B^z}{Z}$ is often omitted, as it only serves as normalization.

However, with this approximation one can now further simplify the calculation using the identity

$$P_{\alpha }^x{\hat{I}}^z(P_{\alpha }^x{)}^{†}=e^{-i\alpha {\hat{I}}^x}{\hat{I}}^z{e}^{i\alpha {\hat{I}}^x}=\cos (\alpha ){\hat{I}}^z+i\sin (\alpha )[{\hat{I}}^z,{\hat{I}}^x]$$

which in our case of a $\pi /2$-pulse reduces to

$$P_{\pi /2}^x{\hat{I}}^z(P_{\pi /2}^x{)}^{†}=-{\hat{I}}^y$$

This leads us to the final expression for the expectation value

$$\begin{array}{rl}⟨{\hat{I}}^x⟩(t)& =⟨{\hat{I}}^x\rho (t)⟩\\ & =⟨{\hat{I}}^x{e}^{-iHt}{P}_{\pi /2}^x\rho (t_0)(P_{\pi /2}^x{)}^{†}{e}^{iHt}⟩\\ & \approx \frac{⟨{\hat{I}}^x{e}^{-iHt}{P}_{\pi /2}^x(1-\beta B^z{\hat{I}}^z)(P_{\pi /2}^x{)}^{†}{e}^{iHt}⟩}{Z}\\ & =\frac{\beta B^z}{Z}⟨{\hat{I}}^x{e}^{-iHt}{\hat{I}}^y{e}^{iHt}⟩\end{array}$$

Frequency domain

Up to now we have only discussed how to calculate the time signal, however what we are actually interested in is the spectrum given by its Fourier transform. Most implementations therefore first calculate a time evolution according to the derivation performed above and then do a Fourier transform of the result. However, as will be explained in the next section, it is much more efficient to do the Fourier transform analytically first and then perform all calculations in the frequency domain. To be able to do this, we have to rewrite the above formulation into a Lehmann representation using sums over the eigenbasis.

$$\begin{array}{rl}⟨{\hat{I}}^x⟩(t)& =\frac{\beta B^z}{Z}⟨{\hat{I}}^x{e}^{-iHt}{\hat{I}}^y{e}^{iHt}⟩\\ & =\frac{\beta B^z}{Z}\sum _n⟨n|{\hat{I}}^x\sum _m|m⟩⟨m|e^{-iHt}{\hat{I}}^y{e}^{iHt}|n⟩\\ & =\frac{\beta B^z}{Z}\sum _{n,m}⟨n|{\hat{I}}^x|m⟩⟨m|{\hat{I}}^y|n⟩e^{i(E_n-E_m)t}\end{array}$$

By introducing a small convergence aiding factor $\eta$ we can now do the Fourier transform.

$$\begin{array}{rl}⟨{\hat{I}}^x⟩(\omega )& =\int e^{i(\omega +i\eta )t}⟨{\hat{I}}^x⟩(t)dt\\ & =\frac{\beta B^z}{Z}\sum _{n,m}⟨n|{\hat{I}}^x|m⟩⟨m|{\hat{I}}^y|n⟩\int e^{i(\omega +i\eta +E_n-E_m)t}dt\\ & =\frac{\beta B^z}{Z}\sum _{n,m}\frac{⟨n|{\hat{I}}^x|m⟩⟨m|{\hat{I}}^y|n⟩}{i(\omega +i\eta +E_n-E_m)}\\ & =-i\frac{\beta B^z}{Z}\sum _{n,m}\frac{⟨n|{\hat{I}}^x|m⟩⟨m|{\hat{I}}^y|n⟩}{\omega +i\eta +E_n-E_m}\end{array}$$

This represents the final form we are evaluating. Note that the factor $\eta$ leads to a broadening of the peaks, giving them generally the shape of a Lorentzian. This broadening has to be chosen to fit the expected resolution, as discussed in the next section.

Energy rediscretization

One problem that arises in NMR is that the peaks are usually very sharp. Therefore, performing calculations in time domain typically requires a lot of time steps, that have to be evaluated. Also in the frequency domain, discretizing the frequency axis with a very fine grid would be computationally ineffective. Therefore, we perform several rounds of computing the spectral function. We start with a linear grid in frequency space and a rather large artificial broadening $\eta$. Then we iteratively rediscretize the frequency space in equal weight partitions of the total spectral function while reducing $\eta$ in each step until we reach the desired broadening. In this way, we obtain a non-linear adaptive grid with accumulation points at the spectral function peaks.