ABCD analysis

by C. Durfee, CSM for PHGN480 Laser Physics

initializations

constants

$TextStyle={FontFamily->"Times",FontSize->14,FontWeight->"Bold"}

m=100;J=1;
cm=(10^(-2))m;mm=(10^(-3))m;um=(10^(-6))m;nm=(10^(-9))m;
uJ=(10^(-6))J;mJ=(10^(-3))J;
m=1;W=1;
cm=(10^(-2))m;mm=(10^(-3))m;um=(10^(-6))m;nm=(10^(-9))m;
J=1;uJ=(10^(-6))J;nJ=(10^(-9))J;mJ=(10^(-3))J;
s=1;ns=(10^(-9))s;fs=(10^(-15))s;
c=3 (10^8);
ϵ0=8.8 (10^(-12));

define abcd elements: 
for free space, using only the definition that includes the refractive index

mFreeSpace[z_,n_]:=(({{1, z/n}, {0, 1}})); (* this is the same result as going through a window *)
mSimpleLens[f_]:=(({{1, 0}, {((-1)/f), 1}}));
mSphericalMirror[R_]:=(({{1, 0}, {((-2)/R), 1}}));
mSimpleLensX[f_,θ_]:=(({{1, 0}, {((-1)/(f Cos[1.8 θ])), 1}}));(* tilt around y-axis, fx changes *)
mSimpleLensY[f_,θ_]:=(({{1, 0}, {((-1)/(f Cos[0.8 θ])), 1}}));(* tilt around x-axis, fy changes *)
mSphericalMirrorX[R_,θ_]:=(({{1, 0}, {((-2)/(R  Cos[θ])), 1}}));
mSphericalMirrorY[R_,θ_]:=(({{1, 0}, {((-2Cos[θ])/R), 1}}));
mGradientLens[γ_,dz_]:=(({{Cos[γ dz], Sin[γ dz]/γ}, {-γ * Sin[γ dz], Cos[γ dz]}}));
mGain[wpump_,gL_,λ0_]:=(({{1, 0}, {-I((λ0 gL)/(π (wpump^2))), 1}}));

additional matrices from Siegman
Tangential plane is where the beam width is measured along the plane of incidence (shortens focal length of tilted mirror).
Sagittal is the other one.

mDielFlatInterface[n1_,n2_]:=(({{1, 0}, {0, (n1/n2)}}));

mFlatInterfaceTang[thInc_,n1_,n2_]:=(({{((Cos[ArcSin[n1 Sin[thInc]/n2]])/(Cos[thInc])), 0}, {0, ((Cos[thInc])/(Cos[ArcSin[n1 Sin[thInc]/n2]]))}}));
mFlatInterfaceSagittal[thInc_,n1_,n2_]:=(({{1, 0}, {0, 1}}));

functions for ABCD raytrace

may have to fix to account for index inside medium

plotting function: plot a propagation segment
y0 = rayIn[[1]] initial height
y' = rayIn[[2]] slope

line to plot: 
y=y0+(z - z0) y' = rayIn[[1]] + (z - z0)  rayIn[[2]]

Don't have beamGradient working, but it should apply to rays

plotRay[rayIn_,zStart_,dz_,zUnit_,wUnit_,pltCol_]:=Plot[(rayIn[[1]]+(z zUnit-zStart)rayIn[[2]])/wUnit,{z,zStart/zUnit,(zStart+dz)/zUnit},PlotRange->All,PlotStyle->pltCol];

(* plotBeamWGradient[qinvIn_,zStart_,dz_,zUnit_,wUnit_,pltCol_]:=Plot[beamW[qinvOptic[qinvIn,GradientLens[alphaGain,(z zUnit-zStart)]]]/wUnit,{z,zStart/zUnit,(zStart+dz)/zUnit},PlotRange->All,PlotStyle->pltCol]; *)

y0=0.6 mm;yPrime=((0.1mm)/(10 mm));
col=Blue;
ray={y0,yPrime};
plotRay[ray,0,300 mm,mm,um,col]

define functions for Gaussian Beam ABCD

calculations are done with inverse of complex beam parameter q
R is radius of wavefront curvature
Z is a Rayleigh range calculated from the local beam size (even though beam may not be collimated there)

