Rapid Adiabatic Passage

Define the Schroedinger Equation and the wave function in the 2-state basis

Clear["Global`*"]
SE[H_,Ψ_]:=I (∂_t)Ψ==H. Ψ;
MatrixForm[Ψ[t_]={e[t],g[t]}];
H=0.5(({{-delta, Omega}, {Conjugate[Omega], delta}}));
MatrixForm[eqs={LogicalExpand[SE[H,Ψ[t]]][[1]],LogicalExpand[SE[H,Ψ[t]]][[2]]}]
(*Initial condition:all atoms are in the ground state at t=0, b: optically excited state, g: final state *)
MatrixForm[Ψ0={0,1}]; 
inits={LogicalExpand[Ψ[0]==Ψ0][[1]],LogicalExpand[Ψ[0]==Ψ0][[2]]}

Solve the DGL in a RAP scenario

OmegaUsed=2*Pi*5*10^3; (* s^-1 *)
deltaMax=2Pi*10^5; (* s^-1 *)

SweepDuration=10^-2; (* s *)

deltaFunction[t_,Duration_]:=deltaMax*(1-2t/Duration);

sols=NDSolve[Union[eqs/.{delta-> deltaFunction[t,SweepDuration],Omega->OmegaUsed},inits],{e,g},{t,0,SweepDuration},MaxSteps->10000][[1]];
Print["Eigenvalues:"];
Eigenvalues[H]
Plot[Eigenvalues[H]/.{Omega->OmegaUsed},{delta,-deltaMax,deltaMax}]
Plot[{Abs[g[t]]^2/.sols,Abs[e[t]]^2/.sols},{t,0,SweepDuration},PlotRange->{{0,SweepDuration},{0,1}},AxesLabel-> {t,Probability}]
Plot[{deltaFunction[t,SweepDuration]},{t,0,SweepDuration},AxesLabel-> {t,delta}]

transfer efficiency of rapid adiabatic passage in dependence of sweep duration

MinDuration=10^-5;
MaxDuration=10^-2;
NumSteps=50;
TransferProbability=Table[0,{NumSteps+1},{2}];
For[i=1,i≤(NumSteps+1),i++,
Duration=MinDuration+(MaxDuration-MinDuration)*(i-1)/NumSteps;
TransferProbability[[i,1]]=Duration;
sols=NDSolve[Union[eqs/.{delta-> deltaFunction[t,Duration],Omega->OmegaUsed},inits],{e,g},{t,0,Duration},MaxSteps->10000][[1]];
TransferProbability[[i,2]]=((Abs[e[Duration]])^2)/.sols;
];
ListPlot[{TransferProbability},Joined->True,PlotRange->{{MinDuration,MaxDuration},{0,1}},AxesLabel-> {"Duration","Transfer probability"}]

Eg[Om_,delta_]:=-0.5delta +0.5*(-delta+Sign[delta]Sqrt[Om^2+delta^2]);
Ee[Om_,delta_]:=0.5delta+0.5*(-delta-Sign[delta]Sqrt[Om^2+delta^2]);


Egd[Om_,delta_]:=Eg[Om,delta]+delta;
Eed[Om_,delta_]:=Ee[Om,delta];

Egdu[Om_,delta_]:=Eg[0,delta]+delta;
Eedu[Om_,delta_]:=Ee[0,delta];

Egs[Om_,delta_]:=If[delta<0,Egd[Om,delta],Eed[Om,delta]];
Ees[Om_,delta_]:=If[delta<0,Eed[Om,delta],Egd[Om,delta]];

GroundStateProbability[Om_,delta_]:=(Normalize[Eigenvectors[(({{-delta, Om}, {Om, delta}}))][[1]]][[1]])^2;
ExcitedStateProbability[Om_,delta_]:=1-GroundStateProbability[Om,delta];

Plotg=Plot[Egs[1,delta],{delta,-3,3},PlotRange->{{-3.4,3.4},{-2,2}},AspectRatio-> 0.8,ColorFunctionScaling->False,ColorFunction->Function[{delta},Blend[ColorConvert[{Blue,Red},"RGB"],GroundStateProbability[1,delta]]],PlotStyle->{Thickness[0.01]},FrameLabel->{time, probability},AxesStyle->{Directive[Thick],Directive[Thick]},TicksStyle->Directive[Thick],LabelStyle->Directive[Large],AxesOrigin->{-3,-1.5}];

Plote=Plot[Ees[1,delta],{delta,-3,3},PlotRange->{{-3.4,3.4},{-2,2}},AspectRatio-> 0.8,ColorFunctionScaling->False,ColorFunction->Function[{delta},Blend[ColorConvert[{Blue,Red},"RGB"],ExcitedStateProbability[1,delta]]],PlotStyle->{Thickness[0.01]},FrameLabel->{time, probability},AxesStyle->{Directive[Thick],Directive[Thick]},TicksStyle->Directive[Thick],LabelStyle->Directive[Large],AxesOrigin->{-3,-1.5}];

Plotb=Plot[{Egdu[1,delta],Eedu[1,delta]},{delta,-3,3} ,PlotRange->{{-3.4,3.4},{-2,2}},AspectRatio-> 0.8,PlotStyle->{{Blue,Dashed,Thickness[0.005]},{Red,Dashed,Thickness[0.005]}},FrameLabel->{time, probability},AxesStyle->{Directive[Thick],Directive[Thick]},TicksStyle->Directive[Thick],LabelStyle->Directive[Large],AxesOrigin->{-3,-1.5}];

