Time reversal paper

home776="/home/qitlab/Paper/timereverse/mathematica_files/281213/";
home780="/home/qitlab/Paper/timereverse/mathematica_files/311213/";

Applying the cavity to photon pairs - Theory and data analysis

Define parameters

R1=0.9410;
R2=0.9984;

ωs=2 π ((3 (10^8))/(776 ));
ωi=2π ((3 (10^8))/780);

Δ=2.7;
τ=2 5.9;

Define the photon wavepacket and cavity transfer function

ψ=A(Sqrt[(2/τ)])(E^(-I (ωs ts  + ωi ti) +((ts-ti)/τ)))   UnitStep[ti-ts] ;

TF=(((Sqrt[R1])-(Sqrt[R2])(e^(i (((ω-ωp-δ)/Δ)))))/(1-(Sqrt[R1 R2]) (e^(i(((ω-ωp-δ)/Δ)))))). 
Near cavity resonance  (ω-ωp-δ) << Δ, series expansion of exponential upto first-order gives,

TF=(((Sqrt[R1])-(Sqrt[R2])(1+I ((ω- ωp-δ)/Δ)))/(1-(Sqrt[R1 R2]) (1+I ((ω-ωp-δ)/Δ))));

Inverse - Fourier transform TF to convolve with ψ ,

TIF= InverseFourierTransform[TF,ω,t,Assumptions->Element[{R1,R2,ω,ωp,δ,Δ,ts,ti,τ,A},Reals]&& 0< R1< 1 && 0< R2< 1 && Δ>0 && τ>0 ]//FullSimplify;

Results - 776 photon through cavity

xmin=292;xmax=452;xpmin=960;xpmax=1040;

Needs["ErrorBarPlots`"]

Δ = 0

file = "18"; x0=186;

expData=Import[home776<>"combined_"<>file<>".dat","Table"];
normcnts=(490.18/16);

acc=((Total[expData[[801;;1200]][[All,2]]])/400)//N;

expData1=Transpose[{expData[[All,1]]-x0,((expData[[All,2]])/normcnts),((Sqrt[expData[[All,2]]])/normcnts)}][[xmin;;xmax]];

dataplot=ErrorListPlot[expData1,PlotRange->All,PlotStyle->{RGBColor[1,0,0]}];

fps=(1/(Sqrt[2π]))Convolve[ψ,TIF/.{t->ts,δ->2 π 0,ωp->ωs},ts,tsp]//FullSimplify;

fitfun= (acc/normcnts)+((Abs[fps])^2)/.{tsp->0,A->14}//FullSimplify;

fitplotn=Plot[fitfun/.{ts->0},{ti,-40,0},PlotRange->All,PlotStyle->{Black,Thickness[0.0035],Dashed}];
fitplotp=Plot[fitfun/.{ts->0},{ti,1,40},PlotRange->All,PlotStyle->{Black,Thickness[0.0035],Dashed}];

plot1=Show[{dataplot, fitplotn,fitplotp},AspectRatio->1/3,AxesOrigin->{-40,0},PlotRange->{{-40,40},{0,35}},TicksStyle->Directive[FontFamily->"Helvetica",FontSize->9,PrivateFontOptions->{"OperatorSubstitution"->False}]]