alternate definitions that include refractive index should replace ones without n throughout code
sign convention may be inconsistent with the -i ω t convention: 
may need to be: 
(1/q)=(1/R)+(i/z)

qinv[R_,Z_]:=(1/R)-I (1/Z);
z0[w_,n_]:=((n π (w^2))/lam0);

functions to extract beam info from qinv

beamWQ[n_,qinvIn_]:=(Sqrt[-(lam0/(n π Im[qinvIn]))]);
beamRQ[qinvIn_]:=If[Re[qinvIn]==0,∞,1/Re[qinvIn]];

beam prop equations based on distance from waist (z) and waist size (w0).

beamWz[z_,w0_,n_]:=w0(Sqrt[1+(z^2)/((z0[w0,n])^2)]);
beamRz[z_,w0_,n_]:=z((((((z0[w0,n])/z)))^2)+1);
beamηz[z_,w0_,n_]:=ArcTan[z/z0[w0,n]];

these expressions give the size and location of the waist, given the size and radius of curvature at a given position

Clear[lam0]
Solve[{w1 == beamWz[z, w0, n], r0 == beamRz[z, w0, n]}, {z, w0}]

zWaist[w1_,r1_,n_]:=(((π^2) r1 (w1^4))/(((((lam0/n)))^2) (r1^2)+(π^2) (w1^4)));
w0Waist[w1_,r1_,n_]:=(((lam0/n) Abs[r1] w1)/(Sqrt[((((lam0/n)))^2) (r1^2)+(π^2) (w1^4)]));

propagation

for an optic: qinv = ((C + D qinv)/(A + B qinv))

qinvZ[z_,n_,qinvIn_]:=(((Subscript[mFreeSpace[z,n], [[2,1]]])+(Subscript[mFreeSpace[z,n], [[2,2]]]) qinvIn)/((Subscript[mFreeSpace[z,n], [[1,1]]])+(Subscript[mFreeSpace[z,n], [[1,2]]]) qinvIn)) ;
qinvOptic[qinvIn_,MatOptic_]:=(((Subscript[MatOptic, [[2,1]]])+(Subscript[MatOptic, [[2,2]]]) qinvIn)/((Subscript[MatOptic, [[1,1]]])+(Subscript[MatOptic, [[1,2]]]) qinvIn)) ;

qinvZ[z,1,q1]

zWaistQ[qinvIn], w0WaistQ[qinvIn] give the size and location of the waist, given the value of q at a given position
(q can be calculated from the known size and radius of curvature)

qinvZ[z, n, qReIn + I qImIn]

waist is at z where real part of this expression = 0

Expand[(i qImIn+qReIn)(1+(((-i qImIn+qReIn) z)/n))]

Solve[qReIn+(((qImIn^2) z)/n)+(((qReIn^2) z)/n)==0,z]

zWaistQ[n_,qinvIn_]:=-((n Re[qinvIn])/((Abs[qinvIn])^2));
w0WaistQ[n_,qinvIn_]:=beamWQ[n,qinvZ[zWaistQ[n,qinvIn],n,qinvIn]];

stability of an oscillator cavity

In Svelto: -1≤((A + D)/2)≤ 1, with edges of this zone being marginally stable. 
Using the parameter s, defined through
s=1-(((((A + D)/2)))^2), 
the cavity is stable if 0 ≤ s ≤ 1, and maximum stability is at s = 1, marginally stable at s = 0.

stability[MatCav_]:=1-((((((Subscript[MatCav, [[1,1]]])+(Subscript[MatCav, [[2,2]]]))/2)))^2);
StableModeSize[MatCavity_,n0_]:=If[Chop[stability[MatCavity]]>0,(Sqrt[(lam0/(π n0))])(((Abs[(Subscript[MatCavity, [[1,2]]])])^(1/2))/((stability[MatCavity])^(1/4))),0];

plotting function: plot a propagation segment