Show[Plotg,Plote,Plotb]

Initial condition: all atoms are in the ground state at t=0

MatrixForm[Ψ0={1,0,0}]
inits={LogicalExpand[Ψ[0]==Ψ0][[1]],LogicalExpand[Ψ[0]==Ψ0][[2]],LogicalExpand[Ψ[0]==Ψ0][[3]]};

Now we solve the DGL in a STIRAP scenario explicitly. Here is the place where we put numerical values for various parameters.

Omega1Max=2*Pi*110*10^3;
Omega2Max=2*Pi*110*10^3;
DeltaT=100*10^-6;
Intensity1[t_]:=If[t<0,0,If[t>DeltaT,1,0.5(1-Cos[Pi t/DeltaT])]];
Intensity2[t_]:=If[t<0,1,If[t>DeltaT,0,0.5(1+Cos[Pi t/DeltaT])]];
Omega1Function[t_]:=Omega1Max*Sqrt[Intensity1[t]];
Omega2Function[t_]:=Omega2Max*Sqrt[Intensity2[t]];
deltaFunction[t_]:=Intensity1[t]*50000-Intensity2[t]*25000;

sols=NDSolve[Union[eqs/.{delta-> deltaFunction[t], Delta->0,Γa->0/(100*10^-6), Γg->0/(100*10^-6),Γb->1/(10*10^-6),Omega1->Omega1Function[t],Omega2->Omega2Function[t]},inits],{a,b,g},{t,0,DeltaT},MaxSteps->10000][[1]];

Now we plot the solution. Here we show the probability of finding the system in the excited state.

Plot[{((Abs[b[t]])^2)/.sols,((Abs[g[t]])^2)/.sols,((Abs[a[t]])^2)/.sols},{t,0,DeltaT},PlotRange->{{0,DeltaT},{0,1}},AxesLabel-> {t,Probability}]
Plot[{Omega1Function[t],Omega2Function[t]},{t,0,DeltaT},AxesLabel-> {t,Omega}]

Now we solve the DGL in a dark resonance scenario explicitly.

(* first we cross check the DGL solution with analytical solution for special case that Omega2=0 and we treat a free-free transition *)
(* constants *)
kB=1.3806503*10^-23;
a0=0.052917721086*10^-9;
hbar=1.05457148*10^-34;
h=2*Pi*hbar;
massunit=1.66053886*10^-27;
massLi=6*massunit;
massK=40*massunit;
massSr=87*massunit;
c=299792458;
kb=1.38 10^-23;
(* analytical description of scattering, equations from Metcalf *)
PowerLaser=0.00007*10^-6;  (* W *)
Waist=100*10 ^-6;
Intensity=N[2*PowerLaser/(Pi*Waist^2)];
Lambda=689.9267499282614*10^-9;(*84Sr*)
LifetimeExcitedState=21.2*10^-6;
Linewidth=1/LifetimeExcitedState;
SaturationIntensity[Lambda_,Linewidth_]:=Pi*h*c/(3*Lambda^3/Linewidth);
Detuning[Deltaf_]:=2*Pi*Deltaf;
LambdaLaser[Deltaf_]:=c/(c/Lambda+Deltaf);
s0=Intensity/SaturationIntensity[Lambda,Linewidth];
ScatteringRate[Detuning_ ]:=s0*(Linewidth/2)/(1+s0+(2*Detuning/Linewidth)^2);
RabiFrequency=Sqrt[s0/2]*Linewidth;
Print["Linewidth = 2 Pi * ",N[Linewidth/(2*Pi)]," Hz"];
Print["RabiFrequency = 2 Pi * ",N[RabiFrequency/(2*Pi)]," Hz"];
Print["Intensity = ",N[Intensity*1000/10000]," mW/cm^2 = ",Intensity/SaturationIntensity[Lambda,Linewidth] ," ISat"];
Print["ScatteringRate on resonance = ",N[ScatteringRate[0]]," s^-1"];

(* compare analytic result and DGL solution: loss spectrum *)
Gammab=Linewidth;
Omega1Max=RabiFrequency;
Omega2Max=0;
Duration=10*10^-3;
DetuningRange=240*10^3;

FrequencyUnit=1000;
NumSteps=200;
MatrixForm[Ψ0={1,0,0}];
inits={LogicalExpand[Ψ[0]==Ψ0][[1]],LogicalExpand[Ψ[0]==Ψ0][[2]],LogicalExpand[Ψ[0]==Ψ0][[3]]};

RemainingAtoms=Table[0,{NumSteps+1},{2}];
RemainingAtomsAnalytical=Table[0,{NumSteps+1},{2}];
For[i=1,i≤(NumSteps+1),i++,
Deltaf=-DetuningRange/2+DetuningRange*(i-1)/NumSteps;
RemainingAtoms[[i,1]]=Deltaf/FrequencyUnit;
RemainingAtomsAnalytical[[i,1]]=RemainingAtoms[[i,1]];
RemainingAtomsAnalytical[[i,2]]=Exp[-ScatteringRate[Detuning[Deltaf]]*Duration];

sols=NDSolve[Union[eqs/.{delta-> 2*Pi*Deltaf, Delta->0,Γa->0, Γg->0,Γb->Gammab,Omega1->Omega1Max,Omega2->Omega2Max},inits],{a,b,g},{t,0,Duration},MaxSteps->100000][[1]];
RemainingAtoms[[i,2]]=((Abs[a[Duration]])^2)/.sols;
];
ListPlot[{RemainingAtoms,RemainingAtomsAnalytical},Joined->True,PlotRange->{{-0.5*DetuningRange/FrequencyUnit,0.5*DetuningRange/FrequencyUnit},{0,1}}]

