Pulse propagation and compression (C. Durfee)

initialization:

constants and unit definitions: internally, all numbers in SI, except pressure in atm
angles in Mathematica are internally in radians, my variables are in degrees

ps=(10^(-12));fs=(10^(-15));
mm=(10^(-3));cm=(10^(-2));cm3=(10^6);um=(10^(-6));nm=(10^(-9));angs=(10^(-10));pm=(10^(-12));
mJ=(10^(-3));uJ=(10^(-6));W=1; V=1;
c=2.9979 10^8;                          (*speed of light,in m/s*)
eps0=8.8542 (10^(-12));              (*permittivity of free space*)
re=2.82 (10^(-15));	              (*m classical electron radius*)
e=1.6 (10^(-19));	                   (*Coul electron charge*)
Natm=760 3.3 (10^16)/(cm^3); (*no.particles in 1 atm*)
Torr=1/760;
eV=1.6 (10^(-19));
hbar=(10^(-34));

refractive index dispersion definitions

define indices

refractive index for BK7: nbk7

nbk7[lambda_]:=√(2.2718929-(1.0108077 (lambda/um)^2)/10^2+1.0592509/(10^2 (lambda/um)^2)+2.0816965/(10^4 (lambda/um)^4)-7.6472538/(10^6 (lambda/um)^6)+4.9240991/(10^7 (lambda/um)^8));

Define the phase and group delay terms

phase shift due to materials

improve: include second set for tilted materials

Clear[nmat,lmat,λ];
nmat[λ_]=nbk7[λ];
ϕmat[λ_,lmat_]:=(2 π/λ)lmat nmat[λ];
ϕ1mat[λ_,lmat_]=(-((λ^2)/(2 π c))) D[ϕmat[λ,lmat],λ];
ϕ2mat[λ_,lmat_]=(-((λ^2)/(2 π c))) D[ϕ1mat[λ,lmat],λ];
ϕ3mat[λ_,lmat_]=(-((λ^2)/(2 π c))) D[ϕ2mat[λ,lmat],λ];
ϕ4mat[λ_,lmat_]=(-((λ^2)/(2 π c))) D[ϕ3mat[λ,lmat],λ];
ϕ5mat[λ_,lmat_]=(-((λ^2)/(2 π c))) D[ϕ4mat[λ,lmat],λ];

taylor expanded phase

0th and 1st order terms are not included
1/n! terms should be included in coeff

Clear[kk2,kk3,kk4,kk5,Δω];
Δω[λ_]:=2 π c(-1/lambda0+1/λ);
ϕTaylor[λ_,kk2_,kk3_,kk4_,kk5_,kk6_]:=(1/(2!))kk2 ((Δω[λ])^2)+(1/(3!))kk3 ((Δω[λ])^3)+(1/(4!))kk4 ((Δω[λ])^4)+(1/(5!))kk5 ((Δω[λ])^5)+(1/(6!))kk6 ((Δω[λ])^6);
ϕ1Taylor[λ_,kk2_,kk3_,kk4_,kk5_,kk6_]=-((λ^2)/(2 π c)) D[ϕTaylor[λ,kk2,kk3,kk4,kk5,kk6],λ];
ϕ2Taylor[λ_,kk2_,kk3_,kk4_,kk5_,kk6_]=-((λ^2)/(2 π c)) D[ϕ1Taylor[λ,kk2,kk3,kk4,kk5,kk6],λ];
ϕ3Taylor[λ_,kk2_,kk3_,kk4_,kk5_,kk6_]=-((λ^2)/(2 π c)) D[ϕ2Taylor[λ,kk2,kk3,kk4,kk5,kk6],λ];
ϕ4Taylor[λ_,kk2_,kk3_,kk4_,kk5_,kk6_]=-((λ^2)/(2 π c)) D[ϕ3Taylor[λ,kk2,kk3,kk4,kk5,kk6],λ];

grating pair phase shift

This equation is for single pass configuration. For a grating pair separated by L, the beam normally goes through 2 L

The original formula can be found in:
  C.H.Brito Cruz et al,Optics Letters,v13,p123(1988)
Note:  compressor and stretcher have opposite sign in phase
  
Diffraction equation for a grating:

θdif[θinc_,λ_,d_]:=ArcSin[-(λ/d)+Sin[θinc]];

γLittrow: Littrow angle is the incident angle where the -1 order comes back at the input

γLittrow[λ_,d_]:=ArcSin[(λ/(2 d))];