plotBeamW[n_,qinvIn_,zStart_,dz_,zUnit_,wUnit_,pltCol_]:=Plot[beamWQ[n,qinvZ[(z zUnit-zStart),n,qinvIn]]/wUnit,{z,zStart/zUnit,(zStart+dz)/zUnit},PlotRange->All,PlotStyle->pltCol];
plotBeamR[n_,qinvIn_,zStart_,dz_,zUnit_,wUnit_,pltCol_]:=Plot[1/(beamRQ[qinvZ[(z zUnit-zStart),n,qinvIn]]/wUnit),{z,zStart/zUnit,(zStart+dz)/zUnit},PlotRange->All,PlotStyle->pltCol];
plotBeamWGradient[γ_,qinvIn_,zStart_,dz_,zUnit_,wUnit_,pltCol_]:=Plot[beamWQ[1,qinvOptic[qinvIn,mGradientLens[γ,(z zUnit-zStart)]]]/wUnit,{z,zStart/zUnit,(zStart+dz)/zUnit},PlotRange->All,PlotStyle->pltCol];

Summary of functions:
qinv[R, Z]	 	
  calculate inverse q from radius R and Z=n π (w^2)/λ, where w is local beam size. beam is not collimated unless R=∞
z0[w, n] 		
   z0=n π (w^2)/λ this is Rayleigh range when the beam is at a waist
beamWQ[n, qinv] 	
  calculate the local beam w from given inverse q. Have to specify local refractive index. 
beamRQ[qinv]		
  calculate the local beam R from given inverse q. Index not needed.

beamWz[z, w0, n]	
  beam radial width as function of distance z from waist of radius w0. This is standard definition for Gaussian beam. 
beamRz[z, w0, n]	
  beam radius of curvature as function of distance z from waist of radius w0, using standard definition for Gaussian

zWaist[w, r, n]		
  calculate position of waist given local w, R, and n
w0Waist[w, r, n]	
  calculate size of waist given local w, R, and n

qinvZ[z, n, qinvIn]	
  do ABCD propagation as function of z. gives new qinv
qinvOptic[z, qinvIn]	
  gives new qinv after going through an optical system specified by input ABCD 

zWaistQ[n, qinv]	
  calculate position of waist for given qinv
w0WaistQ[n, qinv]	
  calculate min waist radius for given qinv

stability[cavityABCD]	
  stability parameter for a round-trip cavity. The function is defined so that a stable cavity gives a value between 0 and 1, with 0 being at the boundary of instability
StableModeSize[cavityABCD, n]
  stable mode radius at starting position in cavity. Permute matrices to get mode size at other positions.

plotBeamW[n, qinvIn, zStart, dz, zUnit, wUnit, plotColor]
plotBeamR[n, qinvIn, zStart, dz, zUnit, wUnit, plotColor]
  returns plot of the beam size (or beam radius of curvature) vs z, given initial qinv and index
  To use: calculate qInvIn first. This is the q-parameter for the starting segment at zStart. The function will plot from zStart to zStart + dz. The axes will be scaled with the units chosen in zUnit and wUnit. Pick plotColor as one of the standard built-in colors (Red, Blue, etc.) To make a plot through a system, calculate the q’s at the beginning of each free-space section, put the return quantity of the plotting function into a variable, then use Show[ ] to plot all segments at once. The option GridLines for Show can be useful to mark positions of lenses etc.

examples

test raytrace

single lens - collimated input

nSegments=2;
rayStart={0.5,0};
z1=10mm;z2=30mm;f1=20mm;
ray1=mSimpleLens[f1].mFreeSpace[z1].rayStart;
rayEnd=mFreeSpace[z2].mSimpleLens[f1].mFreeSpace[z1].rayStart;
Print[rayEnd]

pltLst=Range[nSegments];
pltLst[[1]]=plotRay[rayStart,0,z1,mm,mm,col];
pltLst[[2]]=plotRay[ray1,z1,z2,mm,mm,col];
Show[pltLst]

The following defines a collimated set of rays for the input to the system.

z1=10mm;z2=100mm;f1=20mm;

beamWid=2.;
nRays=10;
dy=beamWid/(nRays-1);
rayHts=Table[(i-1) dy-beamWid/2,{i,1,nRays}];
rayAngles=0 Range[nRays];
raySet0=Transpose[{rayHts,rayAngles}];

This defines a point source in a similar way - starting position at a certain height, then define a set of rays.

beamHt=0;
nRays=10;
angleRange=0.1;
dA=angleRange/(nRays-1);
rayHts=ConstantArray[beamHt,nRays];
rayAngles=Table[(i-1) dA-angleRange/2,{i,1,nRays}];
raySet0=Transpose[{rayHts,rayAngles}];