(* now we fit to experimental data of a dark resonance *)
DarkResonanceData12Feb19S65={{-0.100000000000023,0.908439249881415},{-0.089999999999975,0.922179647893438},{-0.0799999999999841,0.898070457859241},{-0.0699999999999932,0.853938983988036},{-0.0600000000000023,0.835849470773183},{-0.0500000000000114,0.736937600717607},{-0.0400000000000205,0.623642959560419},{-0.0299999999999727,0.421305089998592},{-0.0199999999999818,0.195173324148714},{-0.0099999999999909,0.141235808097285},{-0.0099999999999909,0.131526897384606},{-0.00778000000002521,0.164174847972017},{-0.007000000000005,0.213354083812871},{-0.00599999999997181,0.371127248550499},{-0.00556000000000267,0.454769383240778},{-0.00499999999999545,0.638222985388154},{-0.0040000000000191,0.833744412579236},{-0.00333999999998014,0.719659212502997},{-0.0029999999999859,0.719576846906611},{-0.00200000000000955,0.520946648593524},{-0.00112000000001444,0.339457132021479},{-.000999999999976353,0.316574758815619},{0,0.236288416579219},{0.00112000000001444,0.188322558642563},{0.00333999999998014,0.121654609413886},{0.00556000000000267,0.109078066763101},{0.00778000000002521,0.123421121537304},{0.0099999999999909,0.14886563992691},{0.0199999999999818,0.258329628092483},{0.0299999999999727,0.427279517540086},{0.0400000000000205,0.59136970411884},{0.0500000000000114,0.733961449201246},{0.0600000000000023,0.76828588290647},{0.0699999999999932,0.79629211742032},{0.0799999999999841,0.840748579854445},{0.089999999999975,0.822059532346738},{0.100000000000023,0.809662854171869}};
OffsetD=0.1;
ScaleD=1.15;
ShiftF=0.0035;
For[i=1,i≤(Length[DarkResonanceData12Feb19S65]),i++,
DarkResonanceData12Feb19S65[[i,2]]=ScaleD*(DarkResonanceData12Feb19S65[[i,2]]-OffsetD);
DarkResonanceData12Feb19S65[[i,1]]=DarkResonanceData12Feb19S65[[i,1]]+ShiftF;
];
ListPlot[DarkResonanceData12Feb19S65,Joined->True,PlotRange->{{-0.1,0.1},{0,1}}];
ListPlot[DarkResonanceData12Feb19S65,Joined->True,PlotRange->{{-0.01,0.01},{0,1}}];

(* Compare exp data and DGL solution : dark resonance *)
Gammab=2*1/(21*10^-6);
Omega1Max=2*Pi*0.5*10^3;
Omega2Max=2*Pi*15*10^3;
Duration=100*10^-3;
DetuningRange=200*10^3;

FrequencyUnit=1000000;
NumSteps=Length[DarkResonanceData12Feb19S65];
MatrixForm[Ψ0={1,0,0}];
inits={LogicalExpand[Ψ[0]==Ψ0][[1]],LogicalExpand[Ψ[0]==Ψ0][[2]],LogicalExpand[Ψ[0]==Ψ0][[3]]};

RemainingAtoms=Table[0,{NumSteps+1},{2}];
For[i=1,i≤(NumSteps+1),i++,
Deltaf=FrequencyUnit*DarkResonanceData12Feb19S65[[i,1]];
RemainingAtoms[[i,1]]=Deltaf/FrequencyUnit;
sols=NDSolve[Union[eqs/.{delta-> 2*Pi*Deltaf, Delta->0,Γa->0, Γg->0,Γb->Gammab,Omega1->Omega1Max,Omega2->Omega2Max},inits],{a,b,g},{t,0,Duration},MaxSteps->100000][[1]];
RemainingAtoms[[i,2]]=((Abs[a[Duration]])^2)/.sols;
];
ListPlot[{RemainingAtoms,DarkResonanceData12Feb19S65},Joined->True,PlotRange->{{-0.5*DetuningRange/FrequencyUnit,0.5*DetuningRange/FrequencyUnit},{0,1}}]
ListPlot[{RemainingAtoms,DarkResonanceData12Feb19S65},Joined->True,PlotRange->{{-0.05*DetuningRange/FrequencyUnit,0.05*DetuningRange/FrequencyUnit},{0,1}}]

