NMR and RABI FREQUENCY

Maria Daghofer,  Wolfgang von der Linden
Institut fuer Theoretische Physik
Technische Universitaet Graz
A-8010 Graz,Austria
Projekt Multimediale Lehre des
Zentralen Informatikdienstes der TU Graz

2001-12-13

Parameters

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

Strength of the field in z-direction:

Strength of additional periodic field in x-direction:

Frequency of periodic field:

Schrödinger equation

(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

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]];
];

Clear[a, b];
numsol[ω00, ω1, ωf];
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 '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 (fast oscillation). 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. It is necessary to have a quite accurate starting value for the frequency, otherwise, we might obtain a false solution.

(*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["The plot compares the original numerical solution (black) to \
the cosine-fit (red)\n"];
Print[
        TableForm[{{"Mean Probability of measuring |+z⟩: \
", meanpop}, 
{"Amplitude of oscillation: ", 
              ampR}, 
{"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} ω00; ωlist = {0.8, 0.9, 0.95, 0.99, 1, 1.01, 
          1.05, 1.1, 1.2} ω00;
colldata = 
        Table[Table[{0, 0, 0, 
              0 }, {Length[ω1list]}], {Length[ωlist]}];

This may take some time:

Clear[a, b];
imax = Length[ωlist]; 
    jmax = Length[ω1list];
\(\(Do[\(Do[
numsol[ω00, ω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 -> {\*""ω/(Subscript[ω, 0])"", \*""ω1/(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"];

Analytic Approximation

Clear[ω, ωR, a0, b0, a, b, 
        t];

Im[(Subscript[ω, 1])]  ^= 0; 
    Im[Δω]  ^= 0;
a[t_] := a0[t] Exp[-I (Subscript[ω, 0]) t/2]; 
    b[t_] := b0[t] Exp[I (Subscript[ω, 0]) t/2];
fulleq = {D[a[t], t], D[b[t], t]} == -I H .  {a[t], b[t]} // 
            FullSimplify // TrigToExp;
Solve[fulleq, {D[a0[t], t], D[b0[t], t]}][[1]] // 
      ColumnForm

diffeqs = {D[a0[t], 
            t] == -I (Subscript[ω, 1])/
              4 Exp[-I Δω t] b0[t], 
        D[b0[t], 
            t] == -I (Subscript[ω, 1 ])/4 Exp[
              I Δω t] a0[t]}

Mathematica was not able to directly solve this system of equations, we therefore express b0 in terms of a0' and arrive at a second order differential equation:

diffeqb = Solve[diffeqs[[1]], 
          b0[t]][[1]]
diffeqa = 
      FullSimplify[
        diffeqs[[2]] /. 
          D[diffeqb, t], {(Subscript[ω, 1]) ≠ 
            0, (E^(I t Δω)) ≠ 0}]

In this form, the equation is solved and we obtain for a0:

DSolve[diffeqa, a0[t], t] // ExpandAll

This solution contains a phase \(TraditionalForm\`e\^\(\(\(-i\)\ t\ Δω\)\/2\)\)and an oscillation, which can be expressed as a cosine with a phase angle Δ: \(TraditionalForm\`cos(\(t\/4\) \@\(4\ Δω\^2 - ω\_1\%2\)\  + \ Δ)\).  The frequency of the cosine is half the Rabi frequency and we set: \(TraditionalForm\`Ω = \(1\/2\) \@\(Δω\^2 - \((ω\_1\/2)\)\^2\)\) =((Subscript[ω, Rabi])/2).

Clear[a0, b0];
a0[t_] := 
        Evaluate[DSolve[diffeqa, a0[t], t][[1, 1, 
                  2]] /. (Sqrt[-4 Δ(ω^2) - (Subscript[ω, 1]^2)]) -> -4 I Ω // 
            FullSimplify];
Print["a0(t)= ", a0[t]];
b0[t_] := 
        Evaluate[
          diffeqb[[1, 2]] // FullSimplify];
Print["b0(t)= ", b0[t]];

We set the initial conditions a(0)=1 and b(0)=0, this determines the constants C[1] and C[2]:

consts = 
      Solve[{a[0] == 1, b[0] == 0}, {C[1], C[2]}]

The approximate solutions for a(t) and b(t) are thus:

Print["a(t)=  ", 
        a[t] /. consts[[1]] // Simplify];
Print["b(t)=  ", 
        b[t] /. consts[[1]] // Simplify];

We can now calculate the probabilities to find the system in |+z⟩ or in |-z⟩ at a given time:

Im[(Subscript[ω, R])]  ^= 0; 
    Im[t]  ^= 0; (Subscript[ω, R]) > 0;
a2abs = Abs[a[t]]^2 /. consts[[1]] // FunctionExpand // 
            ComplexExpand // Simplify;
\(\(Print[\*""Probab. for |+z>:  |a(t)\!\(\( | \^2\)\)= "", 
        a2abs];\)\)
b2abs = Abs[b[t]]^2 /. consts[[1]] // FunctionExpand // 
            ComplexExpand // Simplify;
\(\(Print[\*""Probab. for |-z>:  |b(t)\!\(\( | \^2\)\)= "", 
        b2abs, "  =  ", 
        b2abs /. Ω -> \(1\/2\) \@\(Δω\^2 + \((ω\_1\/2)\)\^2\)];\)\)

The next plots show the analytic Rabi frequency \(TraditionalForm\`ω\_R = \(2  Ω = \@\(Δω\^2 + \((ω\_1\/2)\)\^2\)\)\) as a function of  ω and (Subscript[ω, 1])and compare it to our numerical results:

\(\(Plot3D[\@\(\((ω00 - ω)\)\^2 + \((ω\_1\/2)\)\^2\), {ω, 0.8\ ω00, 
          1.2  ω00}, {ω\_1, \(-0.3\) ω00, 
          0.3  ω00}, PlotLabel -> "Analytic Rabi Frequency", 
        AxesLabel -> {\*""ω/(Subscript[ω, 0])"", \[RawEscape] \*""(Subscript[ω, 1])/(Subscript[ω, 0])"", \*""(Subscript[ω, R])/(Subscript[ω, 0])""}];\)\)
\(\(plFreq2 = 
        MultipleListPlot[Transpose[colldata[\([All, All, {2, 3}]\)]], 
          PlotStyle -> {GrayLevel[0.1], Hue[0.3], Hue[0.6], Hue[1]}, 
          SymbolStyle -> {GrayLevel[0], Hue[0.3], Hue[0.6], Hue[1]}, 
          PlotJoined -> False, 
          AxesLabel -> {\*""ω/(Subscript[ω, 0])"", \*""(Subscript[ω, R])/(Subscript[ω, 0])""}, 
          DisplayFunction -> Identity];\)\)
Print["\nComparison to numerical results:"];
\(\(Show[{plFreq2, 
          Plot[Evaluate[\(Function[ω1, \@\(\((ω00 - ω)\)\^2 + \((ω1\/2)\)\^2\)]\)[ω1list]], {ω, 0.8\ ω00, 
              1.2\ ω00}, 
            PlotStyle -> {GrayLevel[0], Hue[0.3], Hue[0.6], Hue[1]}, 
            DisplayFunction -> Identity]}, PlotLabel -> "", 
        DisplayFunction -> $DisplayFunction];\)\)

The symbols give the numerical results, the curves the analytic approximation. There are small differences for larger ω - (Subscript[ω, 0])and (Subscript[ω, 1]).

The next plots show the amplitude of  the oscillation in \(TraditionalForm\`\(\(|\)\(a(t)\)\( | \^2\)\)\), i. e. \(TraditionalForm\`1\/2 - Δω\^2\/\(2  ω\_R\^2\)\), and again compare the analytical and the numerical results:

\(Plot3D[
      1/2 - \(1\/2\) \((ω00 - ω)\)\^2/\((\((ω00 - ω)\)\^2 + \((ω\_1\/2)\)\^2)\), {ω, 0.8  ω00, 
        1.2\ ω00}, {ω\_1, 0.0001, 0.3  ω00}, 
      PlotLabel -> "Amplitude of Rabi Oscillation", 
      AxesLabel -> {\*""ω/(Subscript[ω, 0])"", \[RawEscape] \*""(Subscript[ω, 1])/(Subscript[ω, 0])"", "Ampl"}]; 
    plAmpl = MultipleListPlot[Transpose[colldata[\([All, All, {2, 4}]\)]], 
        PlotStyle -> {GrayLevel[0.1], Hue[0.3], Hue[0.6], Hue[1]}, 
        SymbolStyle -> {GrayLevel[0.1], Hue[0.3], Hue[0.6], Hue[1]}, 
        PlotJoined -> False, 
        DisplayFunction -> Identity];\)
Print["\nComparison to numerical results:"];
\(\(Show[{plAmpl, 
          Plot[\ Evaluate[\ \(Function[ω1, 
                  1/2 - \(1\/2\) \((ω00 - ω)\)\^2/\((\((ω00 - ω)\)\^2 + \((ω1\/2)\)\^2)\)]\)[ω1list]\ ], {ω, 0.8\ ω00, 
              1.2\ ω00}, 
            PlotStyle -> {GrayLevel[0], Hue[0.3], Hue[0.6], Hue[1]}, 
            DisplayFunction -> Identity\ ]}, 
        DisplayFunction -> $DisplayFunction, 
        AxesLabel -> {\*""ω/(Subscript[ω, 0])"", "Ampl"}];\)\)

As before, we notice small differences for larger ω - (Subscript[ω, 0])and (Subscript[ω, 1]).