phase match code

define constants

constants

internal units are um=1,fs=1

sec=(10^15);fs=(10^(-15))sec;ns=(10^(-9))sec;ps=(10^(-12))sec; (* internal fs=1 *) 
meter=(10^6); nm=(10^(-9))meter; mm=(10^(-3))meter; um=(10^(-6))meter; (* internal um=1; *) 
c=3 (10^8)meter/sec;

define dispersion for uniaxial crystals (KDP, BBO)

define a general form for the Sellmeier dispersion equation

These empirical equations are not always formatted this way.
The dispersion equations assume the wavelength is in microns so we convert from SI.

indexSellmeier[λ_,aA_,bB_,cC_,dD_,eE_]:=(Sqrt[aA+((bB (((λ um))^2))/((((λ um))^2)-cC))+(dD/((((λ um))^2)-eE))]);

define the angle-dependent index for uniaxial crystals

Clear[no,ne,neθ];
neθ[λ_,θ_]=(Sqrt[1/((((Cos[θ])^2)/((no[λ])^2))+(((Sin[θ])^2)/((ne[λ])^2)))]);

KD*P: deuterated KDP

aoDKDP=2.2409;
boDKDP=2.2470;
coDKDP=126.9205;
doDKDP=0.0097;
eoDKDP=0.0156;
aeDKDP=2.1260;
beDKDP=0.7844;
ceDKDP=123.4032;
deDKDP=0.0086;
eeDKDP=0.0120;

no[λ_]=indexSellmeier[λ,aoDKDP,boDKDP,coDKDP,doDKDP,eoDKDP];
ne[λ_]=indexSellmeier[λ,aeDKDP,beDKDP,ceDKDP,deDKDP,eeDKDP];

KDP

aoKDP=2.259276;
boKDP=13.00522;
coKDP=400;
doKDP=0.01008956;
eoKDP=0.012942625;
aeKDP=2.132668;
beKDP=3.2279924;
ceKDP=400;
deKDP=0.008637494;
eeKDP=0.012281043;

no[λ_]=indexSellmeier[λ,aoKDP,boKDP,coKDP,doKDP,eoKDP];
ne[λ_]=indexSellmeier[λ,aeKDP,beKDP,ceKDP,deKDP,eeKDP];

BBO Sellmeier equations (from Casix, λ in microns)

negative uniaxial: ne < no

no[λ_]=Sqrt[2.7359+0.01878/((((λ um))^2)-0.01822)-0.01354 (((λ um))^2)];
ne[λ_]=Sqrt[2.3753+0.01224/((((λ um))^2)-0.01667)-0.01516 (((λ um))^2)];

test plot to make sure Type I phase matching is possible

λTst=1.064um;
Print["no at fund: ",no[λTst]];
Print["ne range: ",ne[λTst/2]," to ",no[λTst/2]];
Plot[{no[λ um],ne[λ um]},{λ,.4,1.064},PlotStyle->{Red,Blue,Gray},GridLines->{{λTst/2/um,λTst/um},{no[λTst]}},Frame->True,FrameLabel->{"λ (μm)","ref index"}]

Plot[neθ[λTst/2 ,θ Degree],{θ,0,90},GridLines->{None,{no[λTst/2],ne[λTst/2]}},Frame->True,FrameLabel->{"θ (deg)","ref index"}]

define dispersion for biaxial crystals (LBO)

LBO Sellmeier equations (from Laser Components datasheet, lambda in microns, theta in degrees)

these are list in order of increasing index value
for type I, call the higher index = no, lower ne

nx[λ_]=√(2.454140+(0.011249/((((λ um))^2)-0.011350))-0.014591 (((λ um))^2)-6.60  (10^(-5))(((λ um))^4));
ny[λ_]=(Sqrt[2.539070+(0.012711/((((λ um))^2)-0.012523))-0.018540 (((λ um))^2)+2.00 (10^(-4))(((λ um))^4)]);
nz[λ_]=(Sqrt[2.586179+(0.013099/((((λ um))^2)-0.011893))-0.017968 (((λ um))^2)-2.26 (10^(-4))(((λ um))^4)]);

Plot[{nx[λ um],ny[λ um],nz[λ um]},{λ,.4,1.064},PlotStyle->{Red,Blue,Gray},Frame->True,FrameLabel->{"λ (μm)","ref index"}]

Set up angle-tuned index for Type I
By analogy to type I in negative uniaxial crystals, we place the fundamental polarization along one of the higher index axes (y or z), then tune the harmonic as a combination of the low axis x and the other. 
Version 1: (Subscript[E, 1]) along y, (Subscript[E, 2]) in xz plane
Version 2: (Subscript[E, 1]) along z, (Subscript[E, 2]) in xy plane 
   this configuration is labeled θ=90, φ=11.3, in spherical coords, this could be the direction of k: θ=90 puts k perp to z, in the xy plane. φ=11.3 would be k measured from the x-axis.