AutlerTownesData12Mar20S5={{-0.199999999999989,0.940589800727273},{-0.196000000000026,0.947047261818182},{-0.192000000000007,0.925991504545455},{-0.187999999999988,0.943233340727273},{-0.184000000000026,0.911101347636364},{-0.180000000000007,0.934031435090909},{-0.175999999999988,0.917597164363636},{-0.168000000000006,0.907003632727273},{-0.163999999999987,0.900238250363636},{-0.160000000000025,0.870202160545455},{-0.156000000000006,0.863466479272727},{-0.151999999999987,0.845782977454546},{-0.148000000000025,0.830959618909091},{-0.144000000000005,0.811246014363636},{-0.139999999999986,0.784414961090909},{-0.136000000000024,0.755340091090909},{-0.132000000000005,0.730036280727273},{-0.127999999999986,0.679576372181818},{-0.124000000000024,0.623967934363636},{-0.120000000000005,0.560284338181818},{-0.115999999999985,0.460128844545455},{-0.112000000000023,0.340173326909091},{-0.108000000000004,0.196623148545455},{-0.103999999999985,0.109249632490909},{-0.100000000000023,0.0686682947090909},{-0.0960000000000036,0.0758762838181818},{-0.0919999999999845,0.148656510472727},{-0.0880000000000223,0.282460934181818},{-0.0840000000000032,0.426433764181818},{-0.0799999999999841,0.556809885636364},{-0.0760000000000218,0.673040230727273},{-0.0720000000000027,0.746171895818182},{-0.0679999999999836,0.806203658727273},{-0.0640000000000214,0.832318002545455},{-0.0599399999999832,0.860039544},{-0.0558800000000019,0.881562866181818},{-0.0518200000000206,0.905923151090909},{-0.0477599999999825,0.914674379454546},{-0.0437000000000012,0.914878020909091},{-0.0396400000000199,0.935634133818182},{-0.0355599999999754,0.937531256363636},{-0.0314999999999941,0.945996434727273},{-0.0274400000000128,0.952906369272727},{-0.0233799999999746,0.946855972181818},{-0.0193199999999933,0.953228433818182},{-0.015260000000012,0.944693160545455},{-0.0111999999999739,0.962087720181818},{-0.0071399999999926,0.962779019090909},{-0.0030800000000113,0.972464050363636},{0.000980000000026848,0.952503135090909},{0.00504000000000815,0.959532109272727},{0.00909999999998945,0.967030720181818},{0.0131600000000276,0.960369719454545},{0.0172400000000152,0.976896597272727},{0.0212999999999965,0.972327459272727},{0.0253599999999778,0.967669529454546},{0.029420000000016,0.974392699272727},{0.0334799999999973,0.979569049636364},{0.0375399999999786,0.962940977454546},{0.0416000000000167,0.960571804363636},{0.045659999999998,0.962702910727273},{0.0497199999999793,0.957174310363636},{0.0537800000000175,0.946768079818182},{0.0578399999999988,0.939004445636364},{0.0618999999999801,0.924957286},{0.0659600000000182,0.916351569454545},{0.0700400000000059,0.923690990545455},{0.0740999999999872,0.902570753272727},{0.0781600000000253,0.887842651272727},{0.0822200000000066,0.854387133636364},{0.0862799999999879,0.843375816363636},{0.0903400000000261,0.800893603454546},{0.0944000000000074,0.740023323090909},{0.0984599999999887,0.641596388909091},{0.106580000000008,0.355410246363636},{0.110639999999989,0.201819004363636},{0.114700000000028,0.0784422361272727},{0.118760000000009,0.0637050141454545},{0.122839999999997,0.107284433672727},{0.126899999999978,0.220040348727273},{0.130960000000016,0.312184838181818},{0.135019999999997,0.430495146363636},{0.139079999999979,0.522828336545455},{0.143140000000017,0.603173072909091},{0.147199999999998,0.672594548181818},{0.151259999999979,0.721936195090909},{0.155320000000017,0.754875480363636},{0.159379999999999,0.801886776545455},{0.16343999999998,0.823649329454545},{0.167500000000018,0.84484513},{0.172000000000025,0.867375642181818},{0.175999999999988,0.886104086545455},{0.180000000000007,0.896723358363636},{0.184000000000026,0.902960345454546},{0.187999999999988,0.894948918363636},{0.192000000000007,0.917788198363636},{0.196000000000026,0.929851369090909},{0.199999999999989,0.929135370545455}};
OffsetD=0.06;
ScaleD=1.1;
ShiftF=-0.009;
For[i=1,i≤(Length[AutlerTownesData12Mar20S5]),i++,
AutlerTownesData12Mar20S5[[i,2]]=ScaleD*(AutlerTownesData12Mar20S5[[i,2]]-OffsetD);
AutlerTownesData12Mar20S5[[i,1]]=AutlerTownesData12Mar20S5[[i,1]]+ShiftF;
];
ListPlot[AutlerTownesData12Mar20S5,Joined->True,PlotRange->{{-0.2,0.2},{0,1}}];

(* Autler-Townes splitting *)
Gammab=2*1/(21*10^-6);
Omega1Max=2*Pi*0.6*10^3;
Omega2Max=2*Pi*218*10^3;
Duration=100*10^-3;
DetuningRange=400*10^3;

FrequencyUnit=1000000;
NumSteps=Length[AutlerTownesData12Mar20S5];
MatrixForm[Ψ0={1,0,0}];
inits={LogicalExpand[Ψ[0]==Ψ0][[1]],LogicalExpand[Ψ[0]==Ψ0][[2]],LogicalExpand[Ψ[0]==Ψ0][[3]]};

RemainingAtoms=Table[0,{NumSteps+1},{2}];
For[i=1,i≤(NumSteps+1),i++,
Deltaf=FrequencyUnit*AutlerTownesData12Mar20S5[[i,1]];
RemainingAtoms[[i,1]]=Deltaf/FrequencyUnit;
sols=NDSolve[Union[eqs/.{delta-> 2*Pi*Deltaf, Delta->0,Γa->0, Γg->0,Γb->Gammab,Omega1->Omega1Max,Omega2->Omega2Max},inits],{a,b,g},{t,0,Duration},MaxSteps->100000][[1]];
RemainingAtoms[[i,2]]=((Abs[a[Duration]])^2)/.sols;
];
ListPlot[{RemainingAtoms,AutlerTownesData12Mar20S5},Joined->True,PlotRange->{{-0.5*DetuningRange/FrequencyUnit,0.5*DetuningRange/FrequencyUnit},{0,1}}]

