Nuclear magnetic resonance and electron paramagnetic resonance: Basics

Introduction

Nuclear Magnetic Resonance (NMR) and Electron Paramagnetic Resonance (EPR) are powerful and versatile techniques that have revolutionized the fields of chemistry, physics, and medicine. These techniques exploit the intrinsic magnetic properties of atomic nuclei and unpaired electrons, respectively, to reveal detailed information about the structure, dynamics, and electronic environment of molecules.

When placed in an external magnetic field, certain nuclei—such as hydrogen or carbon in NMR—or unpaired electrons in EPR, behave like tiny gyroscopes, precessing around the magnetic field at a characteristic frequency known as the Larmor frequency. This frequency is directly related to the magnetic field strength and the specific particle being observed, whether it be a nucleus or an electron.

In NMR, the interaction between spinning nuclei and an applied external magnetic field leads to Zeeman splitting, where the energy levels of the nuclei split according to their spin states. Similarly, in EPR, unpaired electrons experience energy level splitting in a magnetic field. By introducing a radio frequency (RF) pulse (in NMR) or a microwave pulse (in EPR), we can perturb these spin systems, leading to a rich array of resonance phenomena that can be detected and analyzed.

This article delves deep into both the classical and quantum mechanical descriptions of NMR, beginning with the phenomenological Bloch equations, which govern the time evolution of nuclear magnetization under the influence of static and oscillating magnetic fields. We explore the concepts of longitudinal ($T_1$) and transverse ($T_2$) relaxation times, which characterize how the system returns to equilibrium after being disturbed, as well as the effects of a classical driving field on spin polarization.

Building on the classical foundation, we then transition to the quantum mechanical description of NMR and EPR, examining the system dynamics through the lens of quantum operators, the density matrix formalism, and the Lindblad master equation. This comprehensive approach allows us to recover the classical Bloch equations from the quantum mechanical model, providing a deep understanding of the microscopic origins of macroscopic observations.

Through a combination of theoretical models and simulations, this article aims to provide a thorough understanding of how NMR and EPR work and why these techniques remain cornerstones in both research and industry.

Update: Since this article only covers the theory there is now an companion article showing you some real world measurements with a small DIY NMR spectrometer to illustrate how real measurements actually look like.

• Introduction
• Nuclear Zeeman effect - the two level system
• Phenomenological Bloch equations
• Solutions to the Bloch equations
  • A steady state solution without relaxation
  • Longitudinal relaxation ($T_1$) and constant external field
  • Transversal relaxation ($T_2$) and constant external field
  • Bloch equations with an external driving field
    • Transforming vectors into rotating fields, differentials of rotation matrices
    • Bloch equation without relaxation with driving field
    • Bloch equation with driving field and relaxation
  • Numerical simulations
    • Driving with an on resonant field
    • Decay due to relaxation effects
• The quantum mechanical description
  • A short reminder on spin operators and Pauli matrices
  • Our Hamiltonian and its Eigenvalues and Eigenvectors
  • Density matrices and their time evolution
  • Initial state of our system
  • Adding a driving field
  • Moving to the rotating frame, rotating frame approximation
  • On resonance solutions and Rabi oscillations
  • Introducing relaxation processes
    • Longitudinal relaxation T1
    • Transverse relaxation T2
    • Complete equation including relaxation
  • Expectation values of magnetization vectors
    • Unitary part (Hamiltonian)
    • Relaxation (Lindblad terms)
      • Longitudinal relaxation (T1)
      • Transverse relaxation (T2)
    • Unitary and non unitary contributions
• Linewidth and excited state lifetime
  • Limits by transversal relaxation
  • Associated lifetime limit by transversal relaxation

Nuclear Zeeman effect - the two level system

The energy splitting in nuclear states exploited by NMR is caused by the coupling of the nuclear spins to an external magnetic field. Depending on their orientation they have different energy states. From a quantum mechanical view this is in it’s most simple configuration a simple two level system with all of it’s consequences - the quantum mechanical view will predict and explain effects like Rabi oscillations and provide a more in depth view of NMR. In this article we’re going to focus on the classical phenomenological view that utilizes the classical coupling of a magnetic gyroscope to the polarizing field - depending on it’s orientation there can be different energy states.

