NMR

Introduction

In NMR experiments, transitions between different states of a quantum mechanical system are induced via a periodic field, usually much weaker than typical fields of the system itself. Its frequency, however, must be comparable to the energy difference between the states (ℏω ≈ ΔE) to have any effect. If ℏω = ΔE, all particles can be transferred from their original state to another and back again in a periodic cycle. Its frequency is the so called Rabi frequency.
This notebook considers a two-state system, namely a Spin - (1/2)particle in a constant magnetic field. The time dependent Schrödinger equation for this system with an additional periodic magnetic field is solved numerically. The Rabi frequency and the minimal population of the original state are then determined for different strengths and frequencies of the perturbation.

Clear["Global`*"];   
Needs["Statistics`NonlinearFit`"]; 
    Needs["Graphics`MultipleListPlot`"];

Schrödinger equation

The Hamiltonian of a Spin-1/2 particle in a homogenous, constant magnetic field in z-direction is:

(Subscript[H, 0]) =  ((e g)/(2 m)) (Subscript[B, z]) (Subscript[S, z])	,	(Subscript[S, z]) = (ℏ/2) (Subscript[σ, 3])

with (Subscript[ω, 0])=((e g (Subscript[B, z]))/(2  m)), this becomes:

(Subscript[H, 0]) =  (Subscript[ω, 0]) (Subscript[S, z])

The eigenvalues of this Hamiltonian are (Subscript[E, σ])=(σ/2) ℏ(Subscript[ω, 0]), with σ=±1. The probability of finding the system in one of these eigenvectors is time-independent:

|⟨σ | ψ(
                    t)⟩\( | \^2\)=|⟨σ | 
                (e^(((-it)/ℏ) H)) | ψ(
                  t = 0)⟩\( | \^2\)= |(e^((it/ℏ) (Subscript[E, σ])))⟨σ | ψ(
                  0)⟩\( | \^2\)= \(TraditionalForm\`\(\(|\)\(⟨σ | ψ(
              0)⟩\)\( | \^2\)\)\)

In order to obtain a transition between the eigenstates, we need a time-dependent Hamiltonian. One possibility is an additional periodic field in x-direction:

H = (Subscript[H, 0]) + ((e g)/(2 m)) (Subscript[B, x]) cos(ωt) (Subscript[S, x])=(Subscript[ω, 0]) (Subscript[S, z])+(Subscript[ω, 1]) cos(ωt) (Subscript[S, x])

We express H and |ψ(t)⟩ in  the (Subscript[S, z]) - basis, the eigenbasis of (Subscript[H, 0]), the eigenvectors are |+z⟩=(({{1}, {0}})) for σ=+1 and |-z⟩=(({{0}, {1}})) for σ=-1:

(Subscript[σ, 1]) = (({{0, 1}, {1, 0}}));  (Subscript[σ, 2]) = (({{0, -I}, {I, 0}}));  (Subscript[σ, 3]) = (({{1, 0}, {0, -1}}));

(Subscript[H, 0]) = (Subscript[ω, 0]) (Subscript[σ, 3])/2;
H = (Subscript[H, 0]) + ((Subscript[ω, 1])/2) Cos[ω t] (Subscript[σ, 1]); 
    H // MatrixForm

With 	\(\(\(\(|\)\(ψ(
                          t)\)\)\(⟩\)\) = ∑\+σ\)(Subscript[C, σ])(t)|σ⟩=(({{a(t)}, {b(t)}})),

the time-dependent Schroedinger equation becomes:

i (∂/∂)t(({{a(t)}, {b(t)}}))=H (({{a(t)}, {b(t)}}))=(1/2)(({{(Subscript[ω, 0]), (Subscript[ω, 1]) cos(ωt)}, {(Subscript[ω, 1]) cos(ωt), -(Subscript[ω, 0])}}))(({{a(t)}, {b(t)}}))

Numerical Solution

We solve this equation for the coefficients a(t) and b(t) numerically.

numsol[ω0_, ω1_,  ω_] := 
        Block[{sol, t}, 
ωest = 
            Max[(ω - ω0), ω1/
                2]; 
tmax = 
            N[2  π/ωest]; 
sol = 
            NDSolve[{I a'[
                      t] == ω0 a[t]/
                        2 + ω1 Cos[ω t] b[t]/2, 
                I b'[
                      t] == -ω0 b[t]/
                        2  + ω1 Cos[ω t] a[t]/2, 
                a[0] == 1, b[0] == 0}, {a[t], b[t]}, {t, 0,  
                tmax}]; 
pup[t_] = 
            Abs[a[t]]^2 /. sol[[1]]; 
pdown[t_] = 
            Abs[b[t]]^2 /. sol[[1]];
];

numsol[ω0 = 1, ω1 = 0.1, ω = 
          0.95];
plOsc = 
        Plot[{pup[t], pdown[t]}, {t, 0,  tmax}, 
          PlotRange -> {{0, 1.1  tmax}, {0, 1}}, 
          PlotStyle -> {Hue[1], Hue[0.6]}, 
          AxesLabel -> {"t", ""}, 
          DisplayFunction -> Identity];
\(\(Show[plOsc, 
        Graphics[{Text[\*""|a(t)\!\(\( | \^2\)\)"", {1.1\ tmax, 
                pup[tmax]}, TextStyle -> {FontColor -> Hue[1]}], 
            Text[\*""|b(t)\!\(\( | \^2\)\)"", {1.1\ tmax, pdown[tmax]}, 
              TextStyle -> {FontColor -> Hue[0.6]}]}], 
        DisplayFunction -> $DisplayFunction];\)\)

The probabilities to find the state in a particular eigenstate are now time dependent, we obtain a superposition of a fast and a slow oscillation. The frequency of the slow oscillation is called the Rabi frequency.

Rabi Frequency

In order to find the involved frequencies, we perform a Fourier transform on \(TraditionalForm\`\(\(|\)\(a(t)\)\( | \^2\)\)\):

(*Fourier*) 
numfourier[fname_, tmax_] := 
        Block[{t, n00, n0m1, n0p1, nmax = 250, at, lmax, 
            i}, 
func[t_] = 
            fname[t]; 
meanpop = 
            NIntegrate[func[t], {t, 0, tmax}]/tmax; 
at = 
            Table[func[i (tmax/nmax)], {i, 0, nmax}] - 
              meanpop; 
fA = 
            Fourier[at]; 
lmax = 
            Position[Abs[Take[fA, Round[nmax/2]]], 
              Max[Abs[Take[fA, Round[nmax/2]]]]]; 
n00 = 
            Min[Flatten[
                lmax]]; 
 (*to slightly impltove the accuracy of the result:*) 
n0m1 = 
            Max[n00 - 1, 1]; 
n0p1 = 
            n00 + 1; 
n0 = (n00*Abs[fA[[n00]]] + 
                    n0m1*Abs[fA[[n0m1]]] + 
                    n0p1*Abs[fA[[n0p1]]])/(Abs[fA[[n00]]] + 
                    Abs[fA[[n0m1]]] + Abs[fA[[n0p1]]]) - 
              1;
];

numfourier[pup, tmax];
ListPlot[Abs[fA],  PlotJoined -> True, 
        PlotRange -> {{0, 100}, All}, PlotStyle -> {Hue[1]}, 
        AxesLabel -> {"ω", ""}];

We get a large peak at a low frequency (the Rabi frequency) and a small peak at a high frequency. The approximate Rabi frequency is:

\(\(Print[\*""\!\(ω\_R\)= "", ωest\ n0];\)\)

Unfortunately, the frequency resolution is not very good and we need an additional nonlinear fit. We will use the frequency obtained from the Fourier transform as a starring value for this fit:

(*Nonlinear Fit*) 
 (*The starting value for the frequency needs to be quite close!*) 
cosfit[funcname_, offset_, amplitude_, frequency_] := 
        Block[{nmax = 250, afit, i, t, amp1, a0, 
            n1}, 
func[t_] = 
            funcname[t]; 
afit = 
            Table[{N[(( i)/nmax)], 
                func[N[i (( tmax)/nmax)]]}, {i, 0, 
                nmax}]; 
fitres = 
            BestFitParameters /. 
              NonlinearRegress[afit, a0 + amp1 Cos[2 π n1 t], 
                t, {{a0, offset}, {amp1, amplitude}, {n1, frequency}}, 
                RegressionReport -> 
                  BestFitParameters];
];

cosfit[pup, meanpop, 1 - meanpop, n0];
ωR = ωest n1 /. fitres;
ampR = amp1 /. fitres;
meanpop = a0 /. fitres;
plfit = 
        Plot[{meanpop + ampR Cos[ωR t], pup[t]}, {t, 0, tmax}, 
          PlotStyle -> {{Hue[1], Thickness[0.01]}, {GrayLevel[0]}}, 
          DisplayFunction -> Identity];
\(\(Show[plfit, 
        Graphics[{Text["fit", {\ tmax, meanpop}, 
              TextStyle -> {FontColor -> Hue[1]}], 
            Text[\*""|a(t)\!\(\( | \^2\)\)"", {\ tmax, 
                meanpop - ampR/2}]}], AxesLabel -> {"t", ""}, 
        DisplayFunction -> $DisplayFunction];\)\)
Print["Mean Propability of measuring |+z⟩: ", 
        meanpop];
Print["Amplitude of oscillation: ", 
        ampR];
Print["Rabi frequency: ", ωR, "  Period of \
Rabi-oscillation: ", 2 π/ωR];

Parameter Dependence

We want to observe the effect of a variation in ω1 and ω, i.e. the strength and the frequency of the perturbation, on the Rabi frequency and the amplitude.

ω1list = {0.05, 0.075, 0.1, 0.15, 0.2, 0.25, 
        0.3}; ωlist = {0.8, 0.9, 0.95, 0.99, 1, 1.01, 1.05, 1.1, 1.2}; 
    colldata = 
      Table[Table[{0, 0, 0, 
            0 }, {Length[ω1list]}], {Length[ωlist]}];

This may take some time:

imax = Length[ωlist]; 
    jmax = Length[ω1list];
\(\(Do[\(Do[
numsol[
              1, ω1list[\([j]\)], ωlist[\([i]\)]]; 
numfourier[pup, tmax]; 
cosfit[pup, 
              meanpop, 1 - meanpop, 
              n0]; 
ωR = ωest\ n1 /. 
                fitres; 
ampR = 
              amp1 /. fitres; 
 (*\(Print["ω1= ", ω1list[\([j]\)], " \[RawEscape]ω= ", ωlist[\([i]\)], \*""    (Subscript[ω, R])= "", ωR, "  Amplitude: ", ampR];\)*) 
colldata[\([i, 
                  j]\)] = {ω1list[\([j]\)], ωlist[\([i]\)], ωR, ampR};
, {j, 1, 
              jmax}];\)
, {i, 1, imax}];\)\)

The Rabi frequency as a function of the frequency ω of the perturbation and of its strength (Subscript[ω, 1]):

\(\(plFreq = 
        Show[Graphics3D[{Line /@ 
                Transpose[colldata[\([All, All, {2, 1, 3}]\)]], 
              Line /@ colldata[\([All, All, {2, 1, 3}]\)]}\ , 
            Boxed -> True, Axes -> True, 
            AxesLabel -> {\*""ω1/(Subscript[ω, 0])"", \*""ω/(Subscript[ω, 0])"", \*""(Subscript[ω, R])/(Subscript[ω, 0])""}, PlotLabel -> "Rabi frequency"]];\)\)

Save plot to ps-file

SetDirectory["~/mathem"];
Export["rabifreq.eps", plFreq, "EPS"];

The amplitude of the Rabi oscillation:

\(\(plAmp = 
        Show[Graphics3D[{Line /@ 
                Transpose[colldata[\([All, All, {2, 1, 4}]\)]], 
              Line /@ colldata[\([All, All, {2, 1, 4}]\)]}\ , 
            Axes -> True, 
            AxesLabel -> {\*""ω/(Subscript[ω, 0])"", \*""ω1/(Subscript[ω, 0])"", ""}, 
            AxesEdge -> {{\(-1\), \(-1\)}, {1, \(-1\)}, {\(-1\), \(-1\)}}, PlotLabel -> "Amplitude"]];\)\)

We see, that the oscillation is stronger for ω close to (Subscript[ω, 0]). For ω = (Subscript[ω, 0]), all particles go from |+z⟩ at the beginning to |-z⟩ at π/(Subscript[ω, R]). Outside this resonance, some particles remain in the |+z⟩ state.

Save plot to ps-file

SetDirectory["~/mathem"];
Export["rabiamp.eps", plAmp, "EPS"];