(* normalize Rabi frequencies by Sqrt(Intensity) *)
PowerLaser1=7*10^-6;  (* W *)
Waist=105*10 ^-6;
Intensity1=N[2*PowerLaser1/(Pi*Waist^2)];
PowerLaser2=2.8*10^-3;  (* W *)
Waist=100*10 ^-6;
Intensity2=N[2*PowerLaser2/(Pi*Waist^2)];
Print["Omega1/Sqrt[Intensity1] = ",Omega1Max/Intensity1];
Print["Omega2/Sqrt[Intensity2] = ",Omega2Max/Intensity2];

(* decay measurement of Omega 1 *)
Omega1Max=2*Pi*4*10^3; (* for 0.127 mW, 127µs *)
Omega2Max=2*Pi*0*10^3;
Gammab=2*1/(21*10^-6);
DeltaT=0.5*10^-3;
DecayTime=0.127*10^-3;
sols=NDSolve[Union[eqs/.{delta-> 0, Delta->0,Γa->0/(100*10^-6), Γg->0/(100*10^-6),Γb->Gammab,Omega1->Omega1Max,Omega2->Omega2Max},inits],{a,b,g},{t,0,DeltaT},MaxSteps->10000][[1]];
Plot[{((Abs[a[t]])^2)/.sols,Exp[-t/DecayTime]},{t,0,DeltaT},PlotRange->{{0,DeltaT},{0,1}},AxesLabel-> {t,Probability}]

PowerLaser1=127*10^-6;  (* W *)
Waist=105*10 ^-6;
Intensity1=N[2*PowerLaser1/(Pi*Waist^2)];
Print["Omega1/Sqrt[Intensity1] = ",Omega1Max/Intensity1];

(* fitting the Autler-Townes splitting *)
Gammab=2*1/(21*10^-6);
Duration=100*10^-3;
DetuningRange=400*10^3;

NumSteps=Length[AutlerTownesData12Mar20S5];
FrequencyUnit=1000000;
Data=AutlerTownesData12Mar20S5;

MatrixForm[Ψ0={1,0,0}];
inits={LogicalExpand[Ψ[0]==Ψ0][[1]],LogicalExpand[Ψ[0]==Ψ0][[2]],LogicalExpand[Ψ[0]==Ψ0][[3]]};

FitModel[Omega1Max_?NumberQ,Omega2Max_?NumberQ,offset_?NumberQ,scale_?NumberQ,fshift_?NumberQ,Df_?NumberQ]:=(
sols=NDSolve[Union[eqs/.{delta-> 2*Pi*(Df-fshift)*FrequencyUnit, Delta->0,Γa->0, Γg->0,Γb->Gammab,Omega1->Omega1Max,Omega2->Omega2Max},inits],{a,b,g},{t,0,Duration},MaxSteps->100000][[1]];
FitModel[Omega1Max,Omega2Max,offset,scale,fshift,Df]=scale*(offset+((Abs[a[Duration]])^2)/.sols);
)

Fitmodel[Omega1Max_?NumberQ,Omega2Max_?NumberQ,offset_?NumberQ,scale_?NumberQ,fshift_?NumberQ]:=(Fitmodel[Omega1Max,Omega2Max,offset,scale,fshift]=First[x/.NDSolve[{x''[t]+γ x'[t]+a x[t]+b x[t]^3==0,x[0]==2,x'[0]==0},x,{t,20}]]);

MinOmega1=2*Pi*0.1*10^3;
StartOmega1=2*Pi*0.6*10^3;
MaxOmega1=2*Pi*2*10^3;

MinOmega2=2*Pi*210*10^3;
StartOmega2=2*Pi*218*10^3;
MaxOmega2=2*Pi*225*10^3;

Minoffset=-0.1;
Startoffset=0;
Maxoffset=0.1;

Minscale=0.8;
Startscale=1;
Maxscale=1.2;

Minfshift=-10000;
Startfshift=0;
Maxfshift=10000;

FitResult=FindFit[Data,{FitModel[Omega1f,Omega2f,offsetf,scalef,fshift,Df],{MinOmega1≤Omega1f≤MaxOmega1,MinOmega2≤Omega2f≤MaxOmega2,Minoffset≤offsetf≤Maxoffset,Minscale≤scalef≤Maxscale,Minfshift≤fshift≤Maxfshift}},{{Omega1f,StartOmega1},{Omega2f,StartOmega2},{offsetf,Startoffset},{scalef,Startscale},{fshift,Startfshift}},Df];



Result=Table[0,{NumSteps+1},{2}];
For[i=1,i≤(NumSteps+1),i++,
Deltaf=AutlerTownesData12Mar20S5[[i,1]];
Result[[i,1]]=Deltaf;
Result[[i,2]]=FitModel[Omega1,Omega2,offset,scale,fshift,Deltaf]/.FitResult;
];


ListPlot[{Data,Result}]