ϕgrastr[λ_,lstr_,gstr_,dstr_]:=-((2 π)/λ) lstr (Sqrt[1-((((λ/dstr)-Sin[gstr]))^2)]);
ϕgracom[λ_,lcom_,gcom_,dcom_]=-ϕgrastr[λ,lcom,gcom,dcom];

calculate derivatives

ϕ1grastr[λ_,lstr_,gstr_,dstr_]=(-((λ^2)/(2 π c))) D[ϕgrastr[λ,lstr,gstr,dstr],λ];ϕ2grastr[λ_,lstr_,gstr_,dstr_]=(-((λ^2)/(2 π c))) D[ϕ1grastr[λ,lstr,gstr,dstr],λ];
ϕ3grastr[λ_,lstr_,gstr_,dstr_]=(-((λ^2)/(2 π c))) D[ϕ2grastr[λ,lstr,gstr,dstr],λ];
ϕ4grastr[λ_,lstr_,gstr_,dstr_]=(-((λ^2)/(2 π c))) D[ϕ3grastr[λ,lstr,gstr,dstr],λ];

ϕ1gracom[λ_,lcom_,gcom_,dcom_]=(-((λ^2)/(2 π c))) D[ϕgracom[λ,lcom,gcom,dcom],λ];
ϕ2gracom[λ_,lcom_,gcom_,dcom_]=(-((λ^2)/(2 π c))) D[ϕ1gracom[λ,lcom,gcom,dcom],λ];
ϕ3gracom[λ_,lcom_,gcom_,dcom_]=(-((λ^2)/(2 π c))) D[ϕ2gracom[λ,lcom,gcom,dcom],λ];
ϕ4gracom[λ_,lcom_,gcom_,dcom_]=(-((λ^2)/(2 π c))) D[ϕ3gracom[λ,lcom,gcom,dcom],λ];

optimization of compressor settings

1. glass and grating compressor

a. Make a table for the dispersion orders for 1 cm BK7, a grating-based stretcher and compressor (1cm separation, Littrow incidence angle, groove density 1200 gr/mm). The center wavelength should be 800nm.
The dispersion for the nth order is best expressed in units of (fs^n)

λ0=800 nm;
lenBK7=0.1cm;
dstr=dcom=1/1200. mm; 
γstr=γcom=γLittrow[λ0,dcom];
lstr=lcom=1cm;

bk7Orders={ϕ2mat[λ0,lenBK7]/(fs^2),ϕ3mat[λ0,lenBK7]/(fs^3),ϕ4mat[λ0,lenBK7]/(fs^4)};
strOrders={ϕ2grastr[λ0,lstr,γstr,dstr]/(fs^2),ϕ3grastr[λ0,lstr,γstr,dstr]/(fs^3),ϕ4grastr[λ0,lstr,γstr,dstr]/(fs^4)};
comOrders={ϕ2gracom[λ0,lcom,γcom,dcom]/(fs^2),ϕ3gracom[λ0,lcom,γcom,dcom]/(fs^3),ϕ4gracom[λ0,lcom,γcom,dcom]/(fs^4)};

Print["orders of dispersion in (fs^n)"];
Print["\t BK7 \t stretcher \t compressor"];
Print[TableForm[Transpose[{{"ϕ2","ϕ3","ϕ4"},bk7Orders,strOrders,comOrders}]]];

b. From these numbers, calculate the grating separation required to compensate second-order phase for 100cm of BK7 glass. 
- first make a plot of the total ϕ2 as a function of grating separation to see where the zero is
- then use FindRoot[ ] to get the result

