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 the ${\hat{S}}^{\alpha }$ are the usual $SU(2)$ spin operators with $\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.

Measurement of a 1D NMR spectrum

A 1D NMR spectrum can be obtained using a pulse-acquire experiment. To discuss how to calculate it we will first elaborate on how the experiment is performed. 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. The induced magnetic field ${\mathbf{B}}^{x/y}$ is directly linked to the Magnetization $\mathbf{M}$ of the sample as

$$\begin{array}{r}{\mathbf{B}}^{x/y}={\mu }_0{\mathbf{M}}^{x/y}(t)\end{array}$$

which in turn can be linked to the magnetic moment $\mathbf{m}$

$$\begin{array}{r}{\mathbf{m}}^{x/y}(t)=\int {\mathbf{M}}^{x/y}(t)dV={\mathbf{M}}^{x/y}(t)V\end{array}$$

where the integral goes over the volume $V$ of the sample and we are assuming a uniform magnetization within the sample to evaluate the integral. At the same time, the magnetic moment of an individual molecule is given as

$$\begin{array}{r}{\mathbf{m}}^{x/y}(t)=⟨\sum _k{\gamma }_k{\hat{\mathbf{S}}}_k^{x/y}⟩(t)\end{array}$$

Therefore, to calculate the spectrum, we simply have to calculate the time-dependent expectation value of the nuclear spin operators of each spin along the $x$- or $y$-direction and scale them with the corresponding gyromagnetic ratio. Note that in the actual implementation we perform the calculation using the dimensionless nuclear spin operator, so we explicitly calculate

$$\begin{array}{r}⟨{\gamma }_k{\widehat{I_k}}^{x/y}⟩(t)={\hslash }^{-1}⟨\sum _k{\gamma }_k{\hat{\mathbf{S}}}_k^{x/y}⟩(t)\end{array}$$

The sum of these contributions is directly proportional to the measured magnetic field. However, to obtain the absolute measured magnetic field, an additional scaling is required. This is left to the user, as within HQS NMR Tool we do not have information on the sample volume and concentration of the molecule. Often though, this scaling is unnecessary as one is not interested in absolute values. So in HQS NMR Tool the spectrum is by default normalized to integrate to the number of relevant spins, i.e., to those of the same isotope type as the reference isotope.

Calculation of the spectrum

Time domain

Let us now discuss, how to evaluate these expectation values. For this we will only consider the field along the $x$-direction. The time-dependent expectation value can be calculated using the density matrix $\rho (t)$

$$\begin{array}{r}⟨{\gamma }_k{\widehat{I_k}}^x⟩(t)=⟨{\gamma }_k{\hat{I}}_k^x\rho (t)⟩\end{array}$$

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

$$\begin{array}{r}\frac{d\rho }{dt}=\frac{i}{\hslash }[H,\rho ]\end{array}$$

In this formulation the equation assumes a Hamiltonian in units of Joule, however the Hamiltonian as defined above is in units of radians per second, which allows us to drop the $\hslash$ factor. Therefore the time evolution can be written as

$$\begin{array}{r}\rho (t)=e^{-iHt}\rho (t_0)e^{iHt}\end{array}$$

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

$$\begin{array}{r}{P}_{\pi /2}^x=e^{-i\frac{\pi }{2}{\hat{I}}^x}\end{array}$$

This leads to the following density matrix after the pulse

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

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

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

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

$$\begin{array}{r}\rho (t_0)=\frac{e^{-\beta \hslash H}}{Z}\end{array}$$

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$. Note that the factor $\hslash$ is necessary, as the Hamiltonian is given in units of radians per second.

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

$$\begin{array}{r}\rho (t_0)=\frac{e^{-\beta \hslash H}}{Z}\approx \frac{1}{Z}(I-\beta \hslash H)\approx \frac{1}{Z}-\frac{\beta \hslash B^z}{Z}\sum _l{\gamma }_l(1+{\delta }_l){\hat{I}}_l^z\end{array}$$