In thermal equilibrium the population imbalance by the two states (the ground state $<\alpha\mid$ and the excited state $<\beta\mid$ is determined by the temperature of the system according to the Boltzmann law:

[ \begin{aligned} \Delta E &= \hbar \gamma B \\ \frac{N_1}{N_2} &= e^{-\frac{\Delta E}{k_B T}} \end{aligned} ]

The net magnetization caused by the difference between higher and lower state that can be detected is thus described by:

[ \begin{aligned} \Delta N &= N_{lower} - N_{upper} \\ &= N \left( 1 - e^{-\frac{\Delta E}{k_B T}} \right) \end{aligned} ]

As one can see one can influence the population difference either by increasing the energy gap (and thus going to higher magnetic fields - which also yields higher resonance frequencies as we’re going to see shortly) or by massively decreasing the temperature - which is the reason why one sees many NMR and EPR spectrometers in science operating at cryogenic temperatures of liquid nitrogen (77K), liquid helium (4K) or even below.

Phenomenological Bloch equations

In the following sections we’re going to look at a phenomenological non quantum mechanical description of NMR. This view is described by the phenomenological Bloch equations.

They describe the net magnetization in a sample when applying an external magnetic field. The idea starts off with looking at a single spin precessing around a quantization axis set by the external static polarizing magnetic field $B_0$, usually assumed along the $z$ axis Since there is no net force acting on the dipole it’s assumed that the torque will be of interest:

[ \begin{aligned} \frac{\partial \vec{L}}{\partial t} &= \sum{\vec{\tau}} \\ &= \vec{\mu} \times \vec{B} \end{aligned} ]

The magnetic momentum of a spin can be related with the angular momentum using the gyromagnetic ratio:

[ \frac{\partial \vec{\mu}}{\partial t} = \gamma \vec{\mu} \times \vec{B} ]

Assuming that the phenomenological net magnetization is resulting of the sum of all magnetic momenta this yields to:

[ \begin{aligned} \vec{M} &= \sum_{i} \vec{\mu_i} \\ \frac{\partial \vec{M}(t)}{\partial t} &= \gamma \vec{M}(t) \times \vec{B}(t) \end{aligned} ]

Solutions to the Bloch equations

A steady state solution without relaxation

First let’s take a look at the steady state solution without relaxation in a static external magnetic field along the $z$ axis without applying any driving fields. We’re going to see the spins are precessing around the quantization axis.

[ \begin{aligned} \vec{B} &= \begin{pmatrix} 0 \\ 0 \\ B_z \end{pmatrix} \\ \frac{\partial}{\partial t} \begin{pmatrix} M_x \\ M_y \\ M_z \end{pmatrix} &= \gamma \begin{pmatrix} M_x \\ M_y \\ M_z \end{pmatrix} \times \begin{pmatrix} 0 \\ 0 \\ B_z \end{pmatrix} \end{aligned} ]

This yields a very simple set of differential equations when being read component wise:

[ \begin{aligned} \partial_t M_x &= \gamma ( M_y * B_z - M_z * 0 ) \\ \partial_t M_y &= \gamma ( -M_x * B_z + M_z * 0 ) \\ \partial_t M_z &= \gamma * 0 \end{aligned} ]

The last equation shows that without relaxation effects $M_z$ is time independent - the spin precesses around the axis.

The remaining two differential equations are simple to solve assuming that $\partial_t B_z = 0$ (and also assuming the gyromagnetic ratio is constant over time). We do this by introducing the quantity

[ M_{xy} = M_{x} + i M_{y} ]

This will introduce a transversal co-rotating field component. We’ll see this shortly. Adding the two equations with an additional $i$ term:

[ \begin{aligned} \partial_t M_x &= \gamma ( M_y * B_z - M_z * 0 ) \\ \partial_t M_y &= \gamma ( -M_x * B_z + M_z * 0 ) \\ \to \partial_t M_x + i * \partial_t M_y &= M_y B_z - i * M_x B_z \\ \partial_t \underbrace{(M_x + i _My)}_{M_{xy}} &= -i * (i \gamma M_y B_z + \gamma M_x B_z) \\ &= -i \gamma B_z (i M_y + M_x) \\ &= -i \gamma B_z \underbrace{(M_x + i M_y)}_{M_{xy}} \end{aligned} ]

We’ve just reduced the system of two coupled differential equations to a single simply solvable differential equation:

[ \partial_t M_{xy} = -i \gamma B_z M_{xy} ]

Using the classical Ansatz

[ \begin{aligned} M_{xy} &= M_{xy}(0) e^{-i \omega t} \\ \to \partial_t M_{xy} &= -i \omega M_{xy}(t) \end{aligned} ]

we now arrive at a very simple solution:

[ \begin{aligned} \partial_t M_{xy} &= - i \gamma B_z M_{xy} \\ \to -i \omega M_{xy}(t) &= -i \gamma B_z M_{xy} \\ \to \omega = \gamma B_z \end{aligned} ]

As we can see for the steady state solution without relaxation the $M_{xy}$ component precesses with the so called Larmor frequency $\omega = \gamma B_z$ around the bias / quantization ($z$) axis keeping a constant amplitude - and the amplitude of the magnetization in the quantization axis $M_z(t)$ also stays constant.

[ \begin{aligned} M_{xy}(t) &= M_{xy}(0) e^{i \gamma B_z t} \\ M_{z}(t) &= M_{z}(0) \end{aligned} ]

Longitudinal relaxation ($T_1$) and constant external field

Out of experimental data one can assume that net magnetization will return to rest if perturbed. Assuming that one has oriented all spins along the $z$ axis at rest and then has kicked them towards the $xy$ plane one can assume that the spins will relax into their steady state orientation with a velocity that’s proportional to the difference of the current magnetization $M_z(t)$ and the steady state magnetization $M_z(0)$ since it’s a stochastic process:

[ \begin{aligned} \frac{\partial M_z(t)}{\partial t} = \frac{1}{T_1} \left(M_0 - M_z(t)\right) \end{aligned} ]

This is a simple inhomogeneous differential equation of first order. First let’s solve the homogeneous differential equation:

[ \begin{aligned} \frac{dM(t)}{dt} &= -\frac{1}{T_1} M(t) \\ \frac{dM(t)}{M(t)} &= -\frac{1}{T_1} dt \\ \int \frac{1}{M(t)} dM(t) &= -\frac{1}{T_1} \int dt \\ ln(M(t)) &= -\frac{t}{T_1} + c_1 \\ M(t) &= e^{-\frac{t}{T_1}} \underbrace{e^{C_1}}_{c_2} \\ \to M_h(t) &= c_2 * e^{-\frac{t}{T_1}} \end{aligned} ]

Now lets search a particular solution. Since our inhomogeneity is constant ($M_0$) and not time dependent we assume that our homogeneous term will also be constant. We can insert this Ansatz into our differential equation:

[ \begin{aligned} M_p'(t) &= 0 \\ M'(t) &= -\frac{1}{T_1} M(t) + \frac{M_0}{T_1} \\ 0 &= -\frac{1}{T_1} M_p + \frac{M_0}{T_1} \\ \to M_p &= M_0 \end{aligned} ]

As usual assume that our complete solution is formed by a combination of the homogeneous solution that captures the dynamics and the particular solution that captures the initial condition:

[ \begin{aligned} M_h(t) &= c_2 e^{-\frac{t}{T_1}} \\ M_p(t) &= M_0 \\ M(t) &= M_h(t) + M_p(t) \\ M(t) &= c_2 e^{-\frac{t}{T_1}} + M_0 \end{aligned} ]

To determine the constant we use the initial condition:

[ \begin{aligned} M(0) &= M_z(0) \\ \to c_2 + M_0 &= M_z(0) \\ \to c_2 &= (M_z(0) - M_0) \\ \to M(t) &= M_z(0) e^{-\frac{t}{T_1}} + M_0 \left( 1 - e^{-\frac{t}{T_1}} \right) \end{aligned} ]

Transversal relaxation ($T_2$) and constant external field

For transversal relaxation one can assume a decay rate of $\frac{1}{T_2}$ that is proportional to the transversal field. The steady state should approach zero so the equations are simpler than for longitudinal relaxation:

[ \begin{aligned} \partial_t M_x &= \gamma B_0 M_y - \frac{M_x}{T_2} \\ \partial_t M_y &= - \gamma B_0 M_x - \frac{M_y}{T_2} \end{aligned} ]

To solve the equation one uses the trick utilizing the complex plane and adds both equations from each other again:

[ \begin{aligned} M_{xy} &= M_x + i M_y \\ \partial_t M_x + i * \partial_t M_y &= (\gamma B_0 M_y - \frac{M_x}{T_2}) + i * (- \gamma B_0 M_x - \frac{M_y}{T_2}) \\ \partial_t (M_x + i M_y) &= \gamma B_0 (M_y - i M_x) - \frac{M_x + i M_y}{T_2} \\ &= -i \gamma B_0 (i M_y + M_x) - \frac{M_x + i My}{T_2} \\ &= -i \gamma B_0 (M_x + i M_y) - \frac{M_x + i M_y}{T_2} \\ \to \partial_t M_{xy} &= -i \gamma B_0 M_{xy} - \frac{M_{xy}}{T_2} \end{aligned} ]

Now using again the Ansatz

[ \begin{aligned} M_{xy} &= c e^{-i \omega t} \\ \partial_t M_{xy} &= - i \omega M_{xy} \end{aligned} ]

We can solve the equation again:

[ \begin{aligned} -i \omega M_{xy} &= -i \gamma B_0 m_{xy} - \frac{1}{T_2} M_{xy} \\ &= \underbrace{(-i \gamma B_0 - \frac{1}{T_2})}_{-i \omega} M_{xy} \\ \to M_{xy}(t) &= c e^{(-i \gamma B_0 - \frac{1}{T_2}) t} \\ &= M_{xy}(0) e^{-i \gamma B_0 t} e^{-\frac{t}{T_2}} \end{aligned} ]

As one can see the second term yields relaxation back into the polarizing field direction. When one waits indefinitely this should yield a zero field in the transversal direction:

[ \lim_{t \to \infty} M_{xy}(t) = M_{xy}(0) e^{-i \gamma B_0 t} \underbrace{e^{- \frac{t}{T_2}}}_{\to 0} \to 0 ]

Bloch equations with an external driving field

Up until now we only looked at the Bloch equations in a steady state - with only an static external $B_z$ bias field. To drive the system one usually applies an classical radio frequency (RF) drive field in an orthogonal direction. This is usually modeled in the $xy$ plane. Under realistic circumstances one also gets components in $B_z$ direction. Then it gets even harder to solve since decoupling of the equations does not work any more.

Transforming vectors into rotating fields, differentials of rotation matrices

Now let’s first recapitulate a little math on rotation matrices. This will help in solving the equation massively. Let’s assume we rotate a vector around the $z$ axis on the $xy$ plane. We can do this by applying the rotation matrix

[ R_{z} (\omega t) = \begin{pmatrix} cos(\omega t) & -sin(\omega t) & 0 \\ sin(\omega t) & cos(\omega t) & 0 \\ 0 & 0 & 1 \end{pmatrix} ]

To transform a vector from our laboratory frame into a co-rotating frame we just multiply it with the rotation matrix:

[ \vec{v}_R = \begin{pmatrix} cos(-\Omega t) & -sin(-\Omega t) & 0 \\ sin(-\Omega t) & cos(-\Omega t) & 0 \\ 0 & 0 & 1 \end{pmatrix} * \vec{v}_{L} ]

Now lets take a closer look at the time derivative of the vector in the rotating frame:

[ \frac{\partial \vec{v}_R}{\partial _t} = \frac{\partial R_z(-\Omega t)}{\partial t} \vec{v}_L + R_z(-\Omega t) \frac{\partial \vec{v}_L}{\partial t} ]

As one can see this contains a time differential of a rotation matrix. Taking a closer look again:

[ \begin{aligned} \frac{\partial R_z(\Omega t)}{\partial t} &= \frac{\partial}{\partial t} \begin{pmatrix} cos(\omega t) & -sin(\omega t) & 0 \\ sin(\omega t) & cos(\omega t) & 0 \\ 0 & 0 & 1 \end{pmatrix} \\ &= - \Omega \begin{pmatrix} -sin(\omega t) & -sin(\omega t) & 0 \\ cos(\omega t) & -sin(\omega t) & 0 \\ 0 & 0 & 0 \end{pmatrix} \end{aligned} ]

On first glance this doesn’t look very convenient but taking a closer look on the outer product

[ \begin{aligned} \begin{pmatrix} 0 \\ 0 \\ 1 \end{pmatrix} \times R_z &= \begin{pmatrix} 0 \\ 0 \\ 1 \end{pmatrix} \times \begin{pmatrix} cos(\omega t) & -sin(\omega t) & 0 \\ sin(\omega t) & cos(\omega t) & 0 \\ 0 & 0 & 1 \end{pmatrix} \\ &= \begin{pmatrix} -sin(\omega t) & -sin(\omega t) & 0 \\ cos(\omega t) & -sin(\omega t) & 0 \\ 0 & 0 & 0 \end{pmatrix} \end{aligned} ]

One can actually rewrite this expression as

[ \frac{\partial R_z(\Omega t)}{\partial t} = -\Omega \begin{pmatrix} 0 \\ 0 \\ 1 \end{pmatrix} \times R_z(-\Omega t) ]

Bloch equation without relaxation with driving field

First lets take a look at the Bloch equations when adding a driving field when neglecting the relaxation processes. This is a reasonable approximation for many pulsed experiments since the driving pulse durations are often very short compared to relaxation times.

Assume the driving field is operating with a frequency $\omega_{RF}$ in the $x$ plane. Thus our driving field vector in the laboratory frame is

[ B_{RF} = \begin{pmatrix} B_1 \cos(\omega_{RF} t) \\ -B_1 \sin(\omega_{RF} t) \\ 0 \end{pmatrix} ]

We can now model our system again using the Bloch equation in the laboratory frame. In the following section all components in the laboratory frame are denoted by the $L$ index, all in the rotating frame by the $R$ label (i.e. $M_L$ and $M_R$).

[ \frac{\partial M_L}{\partial t} = \gamma \vec{M_L} \times \vec{B_L} ]

As we will see later it’s way easier to move into a co-rotating frame that is rotating around our quantization axis $z$ with the Lamour frequency. In the following we assume the rotation frequency of our co-rotating frame is described by $\Omega$. Thus our Bloch equation in the co-rotating frame can be derived by transforming our Magnetization vector into the co-rotating frame:

[ \frac{\partial M_R}{\partial t} = \frac{\partial (R_Z(-\Omega t) M_L)}{\partial t} ]

Now we just differentiate with the product rule and apply the differential of the rotation matrix from the last subsection:

[ \begin{aligned} \frac{\partial M_R}{\partial t} &= \frac{\partial (R_Z(-\Omega t) M_L)}{\partial t} \\ &= \frac{\partial R_z(-\Omega t)}{\partial t} M_L + R_Z(-\Omega t) \frac{\partial M_L}{\partial t} \\ &= -\Omega \begin{pmatrix} 0 \\ 0 \\ 1 \end{pmatrix} \times M_R + R_Z(-\Omega t) \frac{\partial M_L}{\partial t} \\ &= -\Omega \begin{pmatrix} 0 \\ 0 \\ 1 \end{pmatrix} \times M_R + R_Z(-\Omega t) \underbrace{\left( \gamma M_L \times \begin{pmatrix} B_1 \cos(\omega_{RF} t) \\ B_1 \sin(\omega_{RF} t) \\ B_z \end{pmatrix} \right)}_{\frac{\partial M_L}{\partial t}} \end{aligned} ]

Since the rotation matrix is acting on the whole term $\frac{\partial M_L}{\partial t}$ in the laboratory frame and it’s just composed of two vector components we can transform each of the vectors by the rotation matrix into the rotating frame. We first take a look at the magnetic field components:

[ \begin{aligned} \vec{B}_R &= \begin{pmatrix} \cos(\Omega t) & -\sin(\Omega t) & 0 \\ \sin(\Omega t) & \cos(\Omega t) & 0 \\ 0 & 0 & 1 \end{pmatrix} \begin{pmatrix} B_1 \cos(\omega_{RF} t) \\ -B_1 \sin(\omega_{RF} t) \\ B_z \end{pmatrix} \\ &= \begin{pmatrix} B_1 \cos(\omega_{RF} t) \cos(\Omega t) + B_1 \sin(\omega_{RF} t) \sin(\Omega t) \\ B_1 \cos(\omega_{RF} t) \sin(\Omega t) - B_1 \sin(\omega_{RF} t) \cos(\Omega t) \\ B_z \end{pmatrix} \\ &= \begin{pmatrix} \frac{1}{2} \left( B_1 \cos(\omega_{RF}t + \Omega t) + B_1 \cos(\omega_{RF}t - \Omega t) + B_1 \cos(\omega_{RF}t - \Omega t) - B_1 \cos(\omega_{RF}t + \Omega t) \right) \\ \frac{1}{2} \left( B_1 \sin(\omega_{RF}t + \Omega t) + B_1 \sin(-\omega_{RF}t + \Omega t) - B_1 \sin(\omega_{RF}t + \Omega t) - B_1 \sin(\omega_{RF}t - \Omega t) \right) \\ B_z \end{pmatrix} \\ &= \begin{pmatrix} B_1 \cos((\omega_{RF} - \Omega) t) \\ -B_1 \sin((\omega_{RF} - \Omega) t) \\ B_z \end{pmatrix} \end{aligned} ]

We can now insert this field into our Bloch equation in the rotating frame:

[ \begin{aligned} \frac{\partial M_R}{\partial t} &= -\Omega \begin{pmatrix} 0 \\ 0 \\ 1 \end{pmatrix} \times M_R + \left(\gamma M_R \times \begin{pmatrix} B_1 \cos((\omega_{RF} - \Omega) t) \\ -B_1 \sin((\omega_{RF} - \Omega) t) \\ B_z \end{pmatrix} \right) \\ &= - \Omega \begin{pmatrix} -M_{R,y} \\ M_{R,x} \\ 0 \end{pmatrix} + \gamma \begin{pmatrix} M_{R,x} \\ M_{R,y} \\ M_{R,z} \end{pmatrix} \times \begin{pmatrix} B_1 \cos((\omega_{RF} - \Omega)t) \\ -B_1 \sin((\omega_{RF} - \Omega) t) \\ B_z \end{pmatrix} \\ &= - \Omega \begin{pmatrix} -M_{R,y} \\ M_{R,x} \\ 0 \end{pmatrix} + \gamma \begin{pmatrix} M_{R,y} B_Z + M_{R,z} B_1 \sin((\omega_{RF} - \Omega) t) \\ -M_{R,x} B_z + M_{R,z} B_1 \cos((\omega_{RF} - \Omega)t) \\ - M_{R,x} B_1 \sin((\omega_{RF} - \Omega)t) - M_{R,y} B_1 \cos((\omega_{RF} - \Omega) t) \end{pmatrix} \end{aligned} ]

Bloch equation with driving field and relaxation

As a last step lets now add relaxation again. Recall that in the laboratory frame the phenomenological relaxation is described by the last term in the following expression:

[ \frac{\partial M_L}{\partial t} = \gamma \vec{M}_L \times \vec{B}_{eff} - \begin{pmatrix} \frac{M_{L,x}}{T_2} \\ \frac{M_{L,y}}{T_2} \\ -\frac{M_0 - M_{L,z}}{T_1} \end{pmatrix} ]

To include them in the rotating frame let’s just apply the rotation matrix:

[ \begin{pmatrix} \cos(\Omega t) & -\sin(\Omega t) & 0 \\ \sin(\Omega t) & \cos(\Omega t) & 0 \\ 0 & 0 & 1 \end{pmatrix} \begin{pmatrix} \frac{M_{L,x}}{T_2} \\ \frac{M_{L,y}}{T_2} \\ \frac{M_0 - M_{L,z}}{T_1} \end{pmatrix} = \begin{pmatrix} \frac{1}{T_2} (\cos(\Omega t) M_{L,x} - \sin(\Omega t) M_{L,y}) \\ \frac{1}{T_2} (\sin(\Omega t) M_{L,x} + \cos(\Omega t) M_{L,y}) \\ -\frac{M_0 - M_{L,z}}{T_1} \end{pmatrix} ]

Taking a look at the transformed net magnetization vector $M_L$ in the rotating frame:

[ \begin{aligned} \begin{pmatrix} M_{R,x} \\ M_{R,y} \\ M_{R,z} \end{pmatrix} &= \begin{pmatrix} \cos(\Omega t) & -\sin(\Omega t) & 0 \\ \sin(\Omega t) & \cos(\Omega t) & 0 \\ 0 & 0 & 1 \end{pmatrix} \begin{pmatrix} M_{L,x} \\ M_{L,y} \\ M_{L,z} \end{pmatrix} \\ &= \begin{pmatrix} \cos(\Omega t) M_{L,x} - \sin(\Omega t) M_{L,y} \\ \sin(\Omega t) M_{L,x} + \cos(\Omega t) M_{L,y} \\ M_{L,z} \end{pmatrix} \end{aligned} ]

This shows we can simply apply the prefactors $\frac{1}{T_2}$ and the structure $\frac{M_0 - M_{L,z}}{T_1}$ to the transformed vector as one could have expected:

[ \frac{\partial M_R}{\partial t} = - \Omega \begin{pmatrix} -M_{R,y} \\ M_{R,x} \\ 0 \end{pmatrix} + \gamma \begin{pmatrix} M_{R,y} B_Z + M_{R,z} B_1 \sin((\omega_{RF} - \Omega) t) \\ -M_{R,x} B_z + M_{R,z} B_1 \cos((\omega_{RF} - \Omega)t) \\ - M_{R,x} B_1 \sin((\omega_{RF} - \Omega)t) - M_{R,y} B_1 \cos((\omega_{RF} - \Omega) t) \end{pmatrix} - \begin{pmatrix} \frac{M_{R,x}}{T_2} \\ \frac{M_{R,y}}{T_2} \\ -\frac{M_0 - M_{R,z}}{T_1} \end{pmatrix} ]

Thus we have our set of the 3 coupled differential equations that model a simple system with quantization or polarization axis $z$ and an perpendicular RF field in the XY plane:

[ \begin{aligned} \partial_t M_{R,x}(t) &= \Omega M_{R,y}(t) + \gamma M_{R,y}(t) B_z + \gamma M_{R,z}(t) B_1 \sin((\omega_{RF} - \Omega)t) - \frac{1}{T_2} M_{R,x}(t) \\ \partial_t M_{R,y}(t) &= -\Omega M_{R,x}(t) - \gamma M_{R,x}(t) B_z + \gamma M_{R,z}(t) B_1 \cos((\omega_{RF} - \Omega) t) - \frac{1}{T_2} M_{R,y}(t) \\ \partial_t M_{R,z}(t) &= - \gamma M_{R,x}(t) B_1 \sin((\omega_{RF} - \Omega) t) - \gamma M_{R,y}(t) B_1 \cos((\omega_{RF} - \Omega) t) + \frac{M_0 - M_{R,z}(t)}{T_1} \end{aligned} ]

Numerical simulations

What do these equations actually mean? Keep in mind those are classical phenomenological equations - they don’t describe effects as Rabi oscillations - so for example you will only see saturation of the excited state with longer and longer pulse times. This does not fit the experiment where the largest population of the excited state will be achieved using $\frac{\pi}{2}$ pulses.

The following simulations have been calculated with the following parameters:

• Time range over 25 ms
• Initial magnetization along the $z$ axis
• $B_0$ field of 7 T
• $B_1$ field strength of 1 mT
• $\omega_{RF}$ of 42 MHz
• $\Omega$ of 42 MHz
• $\gamma$ - the gyromagnetic ratio - of $42 \frac{MHz}{T}$
• The $T_1$ and $T_2$ times have been chosen way shorter than real values to keep the simulation time span short - they’ve been shorted by a factor of about 1 million for a free proton
  • $T_1 = 650ns$
  • $T_2 = 150ns$

Driving with an on resonant field

To simulate the driving the initial state is a single spin polarized completely parallel to the quantization axis - into $z$ direction. The driving field, relaxation and the bias field are enabled.

As one can see the driving field in rotational frame is constant while oscillating in the lab frame (plots in the lower rows). This is due to the rotating frame rotating with the Lamour frequency. This slowly tilts the spin from the $z$ plane (the magnetization shown in the right upper plot) into the $xy$ plane (the magnetization shown in the left upper plot). As one can see this yields two oscillating components. One in phase with the driving field (i.e. along the co-rotating $x$ axis) - this is called the absorptive component since it counteracts the effective field that also drives the system. The second signal along the $y$ axis is shifted with respect to the drive signal - this is called dispersive signal due to this phase shift.

Decay due to relaxation effects

Illustrating the decay effect due to relaxation effects is done by doing the same simulation, initially polarizing the spins in the $xy$ plane and disabling the drive field.

As one can see the magnetization along the $xy$ plane decays exponentially while the magnetization relaxes back to it’s initial value along the $z$ plane as expected. Again there is an absorptive and dispersive component in the co-rotating $xy$ plane that yields an exponentially decaying absorptive and dispersive signal.

The quantum mechanical description

As we can see above the description of relaxation matches what we would expect and what we will see later in experimental data - but it fails to describe some well known effects like Rabi oscillations (that we would expect for any quantum mechanical two level system and it does not tell us how the system really behaves on microscopic scale but delivers only a macroscopic phenomenological description of the effects. To dive further into the microscopic behavior we have to move on to the quantum mechanical description of our spins. To do this we first look at a single spin. To do this lets first recall the notation of spin operators and Pauli matrices that we’re going to use to describe spins in our $SU(2)$:

A short reminder on spin operators and Pauli matrices

Let’s recall the spin operators in quantum mechanics. They are designed to behave like angular momentum operators. We usually use Pauli matrices as one of the representations of the generators of the $SU(2)$ Lie group that contains our spin state:

[ \begin{aligned} \hat{\sigma}_x &= \begin{pmatrix} 0 & 1 \\ 1 & 0 \end{pmatrix} \\ \hat{\sigma}_y &= \begin{pmatrix} 0 & -i \\ i & 0 \end{pmatrix} \\ \hat{\sigma}_z &= \begin{pmatrix} 1 & 0 \\ 0 & -1 \end{pmatrix} \end{aligned} ]

Those matrices are designed to satisfy the commutation relations that we expect from an angular momentum operator:

[ \begin{aligned} \lbrack \hat{\sigma}_x, \hat{\sigma}_y \rbrack = 2 i \hat{\sigma_z} \\ \lbrack \hat{\sigma}_y, \hat{\sigma}_z \rbrack = 2 i \hat{\sigma_x} \\ \lbrack \hat{\sigma}_z, \hat{\sigma}_x \rbrack = 2 i \hat{\sigma_y} \end{aligned} ]

The spin operators are designed using the Pauli matrices:

[ \begin{aligned} \hat{S}_x &= \frac{\hbar}{2} \hat{\sigma}_x \\ \hat{S}_y &= \frac{\hbar}{2} \hat{\sigma}_y \\ \hat{S}_z &= \frac{\hbar}{2} \hat{\sigma}_z \end{aligned} ]

Our Hamiltonian and its Eigenvalues and Eigenvectors

Now that we recalled how spin operators work lets set up a Hamiltonian for our system. This will model the coupling between a single spin and an external magnetic field:

[ \begin{aligned} \hat{H}_0 &= -\gamma \vec{B}_0 \hat{S} \\ &= -\gamma B_0 \hat{S_z} \end{aligned} ]

As one can see we have assumed again that the external bias field $B_0$ is oriented along the z Axis - so the inner product of $\vec{B}_0$ and the spin operator $\hat{S}$ only yields a $z$ component.

The Eigenstates $\mid \uparrow >$ and $\mid \downarrow >$ of the $\hat{S}_z$ operator can be derived when solving the Eigenstates of the Hamiltonian:

[ \begin{aligned} \hat{H_0} &= -\gamma B_0 \hat{S}_z \\ &= -\frac{\hbar \gamma B_0}{2} \hat{\sigma}_z \\ &= -\frac{\hbar \gamma B_0}{2} \begin{pmatrix} 1 & 0 \\ 0 & -1 \end{pmatrix} \end{aligned} ]

Let’s solve for Eigenvalues and Eigenvectors of this operator:

[ \begin{aligned} \hat{H}_0 \mid \psi > &= E \mid \psi > \\ -\frac{\hbar \gamma B_0}{2} \begin{pmatrix} 1 & 0 \\ 0 & -1 \end{pmatrix} \begin{pmatrix} a \\ b \end{pmatrix} &= E \begin{pmatrix} a \\ b \end{pmatrix} \end{aligned} ]

Using the method of characteristic polynomials we can determine our Eigenvalues:

[ A = \begin{pmatrix} -\frac{\hbar \gamma B_0}{2} & 0 \\ 0 & \frac{\hbar \gamma B_0}{2} \end{pmatrix} \\ ] [ \begin{aligned} \det(\lambda I - A) &= 0 \\ \det\left(\begin{pmatrix} \lambda & 0 \\ 0 & \lambda \end{pmatrix} - \begin{pmatrix} -\frac{\hbar \gamma B_0}{2} & 0 \\ 0 & \frac{\hbar \gamma B_0}{2} \end{pmatrix}\right) &= 0 \\ \det \begin{pmatrix} \lambda + \frac{\hbar \gamma B_0}{2} & 0 \\ 0 & \lambda - \frac{\hbar \gamma B_0}{2} \end{pmatrix} &= 0 \\ \left(\lambda + \frac{\hbar \gamma B_0}{2}\right) \left(\lambda - \frac{\hbar \gamma B_0}{2}\right) &= 0 \\ \lambda^2 - \left( \frac{\hbar \gamma B_0}{2} \right)^2 &= 0 \\ \lambda^2 &= \left( \frac{\hbar \gamma B_0}{2} \right)^2 \\ \lambda &= \pm \frac{\hbar \gamma B_0}{2} \end{aligned} ]

Inserting into the Eigenvalue equation above yields the two Eigenvectors:

[ \begin{aligned} E_1 = -\frac{\hbar \gamma B_0}{2} &\to \vec{v}_1 = \begin{pmatrix} 1 \\ 0 \end{pmatrix} = \mid \uparrow > \\ E_2 = -\frac{\hbar \gamma B_0}{2} &\to \vec{v}_2 = \begin{pmatrix} 0 \\ 1 \end{pmatrix} = \mid \downarrow > \\ \end{aligned} ]

Thus we have seen the Eigenvalues and Eigenstates of the Spin operator $\hat{S}_z$ are:

[ \begin{aligned} \hat{H_0} \mid \uparrow > &= - \frac{\hbar \gamma B_0}{2} \mid \uparrow > \\ \hat{H_0} \mid \downarrow > &= \frac{\hbar \gamma B_0}{2} \mid \downarrow > \\ \end{aligned} ]

Density matrices and their time evolution

Since we are going to use density matrices to model our system later on let’s quickly recall the formalism used. A density matrix describes all states of a system and their population probability - as well as the overlap between non orthogonal states in the off diagonal terms.

[ \hat{\rho}(t) = \sum_{i} p_i \mid \psi_i(t) > < \psi_i(t) \mid ]

For a single pure state the density matrix thus would be very simple:

[ \hat{\rho}_{pure} = \mid \psi(t) > < \psi(t) \mid ]

The time evolution of a density matrix is described by the Liouville-von Neumann equation. This equation can be derived from the time dependent Schrödinger equation. First we reorder the Schrödinger equation:

[ \begin{aligned} i \hbar \frac{\partial}{\partial t} \mid \psi(t) > &= \hat{H}_0 \mid \psi(t) > \\ \to \frac{\partial \mid \psi(t) >}{\partial t} &= -\frac{i}{\hbar} \hat{H}_0 \mid \psi(t) > \end{aligned} ]

Then we insert the definition of the density matrix into the time derivative of the density matrix:

[ \begin{aligned} \frac{\partial \hat{\rho}}{\partial t} &= \frac{\partial}{\partial t} \mid \psi(t) > < \psi(t) \mid \\ &= \frac{\partial \mid \psi(t) >}{\partial t} < \psi(t) \mid + \mid \psi(t) > \frac{\partial < \psi(t) \mid}{\partial t} \\ &= -\frac{i}{\hbar} \hat{H}_0 \mid \psi(t) > < \psi(t) \mid + \frac{i}{\hbar} \mid \psi(t) > < \psi(t) \mid \hat{H}_0 \\ \frac{\partial \hat{\rho}}{\partial t} &= -\frac{i}{\hbar} \lbrack \hat{H}_0, \hat{\rho} \rbrack \end{aligned} ]

As one can see the dynamics of the density matrix is fully described by the commutator between the Hamiltonian and the density matrix - as usual the quantum mechanical description arises from the fact that the operators do not commute.

Initial state of our system

In thermal equilibrium the initial density matrix can be constructed by taking into account for the population probability of a state that is described by the Boltzmann distribution:

[ \hat{\rho}_{eq} = \frac{1}{ Tr(e^{-\frac{\hat{H}_0}{k_B T}}) } e^{-\frac{\hat{H}_0}{k_B T}} ]

The trace in the denominator provides the normalization of the probabilities of all states. Recall the trace is the sum of all diagonal elements:

[ \begin{aligned} Tr(\hat{A}) &= \sum_i < \psi_i \mid \hat{A} \mid \psi_i > \\ \to Tr(\hat{\rho}_{eq}) &= 1 \end{aligned} ]

Adding a driving field

In contrast to our treatment of the phenomenological Bloch equations we are adding the drive term before taking a look at the solution of the equation. The driving field is modeled in a classical way (i.e. we’re not modeling the exchange of photons but just assume the presence of the field that is coupling to the spin system):

[ \vec{B}_1(t) = 2 B_1 \cos(\omega_{RF} t) \hat{x} ]

We’ve chosen to orient the drive field in $\hat{x}$ direction - this is an arbitrary choice in the $xy$ plane that does not limit validity of our model (i.e. we can just rotate each other system around the $z$ axis so our model fits any driving field in the $xy$ plane - assuming dipole radiation).

The full Hamiltonian for our spin system now is

[ \hat{H}(t) = -\gamma B_0 \hat{S}_z - \gamma B_1 \cos(\omega_{RF} t) \hat{S}_x ]

Moving to the rotating frame, rotating frame approximation

As in the classical description solving the actual time dependent equations gets simpler when we move into a co-rotating frame. To model our rotation we’re utilizing a Unitarian transformation (recall that unitarity preserves length of state vectors - so normalization and also total energy is conserved).

[ \hat{U}(t) = e^{i \omega t \hat{S}_z} ]

This yields the Hamiltonian in the co-rotating field:

[ \begin{aligned} \hat{H}_{rot} &= \hat{U}^\dagger(t) \hat{H}(t) \hat{U}(t) - i \hbar \hat{U}^\dagger(t) \frac{\partial \hat{U}(t)}{\partial t} \\ &= e^{-i \omega t \hat{S}_z} \left( -\gamma B_0 \hat{S}_z - \gamma B_1 \cos(\omega_{RF} t) \hat{S}_x \right) e^{i \omega t \hat{S}_z} - i \hbar \hat{U}^\dagger(t) \frac{\partial \hat{U}(t)}{\partial t} \\ &= \hat{U}^\dagger(t) (-\gamma B_0 \hat{S}_z) \hat{U}(t) - \hat{U}^\dagger(t) (-\gamma B_1 \cos(\omega_{RF} t) \hat{S}_x) \hat{U}(t) - i \hbar \hat{U}^\dagger(t) \left(\frac{i \omega \hat{S}_z}{\hbar} \right) \hat{U}(t) \\ &= -\gamma B_0 \hat{S}_z - \gamma B_1 \cos(\omega_{RF} t) \left( \cos(\omega t) \hat{S}_x + \sin(\omega t) \hat{S}_y \right) - \hbar \omega \hat{S}_z \\ &= -\gamma B_0 \hat{S}_z - \gamma B_1 \cos(\omega_{RF} t) \cos(\omega t) \hat{S}_x - \gamma B_1 \cos(\omega_{RF} t) \sin(\omega t) \hat{S}_y \\ &= -\gamma B_0 \hat{S}_z - \gamma B_1 \left(\frac{1}{2} \cos((\omega_{RF} + \omega) t) + \frac{1}{2} \cos((\omega_{RF} - \omega)t) \right) \hat{S}_x - \gamma B_1 \left(\frac{1}{2} \sin((\omega + \omega_{RF}) t) + \frac{1}{2} \sin((\omega - \omega_{RF})t) \right) \hat{S}_y \end{aligned} ]

Now we perform the rotating wave approximation. We assume that terms $\omega_{RF} + \omega$ are oscillating very fast compared to our systems dynamic so they are assumed to average out over time. The detuning is now labeled $\Delta \omega = \omega - \omega_{RF}$ so we can simplify our Hamiltonian in the rotating frame to:

[ \hat{H}_{rot} = -\gamma B_0 \hat{S}_z - \gamma B_1 \frac{1}{2} \cos(\Delta \omega t) \hat{S}_x - \gamma B_1 \frac{1}{2} \sin(\Delta \omega t) \hat{S}_y ]

On resonance solutions and Rabi oscillations

To get a first impression of the equation and it’s implication we assume driving the system on resonance, i.e. $\omega_{RF} = \omega$ and thus $\Delta \omega = 0$:

[ \hat{H}_{rot}(t) = -\gamma B_0 \hat{S}_z - \frac{1}{2} \gamma B_1 \hat{S}_x ]

Now we can take a look at the time dependent Schrödinger equation:

[ \begin{aligned} i \hbar \frac{\partial}{\partial t} \mid \psi_{rot}(t) > &= \hat{H}_{rot}(t) \mid \psi_{rot}(t) > \\ &= \left(- \frac{\hbar \Delta \omega}{2} \hat{\sigma}_z - \frac{\hbar \gamma B_1}{2} \hat{\sigma}_x \right) \mid \psi_{rot}(t) > \end{aligned} ]

Lets assume our state can only be composed of a mixture of up and down states:

[ \mid \psi_{rot} (t) > = a(t) \mid \uparrow > + b(t) \mid \downarrow > ]

We can insert this into our differential equation:

[ \begin{aligned} i \hbar \frac{\partial}{\partial t} \left( a(t) \mid \uparrow > + b(t) \mid \downarrow > \right) &= \left(-\frac{\hbar \Delta \omega}{2} \hat{\sigma}_z - \frac{\hbar \gamma B_1}{2} \hat{\sigma}_x \right) \left( a(t) \mid \uparrow > + b(t) \mid \downarrow > \right) \\ &= \left(- \frac{\hbar \Delta \omega}{2} \begin{pmatrix} 1 & 0 \\ 0 & -1 \end{pmatrix} - \frac{\hbar \gamma B_1}{2} \begin{pmatrix} 0 & 1 \\ 1 & 0 \end{pmatrix} \right) \left( a(t) \mid \uparrow > + b(t) \mid \downarrow > \right) \\ &= \begin{pmatrix} -\frac{\hbar \Delta \omega}{2} & -\frac{\hbar \gamma B_1}{2} \\ -\frac{\hbar \gamma B_1}{2} & \frac{\hbar \Delta \omega}{2} \end{pmatrix} \left( a(t) \mid \uparrow > + b(t) \mid \downarrow > \right) \end{aligned} ]

Writing everything in vector form:

[ \begin{aligned} i \hbar \begin{pmatrix} a(t) \\ b(t) \end{pmatrix} &= \begin{pmatrix} - \frac{\hbar\Delta\omega}{2} & -\frac{\hbar\gamma B_1}{2} \\ -\frac{\hbar \gamma B_1}{2} & \frac{\hbar \Delta \omega}{2} \end{pmatrix} \begin{pmatrix} a(t) \\ b(t) \end{pmatrix} \\ &= \begin{pmatrix} -\frac{\hbar \Delta \omega}{2} a(t) - \frac{\hbar \gamma B_1}{2} b(t) \\ -\frac{\hbar \gamma B_1}{2} a(t) + \frac{\hbar \Delta \omega}{2} b(t) \end{pmatrix} \end{aligned} ]

This yields two coupled differential equations:

[ \begin{aligned} i \hbar \frac{\partial}{\partial t} a(t) &= -\frac{\hbar \Delta \omega}{2} a(t) - \frac{\hbar \gamma B_1}{2} b(t) \\ i \hbar \frac{\partial}{\partial t} b(t) &= -\frac{\hbar \gamma B_1}{2} a(t) + \frac{\hbar \Delta \omega}{2} b(t) \end{aligned} ]

These can be solved by first diving through the prefactor

[ \begin{aligned} \frac{\partial}{\partial t} a(t) &= i \frac{\Delta \omega}{2} a(t) + i \frac{\gamma B_1}{2} b(t) \\ \frac{\partial}{\partial t} b(t) &= i \frac{\gamma B_1}{2} a(t) - i \frac{\Delta \omega}{2} b(t) \end{aligned} ]

Now one can take the second derivative of the first equation and insert the second:

[ \begin{aligned} \frac{\partial^2}{\partial t^2} a(t) &= i \frac{\Delta \omega}{2} \frac{\partial a(t)}{\partial t} + i \frac{\gamma B_1}{2} \frac{\partial b(t)}{\partial t} \\ &= i \frac{\Delta\omega}{2} + i \frac{\gamma B_1}{2}\left(i \frac{\gamma B_1}{2} a(t) - i \frac{\Delta \gamma}{2} b(t) \right) \\ &= i \frac{\Delta \omega}{2} \left(i \frac{\Delta \omega}{2} a(t) + i \frac{\gamma B_1}{2} b(t) \right) + i \frac{\gamma B_1}{2}\left(i \frac{\gamma B_1}{2} a(t) - i \frac{\Delta \gamma}{2} b(t) \right) \\ &= \left(-\frac{\Delta \omega^2}{4} - \frac{\gamma^2 B_1^2}{4} \right) a(t) - \underbrace{\left( \frac{\Delta \omega \gamma B_1}{2} + \frac{\Delta \omega \gamma B_1}{4} \right)}_{0} b(t) \\ &= \left(-\frac{\Delta \omega^2}{4} - \frac{\gamma^2 B_1^2}{4} \right) a(t) \end{aligned} ]

This now simplifies to a very simple second order differential equation:

[ \begin{aligned} \frac{\partial^2}{\partial t^2} a(t) + \underbrace{\left(\frac{\Delta \omega^2}{4} + \frac{\gamma^2 B_1^2}{4} \right)}_{\frac{\Omega^2}{4}} a(t) &= 0 \\ \frac{\partial^2}{\partial t^2} a(t) + \frac{\Omega^2}{4} a(t) &= 0 \end{aligned} ]

Here we have introduced the term $\Omega = \sqrt{\Delta \omega^2 + \gamma^2 B_1^2}$ that will turn out to be the Rabi frequency very soon

Using the Ansatz $a(t) = c_1 \sin(c_2 t) + c_3 \cos(c_2 t)$ we arrive at

[ \begin{aligned} a(t) &= c_1 \sin(c_2 t) + c_3 \cos(c_2 t) \\ \partial_t a(t) &= c_2 * c_1 \cos(c_2 t) - c_3 c_2 \sin(c_2 t) \\ \partial_t^2 a(t) &= -c_1 * c_2^2 \sin(c_2 t) - c_3 c_2^2 \cos(c_2 t) \\ &= -c_2^2 a(t) \end{aligned} ] [ \begin{aligned} \to c_2^2 a(t) + \frac{\Omega^2}{4} a(t) &= 0 \\ \to c_2^2 &= \frac{\Omega^2}{4} \\ \to c_2 &= \pm \frac{\Omega}{2} \\ \to a(t) &= c_1 \sin\left(\frac{\Omega}{2} t\right) + c_3 \cos\left(\frac{\Omega}{2} t\right) \end{aligned} ]

Calculating $b(t)$ in an analogue way results in the similar expression

[ b(t) = c_4 \sin\left(\frac{\Omega}{2} t\right) + c_5 \cos\left(\frac{\Omega}{2} t\right) ]

The prefactors $c_1, c_3, c_4$ and $c_5$ have to be determined by initial conditions. As we can see the states are oscillating between up and down state.

The probability of locating a system in states $\mid \uparrow >$ and $\mid \downarrow >$ can be determined by taking the absolute values:

[ \begin{aligned} P(\mid \uparrow >) &= \mid a(t) \mid^2 \\ P(\mid \downarrow >) &= \mid b(t) \mid^2 \end{aligned} ]

Let’s assume we start at time $t=0$ fully in state $\mid \uparrow >$. Then we get the following constants:

[ \begin{aligned} \mid \psi(t=0) > &= 1 * \mid \uparrow > + 0 * \mid \downarrow > \\ \to a(t = 0) = 1 &\to a(t) = \cos\left(\frac{\Omega}{2}t\right) \\ \to b(t = 0) = 0 &\to b(t) = \sin\left(\frac{\Omega}{2}t\right) \end{aligned} ]

As for the electric dipole this shows oscillations between excited and ground state due to coherent exchange of energy between the driving field $B_1$ and the spin system.

The oscillation between ground and excited state and back is called Rabi oscillation, the frequency $\Omega$ is the Rabi frequency

Again one can define $\frac{\pi}{2}$ and $\pi$ pulses which are the extrema - for a $\pi$ pulse the system is transitioned fully from ground to excited state, for a $\frac{\pi}{2}$ pulse the system is kicked into the $xy$ plane.

Introducing relaxation processes

Introducing relaxation processes is a little bit harder. Up until now we modeled a unitary system (which means energy is conserved). To model the dynamic using the density matrix formalism one uses the commutator with the Hamiltonian (also called the von Neumann equation):

[ \begin{aligned} \frac{\partial \hat{\rho}}{\partial t} = -\frac{i}{\hbar} \lbrack \hat{H}_0, \hat{\rho} \rbrack \end{aligned} ]

This equation works well for closed systems (due to it’s unitarity) - when coupling the system to a bath or reservoir one needs to introduce some non Unitarian terms. This is usually done in the context of the Lindblad master equation (see the Gorini-Kossakowski-Sudarshan and Linblad theorem)

[ \begin{aligned} \frac{\partial \hat{\rho}(t)}{\partial t} = -\frac{i}{\hbar} \lbrack \hat{H}, \hat{\rho}(t) \rbrack + \sum_i \left( L_i \rho(t) L_i^\dagger - \frac{1}{2} \lbrace L_i^\dagger L_i, \hat{\rho}(t) \rbrace \right) \end{aligned} ]

The Lindblad operators $L_j$ describe the interaction with the environment.

Longitudinal relaxation T1

To introduce longitudinal relaxation one has to account for the relaxation of the $z$ component back into equilibrium state by dissipation of energy to the environment (the bath). This is done by decaying from the excited state to the ground state with a decay rate $\frac{1}{T_1}$. This can be modeled with the ladder operator $\hat{\sigma_{-}}$:

[ \begin{aligned} \hat{\sigma}_{-} &= \begin{pmatrix} 0 & 0 \\ 1 & 0 \end{pmatrix} \\ \hat{L}_1 &= \sqrt{\frac{1}{T_1}} \hat{\sigma}_{-} \end{aligned} ]

Transverse relaxation T2

The introduction of transverse relaxation works similar - this is the loss of coherence in the $xy$ plane (i.e. the decay of off diagonal elements in the density matrix - dephasing).

[ \begin{aligned} \hat{L}_2 &= \sqrt{\frac{1}{2 T_2}} \hat{\sigma}_z \end{aligned} ]

This is dephasing without an change in total energy.

Complete equation including relaxation

Inserting into the master equation yields

[ \frac{\partial \hat{\rho}}{\partial t} = -\frac{i}{\hbar} \lbrack \hat{H}, \hat{\rho} \rbrack + \frac{1}{T_1} \left(\hat{\sigma}_{-} \hat{\rho} \hat{\sigma}_{+} - \frac{1}{2} \lbrace \hat{\sigma}_{+} \hat{\sigma}_{-}, \hat{\rho} \rbrace \right) + \frac{1}{2 T_2} \left( \hat{\sigma}_z \hat{\rho} \hat{\sigma}_z - \hat{\rho} \right) ]

Expectation values of Magnetization vectors

To determine the magnetization of the whole system - this is the quantity we have also described phenomenological before - we apply our spin operators and trace out the density matrix:

[ \begin{aligned} M_x(t) &= Tr\left(\hat{\rho}(t) \hat{S}_x\right) = Tr\left(\hat{\rho}(t) \frac{\hbar}{2} \hat{\sigma}_x\right) \\ M_y(t) &= Tr\left(\hat{\rho}(t) \hat{S}_y\right) = Tr\left(\hat{\rho}(t) \frac{\hbar}{2} \hat{\sigma}_y\right) \\ M_z(t) &= Tr\left(\hat{\rho}(t) \hat{S}_z\right) = Tr\left(\hat{\rho}(t) \frac{\hbar}{2} \hat{\sigma}_z\right) \\ \end{aligned} ]

Now let’s take a look at the time evolution of those operators. Recall that

[ \begin{aligned} \frac{\partial}{\partial t} \langle \hat{O} \rangle &= \frac{\partial}{\partial t} Tr\left(\hat{\rho}(t) \hat{O} \right) \\ &= Tr \left(\frac{\partial \hat{\rho}}{\partial t} \hat{O} \right) \end{aligned} ]

This yields the time derivative for arbitrary operators for the Lindblad equation:

[ \frac{\partial}{\partial t} \langle \hat{O} \rangle = Tr \left( \left( -\frac{i}{\hbar} \lbrack \hat{H}_{rot}, \hat{\rho} \rbrack + \sum_k \left(\hat{L}_k \hat{rho} \hat{L}_k^\dagger - \frac{1}{2} \lbrace \hat{L}_k^\dagger \hat{L}_k, \hat{\rho} \right) \right) \hat{O} \right) ]

As one can see equation is consisting - as expected - of two parts:

• A Unitarian part covering the system without relaxation
• The non Unitarian term described by the Lindblad operators describing the relaxation process

Unitary part (Hamiltonian)

[ \begin{aligned} Tr \left(-\frac{i}{\hbar} \lbrack \hat{H}_{rot}, \hat{\rho} \rbrack \hat{O} \right) &= -\frac{i}{\hbar} Tr \left( \lbrack \hat{H}_{rot}, \hat{\rho} \rbrack \hat{O} \right) \\ \to \frac{\partial}{\partial t} \langle \hat{O} \rangle &= \frac{i}{\hbar} \langle \lbrack \hat{H}_{rot}, \hat{O} \rbrack \rangle \end{aligned} ]

We can now apply this to our spin operators:

For the x component:

[ \begin{aligned} M_x(t) &= Tr \left( \hat{rho}(t) \hat{S}_x \right) \\ \lbrack \hat{H}_{rot}, \hat{S}_x \rbrack &= -\frac{\hbar \Delta \omega}{2} \lbrack \hat{\sigma}_z, \hat{\sigma}_x \rbrack \\ &= i \hbar \frac{\Delta \omega}{2} \hat{\sigma}_y \\ \to \frac{\partial M_x(t)}{\partial t} &= -\Delta \omega M_y(t) \end{aligned} ]

For the y component:

\begin{aligned} M_y(t) &= Tr \left( \hat{\rho}(t) \hat{S}_y \right) \\ \lbrack \hat{H}_{rot}, \hat{S}_y \rbrack &= - \frac{\hbar \gamma B_1}{2} \lbrack \hat{\sigma}_x, \hat{\sigma}_y \rbrack \\ &= i \hbar \frac{\gamma B_1}{2} \hat{\sigma_z} \\ \to \frac{\partial M_y(t)}{\partial t} &= \Delta \omega M_z(t) - \gamma B_1 M_z(t) \end{aligned}

And for the z component:

\begin{aligned} M_z(t) &= Tr \left( \hat{\rho}(t) \hat{S}_z \right) \\ \lbrack \hat{H}_{rot}, \hat{S}_z \rbrack &= - \frac{\hbar \gamma B_1}{2} \lbrack \hat{\sigma}_x, \hat{\sigma}_z \rbrack \\ &= i \hbar \frac{\gamma B_1}{2} \hat{\sigma}_y \\ \to \frac{\partial M_z(t)}{\partial t} &= \gamma B_1 M_y(t) \end{aligned}

This is consistent with the phenomenological model in the rotating frame (note that it has been assumed to drive on resonant - this is the reason why we only see the driving field acting onto the $y$ components derivative).

Relaxation (Lindblad terms)

Longitudinal relaxation (T1)

Recall that our operator is given by

[ \begin{aligned} \hat{L}_1 &= \sqrt{\frac{1}{T_1}} \hat{\sigma}_{-} \\ &= \sqrt{\frac{1}{T_1}} \begin{pmatrix} 0 & 0 \\ 1 & 0 \end{pmatrix} \\ \to \frac{\partial \hat{\rho}}{\partial t} \mid_{T_1} &= \frac{1}{T_1} \left( \hat{\sigma}_{-} \hat{\rho} \hat{\sigma}_{+} - \frac{1}{2} \lbrace \hat{\sigma}_{+} \hat{\sigma}_{-}, \hat{\rho} \rbrace \right) \end{aligned} ]

As we can see this influences the $z$ component via the ladder operators $\hat{\sigma}_{+}$ and $\hat{\sigma}_{-}$ which indirectly influences of course the absorptive and dispersive components $M_x$ and $M_y$ via the unitary interactions.

[ \begin{aligned} \frac{\partial M_z(t)}{\partial t} \mid_{T_1} &= - \frac{M_z(t) - M_0}{T_1} \\ \frac{\partial M_x(t)}{\partial t} \mid_{T_1} &= 0 \\ \frac{\partial M_y(t)}{\partial t} \mid_{T_1} &= 0 \\ \end{aligned} ]

Transverse relaxation (T2)

Recall that the operator for transverse relaxation is given by

[ \begin{aligned} \hat{L}_2 &= \sqrt{\frac{1}{2 T_2}} \hat{\sigma}_z \\ &= \sqrt{\frac{1}{2 T_2}} \begin{pmatrix} 1 & 0 \\ 0 & -1 \end{pmatrix} \\ \to \frac{\partial \hat{\rho}}{\partial t} \mid_{T_2} &= \frac{1}{2 T_2} \left( \hat{\sigma}_z \hat{\rho} \hat{\sigma}_z - \hat{\rho} \right) \end{aligned} ]

As one can see this does not affect the $z$ component directly - but the $x$ and $y$ components:

[ \begin{aligned} \frac{\partial M_x(t)}{\partial t} \mid_{T_2} &= -\frac{M_x(t)}{T_2} \\ \frac{\partial M_y(t)}{\partial t} \mid_{T_2} &= -\frac{M_y(t)}{T_2} \end{aligned} ]

Unitary and non unitary contributions

Assembling the equations for all three axis yields

[ \begin{aligned} \frac{\partial M_x(t)}{\partial t} &= -\Delta \omega M_y(t) - \frac{1}{T_2} M_x(t) \\ \frac{\partial M_y(t)}{\partial t} &= \Delta \omega M_x(t) - \gamma B_1 M_z(t) - \frac{1}{T_2} M_y(t) \\ \frac{\partial M_z(t)}{\partial t} &= \gamma B_1 M_y(t) - \frac{M_z(t) - M_0}{T_1} \end{aligned} ]

As one can see we have just recovered the Bloch equations including relaxation terms for an on resonant system from a microscopic theory of nuclear magnetic resonance (NMR) or electron paramagnetic resonance (EPR).

Linewidth and excited state lifetime

Another interesting property is the expected (minimum) width of the signal as well as the expected minimum lifetime of the excited state. For an spin system that is only affected by the longitudinal and transversal relaxation processes this is mainly limited by the relaxation time $T_2$ (i.e. transversal relaxation). Keep in mind that for practical systems effects like magnetic field inhomogenity can quickly yield way larger contributions.

Limits by transversal relaxation

Taking a look at the spectral function of the transversal magnetization

[ M_{xy}(t) \propto M_{xy}(t) e^{-\frac{t}{T_2}} ] [ \begin{aligned} \mathfrak{F}\lbrace M_{xy}(t) \rbrace &= \int_{-\infty}^{\infty} M_{xy}(t) e^{-i 2 \pi \nu t} dt \\ &= \int_{-\infty}^{\infty} M_{xy}(0) e^{-\frac{t}{T_2}} e^{-i2\pi\nu t} dt \\ &= \int_{-\infty}^{\infty} M_{xy}(0) e^{-\left(\frac{1}{T_2} + i 2 \pi \nu \right)t} dt \end{aligned} ]

Using

[ \int_{0}^{\infty} e^{-\alpha t} = \frac{1}{\alpha} ]

we arrive at

[ \begin{aligned} \mathfrak{F}\lbrace M_{xy}(t) \rbrace &= M_{xy}(0) \frac{1}{\frac{1}{T_2} + i 2 \pi \nu} \\ &= M_{xy}(0) \frac{1}{\frac{1 + i \pi \nu T_2}{T_2}} \\ &= \frac{T_2 M_{xy}(0)}{1 + i 2 \pi \nu T_2} \\ &= T_2 M_{xy}(0) \lbrace \frac{1 - i 2 \pi \nu T_2}{1 + (2 \pi \nu T_2)^2} \rbrace \\ &= \lbrace \frac{T_2 M_{xy}(0)}{1 + (2 \pi \nu T_2)^2} - i \frac{2 \pi \nu T_2^2 M_{xy}(0)}{1 + (2 \pi \nu T_2)^2} \rbrace \\ \to \mid \mathfrak{F}\lbrace M_{xy}(t) \rbrace \mid &= \sqrt{ \left( \frac{T_2 M_{xy}(0)}{1 + (2 \pi \nu T_2)^2} \right)^2 + \left( \frac{2 \pi \nu T_2^2 M_{xy}(0)}{1 + (2 \pi \nu T_2)^2} \right) } \\ &= \sqrt{\frac{ (T_2 M_{xy}(0))^2 (1 + (2 \pi \nu T_2)^2) }{(1 + (2 \pi \nu T_2)^2)^2}} \\ &= T_2 M_{xy}(0) \frac{\sqrt{1 + (2 \pi \nu T_2)^2}}{1 + (2 \pi \nu T_2)^2} \\ &= T_2 M_{xy}(0) \frac{1}{\sqrt{1 + (2 \pi \nu T_2)^2}} \\ &= \frac{T_2 M_{xy}(0)}{\sqrt{(2 \pi T_2)^2 \nu^2 + 1}} \\ &= \frac{\frac{1}{T_2} T_2 M_{xy}(0)}{\frac{1}{T_2} \sqrt{(2 \pi T_2)^2 \nu^2 + 1}} \\ &= \frac{M_{xy}(0)}{\sqrt{(2 \pi \nu)^2 + \left(\frac{1}{T_2}\right)^2}} \end{aligned} ]

The last step has been done to be able to compare this function to a typical Lorentzian distribution (Breit-Wigner function):

[ L(\nu_2) = \frac{A}{(\nu_2 - \nu_0)^2 + \Gamma^2} ]

Comparing factors yields:

• $\nu_0 = 0$
• $\Gamma = \frac{1}{T_2}$
• $\nu_2 = 2 \pi \nu$

Taking a look at the full width at half maximum of the Lorentzian we can determine the natural line width limit by transversal relaxation processes:

[ \begin{aligned} \mathfrak{L}(\nu_{FWHM}) &= \frac{A}{2} \\ &= \frac{A}{\nu_{FWHM}^2 + \Gamma^2} \\ \to A(\nu_{FWHM}^2 + \Gamma^2) &= 2A \\ \nu_{FWHM} &= \pm \Gamma \\ \to \nu_{FWHM} &= \pm \frac{1}{T_2} \\ \to \Delta \nu &= 2 \nu_{FWHM} \\ &= \frac{2}{T_2} \\ &= \frac{1}{\pi T_2} \end{aligned} ]

In the last line we moved to a line width in angular frequency instead of Hz.

Associated lifetime limit by transversal relaxation

The lower lifetime boundary associated with the relaxation time $T_2$ can be derived from the energy-time uncertainty relation (and recalling that $E = h \nu$ and $\hbar = \frac{h}{2 \pi}$):

[ \begin{aligned} \Delta E \Delta t &\geq \frac{\hbar}{2} \\ h \Delta \nu \Delta t &\geq \frac{\hbar}{2} \\ \Delta \nu \Delta t &\geq \frac{1}{4 \pi} \\ \frac{1}{\pi T_2} \tau &\geq \frac{1}{4 \pi} \\ \tau &\geq \frac{T_2}{4} \end{aligned} ]

This article is tagged: Physics, Electrodynamics, Basics, How stuff works, Tutorial, Math, Measurements, Quantum mechanics, Quantum optics