Clear[neθ,no,ne1,ne2];

no[λ_]=nz[λ];
ne1[λ_]=nx[λ];
ne2[λ_]=ny[λ];
neθ[λ_,θ_]=(Sqrt[1/((((Cos[θ])^2)/((ne2[λ])^2))+(((Sin[θ])^2)/((ne1[λ])^2)))]);

λTst=1.064 um;
Print["no at fund: ",no[λTst]];
Print["ne range: ",ne1[λTst/2]," to ",ne2[λTst/2]];

Plot[neθ[λTst/2,θ Degree],{θ,0,90},GridLines->{None,{no[λTst]}},Frame->True,FrameLabel->{"θ (deg)","ref index"}]

phase matching angles and group velocity walkoff

Define k, ϕ and group delay ϕ1

keθ[λ_,θ_]:=((2 π)/λ)  neθ[λ,θ] ;
ko[λ_]:=((2 π)/λ)  no[λ] ;

The group delay is dphi=d(phi)/dw. (Type I)
as defined here, Δτ > 0 if the SH takes longer to leave the crystal than the fundamental

ϕEwave[λ_,θ_,z_]:=keθ[λ,θ] z;
ϕOwave[λ_,z_]:=ko[λ] z;
ϕ1Ewave[λ_,θ_,z_]=-((λ^2)/(2 π c)) D[ϕEwave[λ,θ,z],λ];
ϕ1Owave[λ_,z_]=-((λ^2)/(2 π c)) D[ϕOwave[λ,z],λ];
ΔτgroupDelay[λ_,θ_,z_]:=ϕ1Ewave[λ/2,θ,z]-ϕ1Owave[λ,z];

Type I doubling

In Type I, the fundamental travels as an ordinary wave and the harmonic as an extraordinary wave, so the phases (and group delay) of each should be treated accordingly.

ΔkTypeIdoubling[λ_,θ_]=2ko[λ]-keθ[λ/2,θ ];
θPMtypeI[λ_,θTst_]:=θ/.FindRoot[ΔkTypeIdoubling[λ,θ]==0,{θ,θTst}][[1]];

Find phase-matching angle for negative uniaxial crystal, where no(fundamental)=ne(SH).
when this angle is zero, then the optic axis is perp to the plane of the surface (SH sees (Subscript[n, o]) only, which is higher than (Subscript[n, o]) of the fund) 
for Type I doubling KDP angle is 43.8 degrees for 800nm to 400nm
LBO: 2-1 49.5deg, 3-1 36.9deg, 3-2 no soln

λFund=1.064 um;θTst=30 Degree;
θpmI=θPMtypeI[λFund,θTst];
Plot[{no[λFund],neθ[((λFund)/2),θ Degree]},{θ,0,90},
PlotRange->All,PlotStyle->{Red,Green},GridLines->{{θpmI/Degree},None},Frame->True,FrameLabel->{"θ (deg)","ref index"}]
Print["phase matching angle: ",θpmI/Degree," deg"];

show variation in phase matching angle with wavelength

Plot[θPMtypeI[λFund um,θTst]/Degree,{λFund,.9,1.1},Frame->True,FrameLabel->{"λ input (μm)","phase match angle (deg)"}]

spatial walkoff angle
check: for type 1 doubling of 1064nm should be 2.7Deg for BBO, 28mrad for KDP 
deviation from this in BBO is probably from using slightly different index defs.

αWalkOff[λ_,θ_]:=ArcCos[((((ne[λ])^2)((Cos[θ])^2)+((no[λ])^2)((Sin[θ])^2))/(Sqrt[((ne[λ])^4)((Cos[θ])^2)+((no[λ])^4)((Sin[θ])^2)]))];
Print["walk-off angle: ",αWalkOff[λFund/2,θpmI]," deg = ",αWalkOff[λFund/2,θpmI] 1000," mrad"];

temporal walkoff of pulse envelopes
for pulsed input the fundamental travels at a different group velocity than the harmonic

cryLen=2mm;
Print[" for a crystal length of: ",cryLen/mm," mm, the relative group delay is: ",ΔτgroupDelay[λFund,θpmI,cryLen]/fs," fs"];

Depending on the crystal, there can be an angle where the group velocities match, but this is generally different than where the phase velocities match.

Plot[ΔτgroupDelay[800 nm,θ Degree,cryLen]/fs,{θ,0,90},Frame->True,FrameLabel->{"θ (deg)","group delay Δτ (fs)"}]

type II doubling

for type II, the fundamental is linearly polarized at 45deg between no and ne(θ) axes to give equal amplitude between the two components. The SH is along ne(θ). 
Δk=((2 π)/(Subscript[λ, 1]))(Subscript[n, o])((Subscript[λ, 1]))+((2 π)/(Subscript[λ, 1]))(Subscript[n, e])((Subscript[λ, 1]),θ)-((2 π)/(Subscript[λ, 2]))(Subscript[n, e])((Subscript[λ, 2]),θ)
Δk=((2 π)/(Subscript[λ, 1]))((Subscript[n, o])((Subscript[λ, 1]))+(Subscript[n, e])((Subscript[λ, 1]),θ)-2(Subscript[n, e])((Subscript[λ, 2]),θ))
the minimum wavelength that can be phase matched for doubling is longer for type II