λ0=800 nm;
lenBK7=100cm;
dcom=1/1200. mm; 
γcom=γLittrow[λ0,dcom];
Print["incident angle: ",γcom/Degree//N];

Print["total second-order material phase: ",ϕ2mat[λ0,lenBK7]/(fs^2)," (fs^2)"]
Print["ϕ2: material and compressor "];
Plot[{ϕ2mat[λ0,lenBK7]/(fs^2)+ϕ2gracom[λ0,lcom mm,γcom,dcom]/(fs^2)},{lcom,0,2},PlotStyle->{Red,Blue}]

ans=FindRoot[ϕ2mat[λ0,lenBK7]/(fs^2)+ϕ2gracom[λ0,lcomTst mm,γcom,dcom]/(fs^2)==0,{lcomTst,0.5}];
lcomOpt=lcomTst mm/.ans[[1]];
Print["Optimum grating separation: ",lcomOpt/mm," mm"];

c. For this compressor separation, plot the net group delay vs wavelength over a 50nm range around the center wavelength. Compare the shape of the group delay to what you expect. 
- subtract the group delay at (Subscript[λ, 0]) since that overall time is not meaningful

ΔλIn=50nm;
gdelCtr=ϕ1mat[λ0,lenBK7]+ϕ1gracom[λ0,lcomOpt,γcom,dcom];
Plot[(ϕ1mat[λ nm,lenBK7]+ϕ1gracom[λ nm,lcomOpt,γcom,dcom]-gdelCtr)/fs,{λ,(λ0-ΔλIn/2)/nm,(λ0+ΔλIn/2)/nm}]

Clear[gdelCtr,lcom];
ΔλIn=50nm;
gdelCtr[lcom_]:=ϕ1mat[λ0,lenBK7]+ϕ1gracom[λ0,lcom,γcom,dcom];
Manipulate[Plot[(ϕ1mat[λ nm,lenBK7]+ϕ1gracom[λ nm,lcom mm,γcom,dcom]-gdelCtr[lcom mm])/fs,{λ,(λ0-ΔλIn/2)/nm,(λ0+ΔλIn/2)/nm},AxesLabel->{"λ (nm)","τg (fs)"}],
{lcom,0,50}]

We expect this to be parabolic, since we’ve zeroed out the linear chirp. It does have curvature, but isn’t symmetric because of higher-order phase.

d. Estimate the maximum pulse duration that can be reasonably well compressed: The group delay difference across the bandwidth shouldn’t be more than half the pulse duration. As a way to scale the bandwidth, use the approximation that a 10 fs duration has about 90nm of bandwidth.

Δλ10=90nm;
τ10=10.fs;
ΔλIn=10nm;
Print["duration estimate: ",τ10 Δλ10/ΔλIn /fs," fs"];
gdelCtr=ϕ1mat[λ0,lenBK7]+ϕ1gracom[λ0,lcomOpt,γcom,dcom];
Plot[(ϕ1mat[λ nm,lenBK7]+ϕ1gracom[λ nm,lcomOpt,γcom,dcom]-gdelCtr)/fs,{λ,(λ0-ΔλIn/2)/nm,(λ0+ΔλIn/2)/nm}]

estimate of around 90 fs pulse duration. This would be broadened out to roughly 120fs.

Part 2: Stretching, glass propagation, compression. 
With a grating-based stretcher and compressor, we can make the separation and incident angles of the gratings different, giving you two degrees of freedom to compensate for both 2nd and 3rd order dispersion.

a. Estimate the amount of stretch you get from the 50cm BK7 glass and for a stretcher with an effective grating separation of 20cm and Littrow incident angle. Assume the input bandwidth is 50nm. You can either estimate this by calculating it directly from the group variation across the bandwidth, or from the second order phase of the taylor expansion.

λ0=800. nm;
ω0=2 π c/λ0;
ΔλIn=50nm;
ΔωIn=ω0 ΔλIn/λ0;
lenBK7=50cm;
dstr=dcom=1/1200 mm; 
γstr=γLittrow[λ0,dcom];
γstr=32Degree;

lstr=lcom=20cm;

Print["stretch from glass (direct): ",ϕ1mat[λ0-ΔλIn/2,lenBK7]/ps-ϕ1mat[λ0+ΔλIn/2,lenBK7]/ps," ps"];
Print["stretch from glass (Taylor): ",ϕ2mat[λ0,lenBK7]ΔωIn/ps," ps"];

Print["grating stretch (direct): ",ϕ1grastr[λ0-ΔλIn/2,lstr,γstr,dstr]/ps-ϕ1grastr[λ0+ΔλIn/2,lstr,γstr,dstr]/ps," ps"];
Print["grating stretch (Taylor): ",ϕ2grastr[λ0,lstr,γstr,dstr]ΔωIn/ps," ps"];

b. With a fixed stretcher grating separation and incident angle, you have two degrees of freedom for the compressor. The following is a procedure for optimizing these settings. 
- Calculate the ratio of ϕ3/ϕ2 for the glass, for the stretcher and for the combination of the stretcher and the glass. (Notice that the first two ratios are independent of the glass thickness or grating separation). 
- For the stretcher and compressor combination: find the grating incident angle for the compressor that has the same absolute value of ϕ3/ϕ2. 
- For this optimum incident angle, find the compressor grating separation that compresses the stretched pulse.

ϕ3ϕ2mat=((ϕ3mat[λ0,lenBK7]/(fs^3))/(ϕ2mat[λ0,lenBK7]/(fs^2)));
ϕ3ϕ2str=((ϕ3grastr[λ0,lstr,γstr,dstr]/(fs^3))/(ϕ2grastr[λ0,lstr,γstr,dstr]/(fs^2)));
ϕ3ϕ2both=((ϕ3grastr[λ0,lstr,γstr,dstr]/(fs^3)+ϕ3mat[λ0,lenBK7]/(fs^3))/(ϕ2grastr[λ0,lstr,γstr,dstr]/(fs^2)+ϕ2mat[λ0,lenBK7]/(fs^2)));
ϕ3ϕ2com[γcom_]:=((ϕ3gracom[λ0,lcom,γcom,dcom]/(fs^3))/(ϕ2grastr[λ0,lcom,γcom,dcom]/(fs^2)));

Print["ratio of ϕ3/ϕ2 for material: ",ϕ3ϕ2mat," /fs"];
Print["ratio of ϕ3/ϕ2 for stretcher: ",ϕ3ϕ2str," /fs"];
Print["ratio of ϕ3/ϕ2 for both: ",ϕ3ϕ2both," /fs"];
Plot[ϕ3ϕ2com[γcom Degree]+ϕ3ϕ2both,{γcom,30,40},PlotStyle->{Red,Blue}]

first find the AOI that gives the right ratio of ϕ3/ϕ2

ans=FindRoot[ϕ3ϕ2com[γcomTst Degree]+ϕ3ϕ2both==0,{γcomTst,33}];
γcomOpt=γcomTst Degree/.ans[[1]];
Print["Optimum grating angle: ",γcomOpt/Degree," deg"];

then at this angle, find the compressor separation that eliminates ϕ2

Plot[(ϕ2grastr[λ0,lstr,γstr,dstr]+ϕ2mat[λ0,lenBK7]+ϕ2gracom[λ0,lcom cm,γcomOpt,dcom])/(fs^2),{lcom,15,25}]

ans=FindRoot[(ϕ2grastr[λ0,lstr,γstr,dstr]+ϕ2mat[λ0,lenBK7]+ϕ2gracom[λ0,lcomTst cm,γcomOpt,dcom])/(fs^2)==0,{lcomTst,30}];
lcomOpt=lcomTst cm/.ans[[1]];
Print["optimum compressor separation: ",lcomOpt/cm," cm"];

c. Plot the net group delay vs wavelength over a 50nm range around the center wavelength. Compare the shape of the group delay to what you expect.

gdelCtr=ϕ1grastr[λ0,lstr,γstr,dstr]+ϕ1mat[λ0,lenBK7]+ϕ1gracom[λ0,lcomOpt,γcomOpt,dcom];
Plot[(ϕ1grastr[λ nm,lstr,γstr,dstr]+ϕ1mat[λ nm,lenBK7]+ϕ1gracom[λ nm,lcomOpt,γcomOpt,dcom]-gdelCtr)/fs,{λ,750,850}]

Clear[gdelCtr];
gdelCtr[lcom_,γcom_]:=ϕ1grastr[λ0,lstr,γstr,dstr]+ϕ1mat[λ0,lenBK7]+ϕ1gracom[λ0,lcom,γcom,dcom];
Manipulate[
Plot[(ϕ1grastr[λ nm,lstr,γstr,dstr]+ϕ1mat[λ nm,lenBK7]+ϕ1gracom[λ nm,lcom cm,γcom Degree,dcom]-gdelCtr[lcom cm,γcom Degree])/fs,{λ,750,850}],
{lcom,15,25},{γcom,29,35}]

d. Estimate the maximum pulse duration that can be reasonably well compressed: The group delay difference across the bandwidth shouldn’t be more than half the pulse duration. As a way to scale the bandwidth, use the approximation that a 10 fs duration has about 90nm of bandwidth.

Δλ10=90nm;
τ10=10.fs;
ΔλIn=45nm;
Print["duration estimate: ",τ10 Δλ10/ΔλIn /fs," fs"];
ϕ1system[λ_]:=ϕ1grastr[λ,lstr,γstr,dstr]+ϕ1mat[λ,lenBK7]+ϕ1gracom[λ,lcomOpt,γcomOpt,dcom];
gdelCtr=ϕ1system[λ0];
Plot[(ϕ1system[λ nm]-gdelCtr)/fs,{λ,(λ0-ΔλIn/2)/nm,(λ0+ΔλIn/2)/nm}]

Define FFT functions for temporal domain

FFT[fLst_,npts_,dt_]:=(Sqrt[npts (dt^2)])RotateLeft[Fourier[RotateRight[fLst,npts/2]],npts/2];
InverseFFT[FLst_,npts_,dt_]:=(Sqrt[(1/(npts (dt^2)))])RotateLeft[InverseFourier[RotateRight[FLst,npts/2]],npts/2];

plot the output pulse shape in the time domain

set up grids

frequency variable: ω

ntpts=(2^12);
dt = 5.fs;
λ0=800.nm;
ω0=2π c/λ0;

tSpan=ntpts dt;
tLst=Table[i dt,{i,-ntpts/2,ntpts/2-1}];
dω=2π/(ntpts dt);
ωLst=Table[ω0+iω dω,{iω,-ntpts/2,ntpts/2-1}];
ωSpan=ntpts dω; 

λLst=2π c/ωLst;
dλ=λ0 dω/ω0;
Print["tSpan: ",tSpan/fs," fs, dt: ",dt/fs," fs"];
Print["ωSpan: ",ωSpan fs," rad/fs, dω: ",dω fs," /fs"];
Print["λMin: ",λLst[[ntpts]]/nm," nm. λMax: ",λLst[[1]]/nm," nm, dλ: ",dλ/nm," nm"];

λGratingMax=λ nm/.FindRoot[Sin[θdif[γstr,λ nm,dstr]]==-1,{λ,1000}];
Print["maximum wavelength that can be diffracted from the stretcher grating: ",λGratingMax/nm," nm"];
λGratingMax=λ nm/.FindRoot[Sin[θdif[γcomOpt,λ nm,dcom]]==-1,{λ,1000}];
Print["maximum wavelength that can be diffracted from the compressor grating: ",λGratingMax/nm," nm"];

define pulse shape

τIn and Δω are 1/e2 half widths in intensity

τIn=35fs;
ΔωIn=2/τIn;

AtInLst=Exp[-(tLst^2)/((τIn)^2)];
AωInLst=FFT[AtInLst,ntpts,dt];

Print["time domain"];
tRange=4 τIn;
ListPlot[Transpose[{tLst/fs,((Abs[AtInLst])^2)}],Joined->True,PlotRange->{{-tRange/fs,tRange/fs},All},AxesLabel->{"t (fs)","intensity"}]

Print["spectral domain"];
ωRange=3ΔωIn;
ListPlot[Transpose[{ωLst fs,((Abs[AωInLst])^2)}],Joined->True,PlotRange->{{(ω0-ωRange)fs,(ω0+ωRange)fs},All},AxesLabel->{"ω (rad/fs)","intensity"}]

test sample phase profiles

test with Taylor phase

lambda0=800. nm;

kk2=0;
kk3=0;
kk4=1000000(fs^4);
kk5=kk6=0;
ϕsystem[λ_]:=ϕTaylor[λ,kk2,kk3,kk4,kk5,kk6];

phaseTstLst=ϕsystem[λLst]-(ϕsystem[λ0]);

AtLst=InverseFFT[AωInLst Exp[I phaseTstLst],ntpts,dt];
Print["time domain"];
ListPlot[Transpose[{tLst/fs,Abs[AtLst]}],Joined->True,PlotRange->{{-tRange/fs,tRange/fs},All}]

test with grating settings from above

ϕsystem[λ_]:=ϕgrastr[λ,lstr,γstr,dstr]+ϕmat[λ,lenBK7]+ϕgracom[λ,lcomOpt,γcomOpt,dcom];

phaseTstLst=Re[ϕsystem[λLst]-(ϕsystem[λ0]+ϕ1system[λ0](ωLst-ω0))];

ListLinePlot[Transpose[{λLst/nm,phaseTstLst}]]

ListLinePlot[Transpose[{λLst/nm,Abs[AωInLst Exp[I phaseTstLst]]}],PlotRange->All]

AtLst=InverseFFT[AωInLst Exp[I phaseTstLst],ntpts,dt];
Print["time domain"];
ListPlot[Transpose[{tLst/fs,Abs[AtLst]}],Joined->True,PlotRange->{{-tRange/fs,tRange/fs},All}]