The following lines apply the matrix dot multiplication to each ray in the list in the raySets
the # is a placeholder for the function argument, the & identifies what precedes it as a "pure function", and the /@ performs the Map[ ] function that does the . product on each of the vectors in the list. 

To do a trace, you make a ray set for the beginning of each propagation segment. 
Here raySet0 is the beginning, propagating forward from z = 0 to z = z1
raySet1 results from forward propagation to the lens, then propagation through the lens. 
raySet2 is not needed for the display of the rays, but could be used to numerically inspect what the rays are like at a certain plane after the lens.

z1=10mm;

raySet1=(mSimpleLens[f1].mFreeSpace[z1].#)&/@raySet0;
raySet2=(mFreeSpace[z2].#)&/@raySet1;

pltLst=ConstantArray[0,{nSegments,nRays}];
pltLst[[1]]=plotRay[#,0,z1,mm,mm,col]&/@raySet0;
pltLst[[2]]=plotRay[#,z1,z2,mm,mm,col]&/@raySet1;
Show[pltLst,GridLines->{{z1/mm},None},GridLinesStyle->Red]

ABCD propagation examples

lam0=0.589 um;w0=0.16 mm;R0=∞;col=Blue;
nTst=1;
p1=plotBeamW[nTst,qinv[R0,z0[w0,nTst]],0,500 mm,mm,um,col];
Show[p1,AxesOrigin->{0,0}]

focusing beam

lam0=0.8 um;w0=0.5 mm;R0=∞;col=Blue;
n1=1;
Print["input beam Rayleigh range: ", z0[w0,n1]/m," m"];

fTst=50mm;
z2=150mm;
qInv0=qinv[R0,z0[w0,n1]];
qInvLens=qinvOptic[qInv0,mSimpleLens[fTst]];

pNoDiel=plotBeamW[n1,qInvLens,0,z2,mm,um,col];
Show[pNoDiel,AxesOrigin->{0,0}]
Print["z position of waist (no glass): ",zWaistQ[n1,qInvLens]/mm," mm"];Print["min spot radius (no glass): ",w0WaistQ[n1,qInvLens]/um," microns"];

Class exercises

focusing into a dielectric medium

Suppose a beam is focused into a higher-index medium. How does the focal spot size and the location of the focus change? 

Construct a matrix that represents (in order that the beam sees) a lens, translation by z1 through vacuum, dielectric interface, then translation by z2 through glass.  
- Make a plot of the propagation of the beam waist through the system by making plots for region 1: z = 0 to z1, z = z1 to z2. Then use Show[ ] to display all together. 
- use pre-defined functions to compute the z position of the waist and the spot size inside the medium.

relay imaging: single lens

A single lens can be used to make an image of an object. When it is a Gaussian beam passes through the lens, how are the “object” and “image” planes related?  
The example below shows a simple 1:1 imaging sytem with object and image distances the same. 
- try different starting spot sizes and object distances (z1). How are the spot sizes related? Is the complex field imaged directly?

lam0=0.8 um;
w0=100. um;
R0=∞;
col=Blue;
n1=1;
Print["input beam Rayleigh range: ", z0[w0,n1]/cm," cm"];

fTst=500mm;
z1=1000mm;
z2=((((1/fTst)-(1/z1)))^(-1));
Print["image plane: ",z2/mm," mm from lens"];
Print["magnification: ",mag=z2/z1];
qInv0=qinv[R0,z0[w0,n1]];
qInv1=qinvZ[z1,n1,qInv0];
qInvLens=qinvOptic[qInv1,mSimpleLens[fTst]];

p1=plotBeamW[n1,qInv0,0,z1,mm,um,col];
p2=plotBeamW[n1,qInvLens,z1,z2,mm,um,col];
Show[p1,p2,AxesOrigin->{0,0},GridLines->{{z1/mm},{w0/um,mag w0/um}}]

ABCD analysis of cavities

2 mirror cavity

stability calculation for a 2-mirror cavity:

Print[MatrixForm[mSphericalMirror[R1]],MatrixForm[mFreeSpace[L,1]],MatrixForm[mSphericalMirror[R2]],MatrixForm[mFreeSpace[L,1]]]

mSys=mSphericalMirror[R1].mFreeSpace[L,1].mSphericalMirror[R2].mFreeSpace[L,1];
MatrixForm[mSys]

calculate ((A + D)/2):

eq1=FullSimplify[(mSys[[2,2]]+mSys[[1,1]])/2]

define (Subscript[g, 1])=1-L/(Subscript[R, 1]) and (Subscript[g, 2])=1-L/(Subscript[R, 2])

eq2=eq1/.{R1->L/(1-g1),R2->L/(1-g2)};
FullSimplify[eq2]

since the stability condition is -1≤((A + D)/2)≤1, -1≤2 (Subscript[g, 1])(Subscript[g, 2])-1≤1, or  0≤(Subscript[g, 1])(Subscript[g, 2])≤1
In terms of the shifted stability (s) parameter defined above, s=1-(((((A + D)/2)))^2)=1-(((2(Subscript[g, 1])(Subscript[g, 2]) -1))^2)=1-(((2(Subscript[g, 1])(Subscript[g, 2]) -1))^2)
This shifted parameter is better for plotting since the stability boundaries are always at 0.

FullSimplify[1-(((2 g1 g2-1))^2)]

stabilityg1g2[g1_,g2_]:=-4 g1 g2 (g1 g2-1);

Plot[-4 g1g2 (-1+g1g2),{g1g2,-.1,1.1}]

The following plot is a variation on the standard stability plot. The dark edges are at marginal stability, and the most stable regions are the light bands in the middle.

gMax=3;
pStability=ContourPlot[stabilityg1g2[g1,g2],{g1,-gMax,gMax},{g2,-gMax,gMax},PlotRange->{0,1}]

unequal curved mirrors

Here we define a system matrix that is a function of the curvatures and cavity length

Clear[mSys];
mSys[R1_,R2_,L_]=mSphericalMirror[R1].mFreeSpace[L,1].mSphericalMirror[R2].mFreeSpace[L,1];

lam0=632.8nm;
R1Tst=100 cm;R2Tst=45 cm;
Lmax=Min[R1Tst+R2Tst,3];

LTst=R2Tst/2;
wM1=StableModeSize[mSys[R1Tst,R2Tst,LTst],1];
wM2=StableModeSize[mSys[R2Tst,R1Tst,LTst],1];

Print["red: mode size at M1"];
Print["blue: mode size at M2"];
Plot[{
StableModeSize[mSys[R1Tst,R2Tst,L],1]/mm,StableModeSize[mSys[R2Tst,R1Tst,L],1]/mm,
stability[mSys[R1Tst,R2Tst,L]]},{L,0,Lmax},
PlotRange->{0,1},PlotStyle->{Red,Blue,GrayLevel[0]},GridLines->{{LTst},None}]

Print["for cavity length: ",LTst/mm," mm"];
Print["mode dia at M1: ",wM1/um," um"];
Print["mode dia at M2: ",wM2/um," um"];

The plot above is a slice through the stability map:

pLslice=ParametricPlot[{1-L/R1Tst,1-L/R2Tst},{L,0,Lmax},PlotStyle->Red];
Show[pStability,pLslice]

questions to work through:

a two-mirror resonator is formed by a convex mirror of radius R1= -1m and R2 = 1.5m. What is the maimum possible mirror separation if this is to mremain a stable resonator?

A resonator consists of two plane mirrors with a positive lens instered between the two mirrors. If the focal length of the lens is f, and L1 and L2 are the distances of the lens from the two mirrors, calculate: 
1. the spot size at the lens position, and the spot sizes at each mirror
2. the conditions under which the cavity is stable.
3. To make our Nd:YAG laser operate, we will need to construct a resonator of this type. One end mirror is where the gain medium is, the other is the output coupler. Under the constraints that the cavity length can be no more than 10cm to fit on the cage rails, and that the mode radius should be approximately 75μm 1/e2 radius, find a focal length and cavity length for the cavity.  
4. Make a plot of the mode size through the optical cavity.

sequence of lenses system

sequence of lenses as a model of a resonator

unit cell is lens (f), then propagate (L)
confocal is L = 2 f
concentric is L = 4f

Clear[mPropLens]

If the sequence is stable, the beam waist should be halfway between the lenses. 
  system:  distance = L/2 to lens, then apply lens of focal length f, distance = L/2
 Note that the matrix below is only 1/2 of a round trip.

mPropLens[f_,len_]=mFreeSpace[len/2,1].mSimpleLens[f].mFreeSpace[len/2,1];
MatrixForm[mPropLens[f,len]]

test whether one round trip would be stable

lam0=0.8 um;
w0=0.1 mm;
R0=∞;col=Blue;
nTst=1;
fTst=100.mm;
lTst=50mm;

mRT=mPropLens[fTst,lTst].mPropLens[fTst,lTst];
Print["stability parameter: ",stability[mRT]];
wStable=StableModeSize[mRT,1];
Print["stable mode size: ",wStable/mm," mm"];
Print["det: ",Det[mRT]];
Print["eigenvals: ", Eigenvalues[mRT] ];
Print["trace: ",Tr[mRT] ];
qInv0=qinv[R0,z0[w0,nTst]];

plot stability vs lens spacing

lenMin=fTst/10;
lenMax=4.5fTst;
Plot[stability[mPropLens[fTst,l mm].mPropLens[fTst,l mm]],{l,lenMin/mm,lenMax/mm}]

Plot evolution of the beam through this system
for plotting purposes, it is easiest to have the boundaries of the propagation sequence end just after the lens.

mPropLens[f_,len_]=mSimpleLens[f].mFreeSpace[len,1];

define input beam parameters

w0=0.13 mm;
R0=∞;col=Blue;

qInv0=qinv[R0,z0[w0,nTst]];

set up plotting segments

nRT=3;
plots=qInvLst=0 Range[2 nRT+1];
qInvLst[[1]]=qInv0;
plots[[1]]=plotBeamW[nTst,qInvLst[[1]],0,lTst/2,mm,um,col];

Show[plots[[1]]]

For[i=2,i≤2nRT+1,i++,
If[i==2,
qInvLst[[i]]=qinvOptic[qInvLst[[1]],mSimpleLens[fTst].mFreeSpace[lTst/2,1] ],
qInvLst[[i]]=qinvOptic[qInvLst[[i-1]],mPropLens[fTst,lTst]]
];
plots[[i]]=plotBeamW[nTst,qInvLst[[i]],(i-1)lTst-lTst/2,lTst,mm,um,col];
];

gridLst=(Range[2nRT]-1/2)lTst/mm;
Show[plots,GridLines->{gridLst,None}]

exercise: 
vary the input spot size and make some plots to see what happens to the mode size
is the mode confined?

sequence of lenses as a waveguide

For a weak lens, a periodic sequence will change the divergence slowly over a unit cell. A gradient index waveguide can be modeled this way.

Clear[mPropLens]

mPropLens[f_,len_]=mSimpleLens[f].mFreeSpace[len,1];
MatrixForm[mPropLens[f,len]]

w0=0.025 mm;

lam0=0.8 um;
R0=∞;col=Blue;
nTst=1;
fTst=5.mm;
lTst=1mm;

mRT=mPropLens[fTst,lTst].mPropLens[fTst,lTst];
Print["stability parameter: ",stability[mRT]];
qInv0=qinv[R0,z0[w0,nTst]];

nRT=20;
plots=qInvLst=0 Range[2 nRT];
qInvLst[[1]]=qInv0;

For[i=1,i≤2nRT,i++,
If[i>1,
qInvLst[[i]]=qinvOptic[qInvLst[[i-1]],mPropLens[fTst,lTst]]
];
plots[[i]]=plotBeamW[nTst,qInvLst[[i]],(i-1)lTst,lTst,mm,um,col];
];
gridLst=Range[2nRT]lTst/mm;
Show[plots,GridLines->{None,None},AxesOrigin->{0,0}]

in the example above, there is substantial change in the mode size in between the lenses. Change to longer focal lengths to observe the mode size evolution.

quadratic index profile wave guide

γTst=200;
LTst=10cm;
w0=50 um;
Print["guide period: ",π/γTst/mm//N," mm"];
Print["stable mode size: ", (Sqrt[lam0/(π γTst)])/um," μm"];
lam0=0.8 um;
R0=∞;col=Blue;
qInv0=qinv[R0,z0[w0,nTst]];

pGrad=plotBeamWGradient[γTst,qInv0,0,.10,mm,um,col];
Show[pGrad,AxesOrigin->{0,0}]