ΔkTypeIIdoubling[λ_,θ_]=keθ[λ,θ]+ko[λ]-keθ[λ/2,θ];
θPMtypeII[λ_,θTst_]:=θ/.FindRoot[ΔkTypeIIdoubling[λ,θ]==0,{θ,θTst}][[1]];

λFund=1.064 um;θTst=36 Degree;
θpmII=θPMtypeII[λFund,θTst];
Plot[{
no[λFund]+neθ[λFund,θ Degree],
2 neθ[((λFund)/2),θ Degree]
},{θ,0,90},GridLines->{{θpmII/Degree},None},PlotRange->All,PlotStyle->{Red,Green},Frame->True,FrameLabel->{"θ (deg)","ref index"}]
Print["optimum phase matching angle (Type II) ",θpmII/Degree," degrees"];

Type II temporal walkoff

cryLen=1mm;
Print["For a crystal length of: ",cryLen/mm," mm"];
Print["the relative group delay is: ",
(ϕ1Ewave[λFund/2,θpmII,cryLen]-ϕ1Owave[λFund,cryLen])/fs," fs between the E wave harmonic and the O wave fundamental"];
Print["the relative group delay is: ",
(ϕ1Ewave[λFund/2,θpmII,cryLen]-ϕ1Ewave[λFund,θpmII,cryLen])/fs," fs between the E wave fundamental and harmonic"];

Acceptance angle and bandwidth

tolerance to angle detuning (e.g. for diverging beam)

pmConvFactor[Δk_,len_]:=(((Sin[Δk len/2])^2)/(((Δk len/2))^2));

θRange=2 Degree;
lenTst=5 mm;
Plot[ΔkTypeIdoubling[λFund,θpmI+θTst Degree],{θTst,-θRange/Degree,θRange/Degree},PlotRange->All,Frame->True,FrameLabel->{"Δθ (deg)","Δk (/μm)"}]
Plot[pmConvFactor[ΔkTypeIdoubling[λFund,θpmI+θTst Degree],lenTst],{θTst,-θRange/Degree,θRange/Degree},PlotRange->All,Frame->True,FrameLabel->{"Δθ (deg)","relative yield"}]

tolerance to wavelength detuning (e.g. for short pulse)

Calculate the bandwidth such that (sinc^2)(Δk L/2)<50%.
For a pulse with 1/(e^2) half widths of τ, Δω, Δω = 2/τ
In terms of λ: Δλ=((Δω)/(Subscript[ω, 0])) (Subscript[λ, 0])=Δω(( ((Subscript[λ, 0])^2))/(2π c)), so Δλ=(( ((Subscript[λ, 0])^2))/(2π c))(2/τ)

ΔλBW[τ_,λ_]:=((λ^2)/(π c τ));
τPulse=20.fs;
Print["fundamental input pulse duration: ",τPulse/fs," fs. \r wavelength bandwidth: ",ΔλBW[τPulse,λFund]/nm," nm"];
Print["ideal harmonic pulse duration: ",τPulse/(Sqrt[2])/fs," fs. \r wavelength bandwidth: ",ΔλBW[τPulse/(Sqrt[2]),λFund/2]/nm," nm"];

Effect of limited phase matching bandwidth is on the converted spectrum

θTweak changes incident angle away from optimum

gauss[ω_,Δω_,ωCtr_]:= Exp[-2 ((((ωCtr-ω)/Δω))^2)];

ω0=((2 π c)/(λFund));
ΔωBW=((ΔλBW[τPulse,λFund])/(λFund))ω0;
θTweak=0.0 Degree;
λRange=150nm;

cryLen=5mm;
Print[" for a crystal length of: ",cryLen/mm," mm, the relative group delay is: ",ΔτgroupDelay[λFund,θpmI,cryLen]/fs," fs"];
Plot[{
gauss[2 π c/(λ um),ΔωBW,ω0],
gauss[2π c/(λ um),ΔωBW,ω0]pmConvFactor[ΔkTypeIdoubling[λ,θpmI+θTweak],cryLen]
},
{λ,(λFund-λRange)/um,(λFund+λRange)/um},Frame->True,AxesOrigin->{λFund/um,0},FrameLabel->{"λ (μm)","yield"}]

Plot[{Exp[-(((t fs))^2)/((τPulse)^2)],0.5Exp[-(((t fs-50))^8)/(((ΔτgroupDelay[λFund,θpmI,cryLen]/2))^8)]},{t,-150,200},PlotRange->All,PlotStyle->{Red,Blue},Frame->True,FrameLabel->{"time (fs)",None}]