FitModel[Omega1Max_?NumberQ,Omega2Max_?NumberQ,offset_?NumberQ,scale_?NumberQ,fshift_?NumberQ,Deltaf_?NumberQ]:=(
sols=NDSolve[Union[eqs/.{delta-> 2*Pi*(Deltaf-fshift)*FrequencyUnit, Delta->0,Γa->0, Γg->0,Γb->Gammab,Omega1->Omega1Max,Omega2->Omega2Max},inits],{a,b,g},{t,0,Duration},MaxSteps->100000][[1]];
scale*(offset+((Abs[a[Duration]])^2)/.sols)
)
FitModel[StartOmega1,StartOmega2,Startoffset,Startscale,Startfshift,0]

(* test mathematica curve fitting *)
Clear["Global`*"]
Data=Table[{x,x^2},{x,0,20}];
model[a_,o_,exponent_,x_]:=a x^exponent+o;
Mina=0;
Starta=0.1;
Maxa=10;
Mino=-10;
Starto=1;
Maxo=10;
Mine=0.1;
Starte=2;
Maxe=10;
FitResult=FindFit[Data,{model[a,o,exponent,x],{Mina≤a≤Maxa,Mino≤o≤Maxo,Mine≤exponent≤Maxe}},{{a,Starta},{o,Starto},{exponent,Starte}},x]
Result=Table[{x,model[a]/.FitResult},{x,-20,20}];
ListPlot[{Data,Result}]

(* test Mathematica DLG fitting *)
sol=First[x/.NDSolve[{x''[t]+.33 x'[t]+.72x[t]^3==0,x[0]==2,x'[0]==0},x,{t,20}]];
times=N[Range[0,100]/5];
data=Transpose[{times,sol[times]+RandomReal[.1,101]}];
lp=ListPlot[data,PlotRange->All]
model[γ_?NumberQ,a_?NumberQ,b_?NumberQ]:=(model[γ,a,b]=First[x/.NDSolve[{x''[t]+γ x'[t]+a x[t]+b x[t]^3==0,x[0]==2,x'[0]==0},x,{t,20}]]);
fit=FindFit[data,model[γ,a,b][x],{{γ,.1},{a,.1},{b,1}},x,PrecisionGoal->4,AccuracyGoal->4]
{Show[Plot[model[γ,a,b][x]/.fit,{x,0,20}],lp],ListPlot[Transpose[{times,data[[All,2]]-(model[γ,a,b][times]/.fit)}]]}

(* test Mathematica DLG fitting *)
sol=First[x/.NDSolve[{x''[t]+.33 x'[t]+.72x[t]^3==0,x[0]==2,x'[0]==0},x,{t,20}]];
times=N[Range[0,100]/5];
data=Transpose[{times,sol[times]+RandomReal[.1,101]}];
lp=ListPlot[data,PlotRange->All]
model[γ_?NumberQ,a_?NumberQ,b_?NumberQ]:=(1+1;model[γ,a,b]=First[x/.NDSolve[{x''[t]+γ x'[t]+a x[t]+b x[t]^3==0,x[0]==2,x'[0]==0},x,{t,20}]]);
fit=FindFit[data,model[γ,a,b][x],{{γ,.1},{a,.1},{b,1}},x,PrecisionGoal->4,AccuracyGoal->4]
{Show[Plot[model[γ,a,b][x]/.fit,{x,0,20}],lp],ListPlot[Transpose[{times,data[[All,2]]-(model[γ,a,b][times]/.fit)}]]}

(* compare analytic result and DGL solution: loss spectrum *)
Gammab=Linewidth;
Omega1Max=RabiFrequency;
Omega2Max=0;
Duration=10*10^-3;
DetuningRange=240*10^3;

FrequencyUnit=1000;
NumSteps=200;
MatrixForm[Ψ0={1,0,0}];
inits={LogicalExpand[Ψ[0]==Ψ0][[1]],LogicalExpand[Ψ[0]==Ψ0][[2]],LogicalExpand[Ψ[0]==Ψ0][[3]]};

RemainingAtoms=Table[0,{NumSteps+1},{2}];
RemainingAtomsAnalytical=Table[0,{NumSteps+1},{2}];
For[i=1,i≤(NumSteps+1),i++,
Deltaf=-DetuningRange/2+DetuningRange*(i-1)/NumSteps;
RemainingAtoms[[i,1]]=Deltaf/FrequencyUnit;
RemainingAtomsAnalytical[[i,1]]=RemainingAtoms[[i,1]];
RemainingAtomsAnalytical[[i,2]]=Exp[-ScatteringRate[Detuning[Deltaf]]*Duration];

sols=NDSolve[Union[eqs/.{delta-> 2*Pi*Deltaf, Delta->0,Γa->0, Γg->0,Γb->Gammab,Omega1->Omega1Max,Omega2->Omega2Max},inits],{a,b,g},{t,0,Duration},MaxSteps->100000][[1]];
RemainingAtoms[[i,2]]=((Abs[a[Duration]])^2)/.sols;
];
ListPlot[{RemainingAtoms,RemainingAtomsAnalytical},Joined->True,PlotRange->{{-0.5*DetuningRange/FrequencyUnit,0.5*DetuningRange/FrequencyUnit},{0,1}}]