Export["rising_27.dat",Transpose[{expData1[[All,1]],expData1[[All,2]],expData1[[All,3]],Table[Re[fitfun]//N,{ti,-40.5,39.5,0.5}]}]];

Δ = γ

file = "4"; x0=186;

expData=Import[home776<>"combined_"<>file<>".dat","Table"];
normcnts=(240.128/16);

acc=((Total[expData[[801;;1200]][[All,2]]])/400)//N;

expData1=Transpose[{expData[[All,1]]-x0,((expData[[All,2]])/normcnts),((Sqrt[expData[[All,2]]])/normcnts)}][[xmin;;xmax]];

dataplot=ErrorListPlot[expData1,PlotRange->All,PlotStyle->{RGBColor[1,0,0]}];

fps=(1/(Sqrt[2π]))Convolve[ψ,TIF/.{t->ts,δ->2 π 0.027,ωp->ωs},ts,tsp]//FullSimplify;

fitfun= (acc/normcnts)+((Abs[fps])^2)/.{tsp->0,A->15}//FullSimplify;

fitplot=Plot[fitfun/.{ts->0},{ti,-40,40},PlotRange->All,PlotStyle->{Black,Thickness[0.0035],Dashed}];

plot2=Show[{dataplot, fitplot},AspectRatio->1/3,AxesOrigin->{-40,0},PlotRange->{{-40,40},{0,32}},TicksStyle->Directive[FontFamily->"Helvetica",FontSize->9,PrivateFontOptions->{"OperatorSubstitution"->False}]]

Export["rising_27.dat",Transpose[{expData1[[All,1]],expData1[[All,2]],expData1[[All,3]],Table[Re[fitfun]//N,{ti,-40.5,39.5,0.5}]}]];

Δ >> γ

file = "0"; x0=187;

expData=Import[home776<>"combined_"<>file<>".dat","Table"];
normcnts=(480.285/16);

acc=((Total[expData[[801;;1200]][[All,2]]])/400)//N;

expData1=Transpose[{expData[[All,1]]-x0,((expData[[All,2]])/normcnts),((Sqrt[expData[[All,2]]])/normcnts)}][[xmin;;xmax]];

dataplot=ErrorListPlot[expData1,PlotRange->All,PlotStyle->{RGBColor[1,0,0]}];

fps=(1/(Sqrt[2π]))Convolve[ψ,TIF/.{t->ts,δ->2 π 0.120,ωp->ωs},ts,tsp]//FullSimplify;

fitfun= (acc/normcnts)+((Abs[fps])^2)/.{tsp->0,A->15}//FullSimplify;

fitplotn=Plot[fitfun/.{ts->0},{ti,-40,0},PlotRange->All,PlotStyle->{Black,Thickness[0.0035],Dashed}];
fitplotp=Plot[fitfun/.{ts->0},{ti,1,40},PlotRange->All,PlotStyle->{Black,Thickness[0.0035],Dashed}];

plot3=Show[{dataplot, fitplotn,fitplotp},AspectRatio->1/3,AxesOrigin->{-40,0},PlotRange->{{-40,40},{0,33}},TicksStyle->Directive[FontFamily->"Helvetica",FontSize->9,PrivateFontOptions->{"OperatorSubstitution"->False}]]

Export["rising_27.dat",Transpose[{expData1[[All,1]],expData1[[All,2]],expData1[[All,3]],Table[Re[fitfun]//N,{ti,-40.5,39.5,0.5}]}]];

Results - 780 photon through cavity

xmin=620;xmax=780;xpmin=1070;xpmax=1150;

Δ = 0

file = "6";x0=324.5;

expData=Import[home780<>"combined_"<>file<>".dat","Table"];
normcnts=(378.97/16);

acc=((Total[expData[[801;;1200]][[All,2]]])/400)//N;

expData1=Transpose[{expData[[All,1]]-x0,((expData[[All,2]])/normcnts),((Sqrt[expData[[All,2]]])/normcnts)}][[xmin;;xmax]];

dataplot=ErrorListPlot[expData1,PlotRange->All,PlotStyle->{RGBColor[1,0,0]}];

fps=(1/(Sqrt[2π]))Convolve[ψ,TIF/.{t->ti,δ->2 π 0,ωp->ωi},ti,tip]//FullSimplify;

fitfun= (acc/normcnts)+((Abs[fps])^2)/.{ts->0,A->13}//FullSimplify;

fitplotn=Plot[fitfun/.{ts->0},{tip,-20,0},PlotRange->All,PlotStyle->{Black,Thickness[0.0035],Dashed}];
fitplotp=Plot[fitfun/.{ts->0},{tip,0.5,60},PlotRange->All,PlotStyle->{Black,Thickness[0.0035],Dashed}];

plot1p=Show[{dataplot, fitplotn,fitplotp},AspectRatio->1/3,AxesOrigin->{-20,0},PlotRange->{{-20,60},{0,30}},TicksStyle->Directive[FontFamily->"Helvetica",FontSize->9,PrivateFontOptions->{"OperatorSubstitution"->False}]]

Export["decaying_0.dat",Transpose[{expData1[[All,1]],expData1[[All,2]],expData1[[All,3]],Table[Re[fitfun1]//N,{tip,-15.5,64.5,0.5}]}]];

Plots

Cavity in signal

GraphicsRow[{plot1,plot2,plot3}]

Cavity in idler

Show[plot1p]

Estimation of photon number in the cavity

Defining known parameters

η=0.0016; (* Losses per round trip *)
Δη=0.001; (* Error in l *)
trt=((2 0.055)/0.3); (* round trip time in ns. = 1/FSR *)

Import on - resonant ((Subscript[G, on])) and off - resonant ((Subscript[G, off])) coincidence data file

(* Cavity in 776 signal mode *)
offdata=Import[home776<>"offres_776.dat","Table"];
ondata=Transpose[{Import[home776<>"combined_18.dat","Table"][[All,1]]-186.5,Import[home776<>"combined_18.dat","Table"][[All,2]]}];

(* Cavity in 780 idler mode *)
offdatap=Import[home780<>"offres_780.dat","Table"];
ondatap=Transpose[{Import[home780<>"combined_6.dat","Table"][[All,1]]-323,Import[home780<>"combined_6.dat","Table"][[All,2]]}];

Determine accidental coincidences on all data files

offacc=((Total[offdata[[401;;800,2]]])/400)//N;
onacc=((Total[ondata[[601;;1000,2]]])/400)//N;

offaccp=((Total[offdatap[[401;;800,2]]])/400)//N;
onaccp=((Total[ondatap[[301;;600,2]]])/300)//N;

Selecting  ± 50 ns of coincidences from the trigger, subtract accidental

offdata1=Sort[Transpose[{-1 offdata[[868;;1068,1]],offdata[[868;;1068,2]]-offacc}]];
ondata1=Sort[Transpose[{-1 ondata[[275;;475,1]],ondata[[275;;475,2]]-onacc}]];

offdatap1=Sort[Transpose[{offdatap[[1000;;1200,1]],offdatap[[1000;;1200,2]]-offaccp}]];
ondatap1=Sort[Transpose[{ondatap[[571;;771,1]],ondatap[[571;;771,2]]-onaccp}]];

GraphicsRow[{ListLinePlot[{offdata1,ondata1},PlotRange->{{-50,50},{-5,1100}}],ListLinePlot[{offdatap1,ondatap1},PlotRange->All]}]

Compute total pairs in the off - resonant case used for normalising to 1

pairs=Total[offdata1[[All,2]]];

pairsp=Total[offdatap1[[All,2]]];

Compute the photon number n(t). Loss η is inculded to the photon number in the cavity
	n(t)=(∫_(-∞)^t)((((Subscript[G, off])(s) -(Subscript[G, on])(s))/pairs) - n(s) (η/trt)) ds 
	=>((dn)/(dt))= (((Subscript[G, off])(t) -(Subscript[G, on])(t))/pairs) - n(t) (η/trt)
	n(t)= (e^(- η (t/trt))) (∫_(-∞)^t)(((Subscript[G, off])(s) -(Subscript[G, on])(s))/pairs)(e^(η (s/trt))) ds

n=offdata1;error=Table[0,{i,1,Length[offdata1]}];

For[i=1;sum=0;err=0,i≤201,i++,
sum+=(((offdata1[[i,2]]-ondata1[[i,2]]))/pairs) (E^((η offdata1[[i,1]])/trt));
n[[i,2]]=(E^(- ((η offdata1[[i,1]])/trt))) sum;
err+= ((((((offdata1[[i,2]]-ondata1[[i,2]]))/pairs) (E^((η offdata1[[i,1]])/trt))  ((offdata1[[i,1]])/trt) Δη))^2)+(((offdata1[[i,2]]+ondata1[[i,2]])/(pairs^2))(E^((2 η offdata1[[i,1]])/trt)));
error[[i]]=(E^(- ((η offdata1[[i,1]])/trt))) √err;
]

np=offdatap1;errorp=Table[0,{i,1,Length[offdatap1]}];

For[i=1;sum=0;err=0,i≤201,i++,
sum+=(((offdatap1[[i,2]]-ondatap1[[i,2]]))/pairsp) (E^((η offdatap1[[i,1]])/trt));
np[[i,2]]=(E^(- ((η offdatap1[[i,1]])/trt))) sum;
err+= ((((((offdatap1[[i,2]]-ondatap1[[i,2]]))/pairsp) (E^((η offdatap1[[i,1]])/trt))  ((offdatap1[[i,1]])/trt) Δη))^2)+(((offdatap1[[i,2]]+ondatap1[[i,2]])/(pairsp^2))(E^((2 η offdatap1[[i,1]])/trt)));
errorp[[i]]=(E^(- ((η offdatap1[[i,1]])/trt))) √err;
]

nerr=Transpose[{n,Map[ErrorBar,error]}];

nperr=Transpose[{np,Map[ErrorBar,errorp]}];

Plot of mean photon number

Needs["ErrorBarPlots`"]

GraphicsRow[{ErrorListPlot[nperr,PlotRange->{{-20,60},{-0.1,1}},PlotStyle->Red],ErrorListPlot[nerr,PlotRange->{{-40,40},{-0.1,1}}]}]

export the plot points to data file for gnuplot

Export["photon_number_signal.dat",Transpose[{n[[All,1]],n[[All,2]],error}]]

Export["photon_number_idler.dat",Transpose[{np[[All,1]],np[[All,2]],errorp}]]