In the second approximation applied here, the first term with the identity $I$ can be dropped in the following calculations as it does not contribute to the overall expectation value and the second term is approximated with the magnetic field term of the Hamiltonian. This is valid as the contribution from the J-coupling term to the distribution is insignificant in comparison.

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 \hslash 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

$$\begin{array}{r}{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]\end{array}$$

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

$$\begin{array}{r}{P}_{\pi /2}^x{\hat{I}}^z(P_{\pi /2}^x{)}^{†}=-{\hat{I}}^y\end{array}$$

This leads us to the final expression for the expectation value

$$\begin{array}{rl}⟨{\gamma }_k{\widehat{I_k}}^x⟩(t)& =⟨{\gamma }_k{\widehat{I_k}}^x\rho (t)⟩\\ & =⟨{\gamma }_k{\widehat{I_k}}^x{e}^{-iHt}{P}_{\pi /2}^x\rho (t_0)(P_{\pi /2}^x{)}^{†}{e}^{iHt}⟩\\ & \approx \frac{⟨{\gamma }_k{\widehat{I_k}}^x{e}^{-iHt}{P}_{\pi /2}^x(1-\beta \hslash B^z{\sum }_l{\gamma }_l(1+{\delta }_l){\hat{I}}^z)(P_{\pi /2}^x{)}^{†}{e}^{iHt}⟩}{Z}\\ & =\frac{\beta \hslash B^z}{Z}{\gamma }_k\sum _l{\gamma }_l(1+{\delta }_l)⟨{\widehat{I_k}}^x{e}^{-iHt}{\hat{I}}_l^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}⟨{\gamma }_k{\widehat{I_k}}^x⟩(t)& =\frac{\beta \hslash B^z}{Z}{\gamma }_k\sum _l{\gamma }_l(1+{\delta }_l)⟨{\hat{I}}_k^x{e}^{-iHt}{\hat{I}}_l^y{e}^{iHt}⟩\\ & =\frac{\beta \hslash B^z}{Z}{\gamma }_k\sum _l{\gamma }_l(1+{\delta }_l)\sum _n⟨n|{\hat{I}}_k^x\sum _m|m⟩⟨m|e^{-iHt}{\hat{I}}_l^y{e}^{iHt}|n⟩\\ & =\frac{\beta \hslash B^z}{Z}{\gamma }_k\sum _l{\gamma }_l(1+{\delta }_l)\sum _{n,m}⟨n|{\hat{I}}_k^x|m⟩⟨m|{\hat{I}}_l^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}⟨{\gamma }_k{\hat{I}}_k^x⟩(\omega )& =\int e^{i(\omega +i\eta )t}⟨{\gamma }_k{\hat{I}}_k^x⟩(t)dt\\ & =\frac{\beta \hslash B^z}{Z}{\gamma }_k\sum _l{\gamma }_l(1+{\delta }_l)\sum _{n,m}⟨n|{\hat{I}}_k^x|m⟩⟨m|{\hat{I}}_l^y|n⟩\int e^{i(\omega +i\eta +E_n-E_m)t}dt\\ & =\frac{\beta \hslash B^z}{Z}{\gamma }_k\sum _l{\gamma }_l(1+{\delta }_l)\sum _{n,m}\frac{⟨n|{\hat{I}}_k^x|m⟩⟨m|{\hat{I}}_l^y|n⟩}{i(\omega +i\eta +E_n-E_m)}\\ & =-i\frac{\beta \hslash B^z}{Z}{\gamma }_k\sum _l{\gamma }_l(1+{\delta }_l)\sum _{n,m}\frac{⟨n|{\hat{I}}_k^x|m⟩⟨m|{\hat{I}}_l^y|n⟩}{\omega +i\eta +E_n-E_m}\end{array}$$

This represents the final form we are evaluating although for the spectrum we are measuring in NMR typically just the real part is required. However for some postprocessing routines also the imaginary part is required, in which case we will refer to this as a Green's function.

Additionally, 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.