(* now we fit to experimental data of a dark resonance *)
DarkResonanceData12Feb19S65={{-0.100000000000023,0.908439249881415},{-0.089999999999975,0.922179647893438},{-0.0799999999999841,0.898070457859241},{-0.0699999999999932,0.853938983988036},{-0.0600000000000023,0.835849470773183},{-0.0500000000000114,0.736937600717607},{-0.0400000000000205,0.623642959560419},{-0.0299999999999727,0.421305089998592},{-0.0199999999999818,0.195173324148714},{-0.0099999999999909,0.141235808097285},{-0.0099999999999909,0.131526897384606},{-0.00778000000002521,0.164174847972017},{-0.007000000000005,0.213354083812871},{-0.00599999999997181,0.371127248550499},{-0.00556000000000267,0.454769383240778},{-0.00499999999999545,0.638222985388154},{-0.0040000000000191,0.833744412579236},{-0.00333999999998014,0.719659212502997},{-0.0029999999999859,0.719576846906611},{-0.00200000000000955,0.520946648593524},{-0.00112000000001444,0.339457132021479},{-.000999999999976353,0.316574758815619},{0,0.236288416579219},{0.00112000000001444,0.188322558642563},{0.00333999999998014,0.121654609413886},{0.00556000000000267,0.109078066763101},{0.00778000000002521,0.123421121537304},{0.0099999999999909,0.14886563992691},{0.0199999999999818,0.258329628092483},{0.0299999999999727,0.427279517540086},{0.0400000000000205,0.59136970411884},{0.0500000000000114,0.733961449201246},{0.0600000000000023,0.76828588290647},{0.0699999999999932,0.79629211742032},{0.0799999999999841,0.840748579854445},{0.089999999999975,0.822059532346738},{0.100000000000023,0.809662854171869}};
OffsetD=0.1;
ScaleD=1.15;
ShiftF=0.0035;
For[i=1,i≤(Length[DarkResonanceData12Feb19S65]),i++,
DarkResonanceData12Feb19S65[[i,2]]=ScaleD*(DarkResonanceData12Feb19S65[[i,2]]-OffsetD);
DarkResonanceData12Feb19S65[[i,1]]=DarkResonanceData12Feb19S65[[i,1]]+ShiftF;
];
ListPlot[DarkResonanceData12Feb19S65,Joined->True,PlotRange->{{-0.1,0.1},{0,1}}];
ListPlot[DarkResonanceData12Feb19S65,Joined->True,PlotRange->{{-0.01,0.01},{0,1}}];

(* Compare exp data and DGL solution : dark resonance *)
Gammab=2*1/(21*10^-6);
Omega1Max=2*Pi*0.5*10^3;
Omega2Max=2*Pi*15*10^3;
Duration=100*10^-3;
DetuningRange=200*10^3;

FrequencyUnit=1000000;
NumSteps=Length[DarkResonanceData12Feb19S65];
MatrixForm[Ψ0={1,0,0}];
inits={LogicalExpand[Ψ[0]==Ψ0][[1]],LogicalExpand[Ψ[0]==Ψ0][[2]],LogicalExpand[Ψ[0]==Ψ0][[3]]};

RemainingAtoms=Table[0,{NumSteps+1},{2}];
For[i=1,i≤(NumSteps+1),i++,
Deltaf=FrequencyUnit*DarkResonanceData12Feb19S65[[i,1]];
RemainingAtoms[[i,1]]=Deltaf/FrequencyUnit;
sols=NDSolve[Union[eqs/.{delta-> 2*Pi*Deltaf, Delta->0,Γa->0, Γg->0,Γb->Gammab,Omega1->Omega1Max,Omega2->Omega2Max},inits],{a,b,g},{t,0,Duration},MaxSteps->100000][[1]];
RemainingAtoms[[i,2]]=((Abs[a[Duration]])^2)/.sols;
];
ListPlot[{RemainingAtoms,DarkResonanceData12Feb19S65},Joined->True,PlotRange->{{-0.5*DetuningRange/FrequencyUnit,0.5*DetuningRange/FrequencyUnit},{0,1}}]
ListPlot[{RemainingAtoms,DarkResonanceData12Feb19S65},Joined->True,PlotRange->{{-0.05*DetuningRange/FrequencyUnit,0.05*DetuningRange/FrequencyUnit},{0,1}}]

AutlerTownesData12Mar20S5={{-0.199999999999989,0.940589800727273},{-0.196000000000026,0.947047261818182},{-0.192000000000007,0.925991504545455},{-0.187999999999988,0.943233340727273},{-0.184000000000026,0.911101347636364},{-0.180000000000007,0.934031435090909},{-0.175999999999988,0.917597164363636},{-0.168000000000006,0.907003632727273},{-0.163999999999987,0.900238250363636},{-0.160000000000025,0.870202160545455},{-0.156000000000006,0.863466479272727},{-0.151999999999987,0.845782977454546},{-0.148000000000025,0.830959618909091},{-0.144000000000005,0.811246014363636},{-0.139999999999986,0.784414961090909},{-0.136000000000024,0.755340091090909},{-0.132000000000005,0.730036280727273},{-0.127999999999986,0.679576372181818},{-0.124000000000024,0.623967934363636},{-0.120000000000005,0.560284338181818},{-0.115999999999985,0.460128844545455},{-0.112000000000023,0.340173326909091},{-0.108000000000004,0.196623148545455},{-0.103999999999985,0.109249632490909},{-0.100000000000023,0.0686682947090909},{-0.0960000000000036,0.0758762838181818},{-0.0919999999999845,0.148656510472727},{-0.0880000000000223,0.282460934181818},{-0.0840000000000032,0.426433764181818},{-0.0799999999999841,0.556809885636364},{-0.0760000000000218,0.673040230727273},{-0.0720000000000027,0.746171895818182},{-0.0679999999999836,0.806203658727273},{-0.0640000000000214,0.832318002545455},{-0.0599399999999832,0.860039544},{-0.0558800000000019,0.881562866181818},{-0.0518200000000206,0.905923151090909},{-0.0477599999999825,0.914674379454546},{-0.0437000000000012,0.914878020909091},{-0.0396400000000199,0.935634133818182},{-0.0355599999999754,0.937531256363636},{-0.0314999999999941,0.945996434727273},{-0.0274400000000128,0.952906369272727},{-0.0233799999999746,0.946855972181818},{-0.0193199999999933,0.953228433818182},{-0.015260000000012,0.944693160545455},{-0.0111999999999739,0.962087720181818},{-0.0071399999999926,0.962779019090909},{-0.0030800000000113,0.972464050363636},{0.000980000000026848,0.952503135090909},{0.00504000000000815,0.959532109272727},{0.00909999999998945,0.967030720181818},{0.0131600000000276,0.960369719454545},{0.0172400000000152,0.976896597272727},{0.0212999999999965,0.972327459272727},{0.0253599999999778,0.967669529454546},{0.029420000000016,0.974392699272727},{0.0334799999999973,0.979569049636364},{0.0375399999999786,0.962940977454546},{0.0416000000000167,0.960571804363636},{0.045659999999998,0.962702910727273},{0.0497199999999793,0.957174310363636},{0.0537800000000175,0.946768079818182},{0.0578399999999988,0.939004445636364},{0.0618999999999801,0.924957286},{0.0659600000000182,0.916351569454545},{0.0700400000000059,0.923690990545455},{0.0740999999999872,0.902570753272727},{0.0781600000000253,0.887842651272727},{0.0822200000000066,0.854387133636364},{0.0862799999999879,0.843375816363636},{0.0903400000000261,0.800893603454546},{0.0944000000000074,0.740023323090909},{0.0984599999999887,0.641596388909091},{0.106580000000008,0.355410246363636},{0.110639999999989,0.201819004363636},{0.114700000000028,0.0784422361272727},{0.118760000000009,0.0637050141454545},{0.122839999999997,0.107284433672727},{0.126899999999978,0.220040348727273},{0.130960000000016,0.312184838181818},{0.135019999999997,0.430495146363636},{0.139079999999979,0.522828336545455},{0.143140000000017,0.603173072909091},{0.147199999999998,0.672594548181818},{0.151259999999979,0.721936195090909},{0.155320000000017,0.754875480363636},{0.159379999999999,0.801886776545455},{0.16343999999998,0.823649329454545},{0.167500000000018,0.84484513},{0.172000000000025,0.867375642181818},{0.175999999999988,0.886104086545455},{0.180000000000007,0.896723358363636},{0.184000000000026,0.902960345454546},{0.187999999999988,0.894948918363636},{0.192000000000007,0.917788198363636},{0.196000000000026,0.929851369090909},{0.199999999999989,0.929135370545455}};
OffsetD=0.06;
ScaleD=1.1;
ShiftF=-0.009;
For[i=1,i≤(Length[AutlerTownesData12Mar20S5]),i++,
AutlerTownesData12Mar20S5[[i,2]]=ScaleD*(AutlerTownesData12Mar20S5[[i,2]]-OffsetD);
AutlerTownesData12Mar20S5[[i,1]]=AutlerTownesData12Mar20S5[[i,1]]+ShiftF;
];
ListPlot[AutlerTownesData12Mar20S5,Joined->True,PlotRange->{{-0.2,0.2},{0,1}}];

(* Autler-Townes splitting *)
Gammab=2*1/(21*10^-6);
Omega1Max=2*Pi*0.6*10^3;
Omega2Max=2*Pi*218*10^3;
Duration=100*10^-3;
DetuningRange=400*10^3;

FrequencyUnit=1000000;
NumSteps=Length[AutlerTownesData12Mar20S5];
MatrixForm[Ψ0={1,0,0}];
inits={LogicalExpand[Ψ[0]==Ψ0][[1]],LogicalExpand[Ψ[0]==Ψ0][[2]],LogicalExpand[Ψ[0]==Ψ0][[3]]};

RemainingAtoms=Table[0,{NumSteps+1},{2}];
For[i=1,i≤(NumSteps+1),i++,
Deltaf=FrequencyUnit*AutlerTownesData12Mar20S5[[i,1]];
RemainingAtoms[[i,1]]=Deltaf/FrequencyUnit;
sols=NDSolve[Union[eqs/.{delta-> 2*Pi*Deltaf, Delta->0,Γa->0, Γg->0,Γb->Gammab,Omega1->Omega1Max,Omega2->Omega2Max},inits],{a,b,g},{t,0,Duration},MaxSteps->100000][[1]];
RemainingAtoms[[i,2]]=((Abs[a[Duration]])^2)/.sols;
];
ListPlot[{RemainingAtoms,AutlerTownesData12Mar20S5},Joined->True,PlotRange->{{-0.5*DetuningRange/FrequencyUnit,0.5*DetuningRange/FrequencyUnit},{0,1}}]

PowerLaser1=7*10^-6;  (* W *)
Waist=100*10 ^-6;
Intensity1=N[2*PowerLaser1/(Pi*Waist^2)];
PowerLaser2=2.8*10^-3;  (* W *)
Waist=100*10 ^-6;
Intensity2=N[2*PowerLaser2/(Pi*Waist^2)];
Print["Omega1/Sqrt[Intensity1] = ",Omega1Max/Intensity1];
Print["Omega2/Sqrt[Intensity2] = ",Omega2Max/Intensity2];