Paraxial Beam, Gaussian Basics

ECE 5368/6358                  han q le - copyrighted
Use solely for students registered for UH ECE 6358/5368 during courses - DO NOT distributed (copyrighted materials)

Utility

EGauss[{z_,ρ_,t_},{λ_,w_,ω_}]:=Module[{k,Q}
,k=2 π/λ ; Q=z-I  ((π (w^(2 )))/λ) ;
(( w k)/(I Q))(E^(I  (k ( z+  ((  (ρ^2))/(2 Q)))- ω t) ))
] ;

0. Introduction

Why study Gaussian beam?
- Historically, long before the ubiquitous diode lasers, early lasers were often gas lasers
- Gas lasers are built with macroscopic external cavity with concave mirrors 
- Gaussian beam is the output profile these gas lasers

(http://en.wikipedia.org/wiki/Gaussian_beam)

- Understanding Gaussian beam is essential to design these laser cavities (Hermite-Gaussian profiles are the transverse mode of most of these lasers), as well to use of these laser beams outside the cavity.
- More than just for the practical reason dealing with lasers, the mathematics of Gaussian beam is manageable and instructive to gain insight and understanding of the basics of coherent light propagation, hence it is a worthwhile topic in itself. 
- Concepts such as beam waist and beam divergence are basic to any beam and light propagation in general. As a light beam undergoes transformation through optics, especially Fourier optics (since the Fourier transform of a Gaussian is also a Gaussian), it is important to determine its properties, and Gaussian beam provides a very useful approximation even for non-Gaussian beams.
- Nowadays, most common coherent light beams are the outputs of diode or semiconductor-waveguide laser or optical fiber. They have transverse mode profiles that are not necessarily Gaussian, but simple Gaussian-based formula can still be used as a rough approximation or quick rule-of-thumb estimation. The calculation with Gaussian beam is simple and easy to do.

1. Paraxial approximation for Gaussian beam (can skip)

Earlier, we review 2D Gaussian beam. Here, we will study 3D (most common) Gaussian beam.

A beam coming out from a small spot acts like a spherical wave. A beam coming out from a large spot acts like a plane wave. How do we account for these extreme behaviors with a single, simple formulation? We don't want to have to use a different "approximate" for each situation.
Beam propagating along an axis: paraxial beam. The solution can be separated in two parts: traveling logngitudinal part (along the axis) and lateral part. In free space, an important type of paraxial beam is Gaussian beam.

Hemholtz's equation:
(∇^2)U + ((Subscript[k, 0])^2)ϵ(r) U = 0  ==>  ((Subscript[∇, T])^2)U +  (∂_(z,z))U  +  ((Subscript[k, 0])^2)ϵ(r) U = 0 
Separate by part would suggest U -> (e^(i k z )) F . But we have to remember that at infinity, the plane wave approx. is good over a small region (compared with the beam spread), but a spherical wavefront is the correct approximate over the entire beam.

1.1 Spherical wave review

<<"VectorAnalysis`"

Expand[Laplacian[U[r,θ,ϕ],{r,θ,ϕ},"Spherical"]]

For spherical wave without angular dependence:

Expand[Laplacian[U[r],{r,θ,ϕ},"Spherical"]]

The Helmholtz Eq.:

DSolve[((2 U'[r])/r)+U''[r]+ ((Subscript[k, 0])^2)ϵ U[r]==0,U[r],r]

There are two solutions:
(1/r)(e^(i k r )) or (1/r)(e^(-i k r))
for radiating wave or collapsing wave: these wave have singularity at origin.
((2 U'[r])/r)+U''[r]+ ((Subscript[k, 0])^2)ϵ U[r]==0

1.2 In cylindrical coordinate

We can solve the Helmholtz Eq. by part for the cylindrical geometry		
Let             U[ρ,z ]=(e^(i k z))A[ρ ]

Expand[Laplacian[(E^(I k z))A[ρ ] ,{ρ,θ,z},"Cylindrical"]]

1.2.1 Exact solution: Bessel J0

The equation becomes:
- A(ρ) (k^2)+(( A'(ρ))/ρ)+ A''(ρ)+ ϵ (Subscript[k, 0]^2) A[ρ]==0  
(ρ^2) A''[ρ] +ρ A'[ρ]  +  (ϵ (Subscript[k, 0]^2)-(k^2)) (ρ^2) A[ρ]==0 
Change variable: for  (ϵ (Subscript[k, 0]^2)-(k^2)) > 0:
(Sqrt[ϵ (Subscript[k, 0]^2)-(k^2)]) ρ = q ρ = ξ  
(ξ^2) A''[ξ] +ξ A'[ξ]  +  (ξ^2) A[ξ]==0

DSolve[(ξ^2) A''[ξ] +ξ A'[ξ]  +  (ξ^2) A[ξ]==0,A[ξ],ξ]

1.2.2 Review Bessel function

Plot[BesselJ[0,x],{x,0,100}, PlotRange->All]

The solution is Bessel (Subscript[J, 0])[q ρ] that is not singular at 0, (Subscript[Y, 0])[q ρ] diverges at zero.
Bessel beams are the only type of beam that does not change shape as it travels!  But a pure Bessel beam extends to infinity and goes down very slowly (1/(Sqrt[ρ])). These are not common beams.

BesselY[0,0]

1.2.3 Linear superposition solution

General solution is then:
∫c[q] (e^(i k z))(Subscript[J, 0])[ q ρ]dq where k=(Sqrt[ϵ (Subscript[k, 0]^2)-(q^2)])
What if  ϵ (Subscript[k, 0]^2)-(k^2) = 0?

DSolve[Derivative[1][A](ρ)+ρ Derivative[2][A](ρ)==0, A[ρ],ρ]

The non-plane wave solution (Subscript[c, 1]) log(ρ) is neither regular at zero or infinity: this is not a physical meaningful solution for a finite beam at ρ=0.

What if  (ϵ (Subscript[k, 0]^2)-(k^2))<0? We define (q^2) =(k^2) -ϵ (Subscript[k, 0]^2)
The equation:
(ξ^2) A''[ξ] +ξ A'[ξ]  -  (ξ^2) A[ξ]==0 
has solution Bessel I or Bessel K: neither one is regular: I is infinite at infinity, K is infinite at zero.

DSolve[(ξ^2) A''[ξ] +ξ A'[ξ]  -  (ξ^2) A[ξ]==0,A[ξ],ξ]

BesselI[0,100.]
BesselK[0,0]

This is a problem: what if the wavefront at the plane z=0 has q component up to ∞? We notice that if q > (Sqrt[ϵ ((Subscript[k, 0])^2)]), the only way the Helmholtz equation is satisfied is that the z-dependence plane wave is: (e^(±κ z))  where κ=(Sqrt[(q^2) -ϵ ((Subscript[k, 0])^2)]). This is an acceptable solution in either + z or - z direction as long as it goes to zero at infinity: We call this evanescence wave. It is meaningful only for a short distance of z. Thus, for a propagating wave solution, this component vanishes eventually and can be ignored. But it can be important at near field solution, which is the basis for near field spectroscopy.

1.3 Gaussian beam paraxial propagating solution

1.3.1 Gaussian

The general propagating wave solution is a linear superposition of Bessel beams.
How do we form a linear superposition of Bessel beams?  
                      (∫_0^((Subscript[k, 0])(Sqrt[ϵ ])))c[q] (e^(i k z))(Subscript[J, 0])[ q ρ]dq  where k=(Sqrt[ϵ (Subscript[k, 0]^2)-(q^2)])            (1.3.1)
Notice that we have a cut-off at (Subscript[k, 0])(Sqrt[ϵ ]). Near this value, (e^(i k z)) is nearly constant: the wave is not propagating in the z direction. The paraxial solution for typical laser beam should have k~(Subscript[k, 0])(Sqrt[ϵ  ]) and q very small.

We have formed a linear superposition of plane wave. We'll do the same for Bessel waves; but we need to assume certain function c[q] . Obviously, c[q] can not be finite at q=0: we have singular solution! so we'll choose c[q]= q (e^(-(q^2))).

Fs[ρ_,z_, k1_,Cfunc_]:=NIntegrate[Cfunc[q]*Exp[ I z (Sqrt[k1^2-q^2])]*BesselJ[0,q ρ],{q,0,k1}]

Cfunc[q_]:=2*q*Exp[-q^2];
Plot[ Cfunc[q],{q,0,5}]

What we find, interestingly is that by choosing c[q]= q (e^(-(q^2))), along with paraxial approximation, we get a Gaussian profile.

q= .; w= .;

(1/(  (Sqrt[2π])))(∫_0^∞)  w q (E^(- (q^2) (w^2)/4 ))BesselJ[0, ρ q]\[DifferentialD]q

(2/(π (w^2)))(∫_0^∞)2 π  (E^(-((2(ρ^2))/(w^2))))ρ \[DifferentialD]ρ

(∫_0^∞)(E^(-(1/4) (q^2) (w^2))) (E^(I (Sqrt[((Subscript[k, 1])^2) - (q^2)]) z))  q  (Subscript[J, 0])(q ρ)\[DifferentialD]q

Thus:              (w/(  (Sqrt[2π])))(∫_0^∞)(e^(-(1/4) (q^2) (w^2))) q  (Subscript[J, 0])(q ρ)dq =(Sqrt[(2/(π ))]) (1/w)(e^(-((ρ^2)/(w^2))))
We see that with the approximation (∫_0^∞)(e^(-(1/4) (q^2) (w^2))) q (w^2) (Subscript[J, 0])(q ρ)dq with the integration limit to infinity instead of (Subscript[k, 0])(Sqrt[ϵ ]), we have a closed-form Gaussian solution. This is the Gaussian profile at z=0 that we want! To extend this approximation to any z, we write:  
(Sqrt[((Subscript[k, 1])^2) - (q^2)]) = (Subscript[k, 1]) (Sqrt[1 - ((q^2)/((Subscript[k, 1])^2))])≈  (Subscript[k, 1])-((q^2)/(2 (Subscript[k, 1])  )) : this is the paraxial approximation.
(e^(i  (Subscript[k, 1]) z))  (w/(  (Sqrt[2π])))∫(e^(-(1/4) (q^2) (w^2)-i ((q^2)/(2 (Subscript[k, 1]))) z)) q  (Subscript[J, 0])[ q ρ]dq
= (e^(i  (Subscript[k, 1]) z))  (w/(  (Sqrt[2π])))∫(e^(-(1/4) (q^2) ((w^2)+i ((2 z)/( (Subscript[k, 1]))) ))) q  (Subscript[J, 0])[ q ρ]dq

(w/(  (Sqrt[2π])))(∫_0^∞)  q (E^(-(1/4) (q^2) ((w^2)+I ((2 z)/( (Subscript[k, 1]))) )))BesselJ[0, ρ q]\[DifferentialD]q

The approximated solution is:
             E[ρ,θ,z]=(1/(Sqrt[2 π]))((w  (Subscript[k, 1]))/(i ( z-i  (Subscript[k, 1]) (w^2)/2)))(e^(i  (Subscript[k, 1]) z)) (e^(i  (( (Subscript[k, 1]) (ρ^2))/(2(z-i  (Subscript[k, 1]) (w^2)/2))))) 
                              = (1/(Sqrt[2 π]))(( w (Subscript[k, 1]))/(i Q[z]))(e^(i (Subscript[k, 1]) (z+((  (ρ^2))/(2 Q[z]))) ))                                    (1.3.2a)
where                                                     Q[z]=z-i  (((Subscript[k, 1]) (w^2))/2)               (1.3.2b)
What is the meaning of this Q[z] term?  (((Subscript[k, 1])(w^(2 )))/2)=((π (w^(2 )))/λ)= (Subscript[ζ, 0]) : defined as Rayleigh range.  So Q[z] is a quantity of unit length, but has an imaginary component: Q[z]=z-i (Subscript[ζ, 0]). 
To include the time:
                                         E[ρ,θ,z,t]=(1/(Sqrt[2 π]))(( w (Subscript[k, 1]))/(i Q[z]))(e^(i ((Subscript[k, 1])z-ω t) ))  (e^(i (Subscript[k, 1])((  (ρ^2))/(2 Q[z])) ))    (1.3.3)

Let's plot electric field (or H field):

λ=0.6328 ;ω = 2 π; n=1;w= 0.5 ;(Subscript[k, 1])=((2π n)/λ);y=0; 
t=0.; 
beamplot=Plot3D[ Re[((w  (Subscript[k, 1]))/(I ( z-I  (Subscript[k, 1]) (w^2)/2))) (E^(I  ((Subscript[k, 1]) z - ω t)))(E^(I (Subscript[k, 1])  (((x^2)+(y^2))/(2(z-I  (Subscript[k, 1]) (w^2)/2)))))],{x,-5,5},{z,0,5},PlotPoints->{50,100},
		PlotRange->All, Mesh->False (*, FaceGrids->All  *), BoxRatios->{2,1,2}, ViewPoint->{3,-1,2}]

λ=0.6328 ;ω = 2 π; n=1;w= 0.7 ;(Subscript[k, 1])=((2π n)/λ);y=0; 
t=0.; 
Animate[Plot3D[ Re[((w  (Subscript[k, 1]))/(I ( z-I  (Subscript[k, 1]) (w^2)/2))) (E^(I  ((Subscript[k, 1]) z - ω t)))(E^(I (Subscript[k, 1])  (((x^2)+(y^2))/(2(z-I  (Subscript[k, 1]) (w^2)/2)))))],{z,0,5},{x,-5,5},PlotPoints->{80,120},
		PlotRange->{{0,5},{-3.5,3.5},{-4,4}}, Mesh->False (*, FaceGrids->All  *), BoxRatios->{5,7,4}, ViewPoint->{3,-1,2}],{t,0,1,0.025}
,DisplayAllSteps->True,AnimationRunning->False]

Below is a code to use if the computation is slow:

λ=0.6328 ;ω = 2 π; n=1;w= 0.7 ;(Subscript[k, 1])=((2π n)/λ);y=0; 
gaussplot={} ;
For[i=0, i<40, i++,t=i/40. ;
AppendTo[gaussplot,
Plot3D[ Re[((w  (Subscript[k, 1]))/(I ( z-I  (Subscript[k, 1]) (w^2)/2))) (E^(I  ((Subscript[k, 1]) z - ω t)))(E^(I (Subscript[k, 1])  (((x^2)+(y^2))/(2(z-I  (Subscript[k, 1]) (w^2)/2)))))],{z,0,5},{x,-5,5},PlotPoints->{50,80},
		PlotRange->{{0,5},{-3.5,3.5},{-4,4}}, Mesh->False (*, FaceGrids->All  *), BoxRatios->{5,7,4}, ViewPoint->{3,-1,2}] ]
] ;

Manipulate[gaussplot[[i]],{i,1,Length[gaussplot],1} ]

Gaussian beam solution (in spite of being a paraxial approximation) is very good in describing circularly symmetric laser beam behavior: from tightly focussed, highly divergent beam to plane-wave-like beam.

Even if the beam is not circularly symmetric, but elliptic with x, y axis, it can still be used with separate x and y profile and all the key properties apply.

1.3.2 Phasefront

From the Gaussian beam description: 
                                         E[ρ,θ,z]=(1/(Sqrt[2 π]))(( w (Subscript[k, 1]))/(i Q[z]))(e^(i (Subscript[k, 1])(z+((  (ρ^2))/(2 Q[z]))) ))              (1.3.2a)
where                                         Q[z]=z-i  (((Subscript[k, 1]) (w^2))/2)                                  (1.3.2b)

We write separate the phase and the amplitude:
                        ((  w (Subscript[k, 1]))/(i Q[z]))(e^(i  (Subscript[k, 1]) (z+ (( (ρ^2))/(2Q[z]))) ))= (( w (Subscript[k, 1]))/(Abs[Q[z]] i))(e^(-(Subscript[k, 1]) (ρ^2) Im[(1/(2Q[z]))]  )) (e^(i  ( (Subscript[k, 1])(z+ (ρ^2) Re[(1/(2 Q[z]))])- Arg[Q[z]])  ))
Amplitude:                                                  (e^(-(Subscript[k, 1]) (ρ^2) Im[(1/(2Q[z]))]  ))            (1.3.3)
 Phase:                                             (e^(i  ( (Subscript[k, 1])(z+ (ρ^2) Re[(1/(2 Q[z]))])- Arg[Q[z]])  ))  (1.3.4)

The constant phase front (cophasar surface) is:
                                      (Subscript[k, 1])(z+ (ρ^2) Re[(1/(2 Q[z]))])- Arg[Q[z]]== constant        (1.3.5)
At z = 0, Re[(1/(Q[z]))]=0,  Arg[Q[z]]= constant, so the phase front is a plane: it is independent of ρ.

What about at farfield, where z-> ∞?

Using:                             Re[(1/(z-i  (Subscript[ζ, 0])))] = (z/((z^2)+  ((Subscript[ζ, 0])^2) )) ; Im[(1/(z-i  (Subscript[ζ, 0])))] = ((Subscript[ζ, 0])/((z^2)+  ((Subscript[ζ, 0])^2)))
                                         (Subscript[k, 1]) z +((  (Subscript[k, 1]) (ρ^2) z)/(2((z^2)+  ((Subscript[ζ, 0])^2)))) - ArcTan[((-(Subscript[ζ, 0]))/z)]== constant

As z-> ∞, ArcTan[((-(Subscript[ζ, 0]))/z)]-> -0, the phase front is given by 
                      ~ z + (( (ρ^2))/(2 z))= (Sqrt[(z^2)+(ρ^(2 ))])=r = constant: this is nearly a spherical wavefront.

Gaussian beam can describe both planewave-like behavior at near field and give spherical-wave-like at farfield. It is a highly versatile model for coherent optical beams (laser beam) that can be used in many calculation.

z= .; (Subscript[k, 1])=2.*Pi ; w=1;
ContourPlot[Mod[(Subscript[k, 1]) ( z+((  (ρ^2)  z)/(2((z^2)+  ((Subscript[k, 1])^2) (w^4) /4))))
			                             - ArcTan[z, ((-  (Subscript[k, 1])  (w^2))/(2 ))],2*Pi]==Pi,
	  {z,0,5} ,   {ρ,0,5  } , PlotPoints->20]

1.3.2 Compared with spherical wave

Gaussian beam:            (( (w^2) (Subscript[k, 1]))/(i Q[z]))(e^(i  (Subscript[k, 1]) (z+ (( (ρ^2))/(2 Q[z]))) )) where:
                             Q[z]=z-i  (Subscript[ζ, 0]) ; and  (Subscript[ζ, 0]) = (Subscript[k, 1]) (w^(2 ))/2 = ((π (w^(2 )) (Subscript[n, 1 ]))/λ)
Spherical wave: (1/r)(e^(i (Subscript[k, 1])r ))

(Subscript[k, 1])= .; r= .;

(1/r)(E^(I (Subscript[k, 1])r)) /. r->(Sqrt[(ρ^2)+ (z^2)])

((E^(I (Sqrt[(z^2)+(ρ^2)]) (Subscript[k, 1])))/(Sqrt[(z^2)+(ρ^2)])) /. (Sqrt[(z^2)+(ρ^2)])-> z +((ρ^2)/(2 z))

Note that:
Spherical wave in paraxial approx ((e^(i (((ρ^2)/(2 z))+z) (Subscript[k, 1])))/(Sqrt[(z^2)+(ρ^2)])) ≈ (1/z)(e^(i (((ρ^2)/(2 z))+z) (Subscript[k, 1]))) = (e^(i  (Subscript[k, 1])z))(1/z)(e^(i (((Subscript[k, 1]) (ρ^2))/(2 z ))  ))
Gaussian wave                     (( (w^2) (Subscript[k, 1]))/(i Q[z]))(e^(i  (Subscript[k, 1]) (z+ (( (ρ^2))/(2 Q[z]))) ))
If we change z-> Q[z]=z-i  (Subscript[z, 0 ]) for the paraxial spherical wave, we obtain:
                                   (e^(i  (Subscript[k, 1])z))(1/(Q[z]))(e^(i (((Subscript[k, 1]) (ρ^2))/(2 Q[z] ))))
which is the Gaussian beam solution. The term  -i  (Subscript[ζ, 0]) =-i (( π (w^2))/λ)  acts as if it is a virtual point of origin, instead of z=0, it is at -i  (Subscript[z, 0 ]).

From this perspective, Gaussian beam acts like a spherical wavefront originating from a virtual point. It should be remembered that Gaussian beam is only an approximated solution. We obtain it by solving the exact Helmholtz equation. Then, we can also ask: is there an approximated Helmholtz equation that gives Gaussian beam as an exact solution?

1.3.3 Paraxial Helmholtz equation

(∂_(z,z))U +(((∂_ρ) U)/ρ)+(∂_(ρ,ρ))U+ ((Subscript[k, 0])^2)ϵ U==0

(∂_(z,z))((E^(I k z))A[ρ,z ]) +(((∂_ρ) ((E^(I k z))A[ρ,z ]))/ρ)+(∂_(ρ,ρ))((E^(I k z))A[ρ,z ])+ ((Subscript[k, 0])^2)ϵ ((E^(I k z))A[ρ,z ])==0

FullSimplify[-(E^(I k z)) A(ρ,z) (k^2)+2 I (E^(I k z)) Derivative[0,1][A](ρ,z) k+(E^(I k z)) ϵ A(ρ,z) (Subscript[k, 0]^2)+(E^(I k z)) Derivative[0,2][A](ρ,z)+(((E^(I k z)) Derivative[1,0][A](ρ,z))/ρ)+(E^(I k z)) Derivative[2,0][A](ρ,z)]

Derivative[1,0][A](ρ,z)+ρ (A(ρ,z) (ϵ (Subscript[k, 0]^2)-(k^2))+2 I k Derivative[0,1][A](ρ,z)+Derivative[0,2][A](ρ,z)+Derivative[2,0][A](ρ,z))==0

Remember that we have the equation: 
Derivative[1,0][A](ρ,z)+ρ (A(ρ,z) (ϵ (Subscript[k, 0]^2)-(k^2))
					 +2 i k Derivative[0,1][A](ρ,z)+Derivative[0,2][A](ρ,z)+Derivative[2,0][A](ρ,z))==0 
Now, we will let ϵ (Subscript[k, 0]^2)-(k^2)=0:
Derivative[1,0][A](ρ,z)+ρ ( 2 i k Derivative[0,1][A](ρ,z)+Derivative[0,2][A](ρ,z)+Derivative[2,0][A](ρ,z))==0
(((∂^2) A)/(∂ (ρ^2)))+ (1/ρ)((∂ A)/(∂ ρ)) +2 i k ((∂ A)/(∂ z)) + (((∂^2) A)/(∂ (z^2)))==0 .
A[ρ,z] is the envelop function superimposed on the plane wave (e^(i (Subscript[k, 1]) z)). We assume that the wave is very close to plane wave, then the function A[ρ,z] is slowly varying with respect to z. If so, we can drop (((∂^2) A)/(∂ (z^2))) and obtain:
(((∂^2) A)/(∂ (ρ^2)))+ (1/ρ)((∂ A)/(∂ ρ)) +2 i k ((∂ A)/(∂ z))==0. Here we have reduced the equation to the first order of z.
 (((∂^2) A)/(∂ (ρ^2)))+ (1/ρ)((∂ A)/(∂ ρ)) = - 2 i k ((∂ A)/(∂ z)). Now we can solve by part: 
 R[z] ((((∂^2) S)/(∂ (ρ^2)))+ (1/ρ)((∂ S)/(∂ ρ)) )= -S[ρ] 2 i k ((∂ R)/(∂ z)) or:
 (1/(S[ρ]))((((∂^2) S)/(∂ (ρ^2)))+ (1/ρ)((∂ S)/(∂ ρ)) )= - 2 i k (1/(R[z]))((∂ R)/(∂ z))= constant

DSolve[ - 2 I k (1/(R[z]))(∂_z) R[z] == -(K^2) ,R[z],z]

DSolve[(ρ^2) (∂_(ρ,ρ))S[ρ] + ρ(∂_ρ) S[ρ] + (K^2) (ρ^2)S[ρ]==0,S[ρ],ρ]

If K is real positive, the general solution is:
(∫_0^∞) B[K] (e^(-((i (K^2) z)/(2 k))))(Subscript[J, 0])(K ρ)dK
Choose: B[K]= (w^2) K (e^(-(w^2) (K^2)/2))

(∫_0^∞)(w^2) K (E^(-(w^2) (K^2)/2))(E^(-((I (K^2) z)/(2 k))))BesselJ[0, K ρ]\[DifferentialD]K

We have the exact solution: (((e^((k (ρ^2))/(-2 k (w^2)-2 i z))) k (w^2))/(k (w^2)+i z))= ((-i k (w^2))/(Q[z])) (e^(-i k (z - (( (ρ^2))/(2 Q[z]))) ))

Summary

Gaussian beam acts like a the paraxial approximation of a spherical wavefront: (e^(i  (Subscript[k, 1])z))(1/(Subscript[z, V]))(e^(i (((Subscript[k, 1]) (ρ^2))/(2 (Subscript[z, V])))  )), but originating from a virtual point instead of a real point: (Subscript[z, V])->Q[z]=z-i  (Subscript[z, 0 ]), to give:      (e^(i  (Subscript[k, 1])z))(1/(Q[z]))(e^(i (((Subscript[k, 1]) (ρ^2))/(2 Q[z] )))).

Gaussian beam is an exact solution to the paraxial approximation of the Helmholtz Eq. and it is an approximated solution to the exact Helmholtz Eq.

2. Important properties of Gaussian beam

2.1 Beam description

From the above (we drop the subscript 1 from (Subscript[k, 1]) for simplicity):
                                 E[z,ρ]=(1/(Sqrt[2 π]))(( w k)/(i Q[z]))(e^(i k (z+((  (ρ^2))/(2 Q[z]))) ))                    (2.1.1a)
where:                      Q[z]=z-i  (Subscript[ζ, 0])      ;     (Subscript[ζ, 0]) =((k (w^(2 )))/2) = ((π (w^(2 )))/λ)         (2.1.1b)

Note: we have E[z,ρ] in the unit of (1/length) so that the normalized profile:
                                         (∫_0^(2 π))(∫_0^∞)((|E[z,ρ]|)^2)ρdρdθ= 1  is unitless.
Real E has real physical unit, of course. The formula is only for mathematical convenience. We can also drop the normalization factor (1/(Sqrt[2 π])) for convenience. It doesn’t have any special meaning.

We can also write:
                       E[z,ρ]=(( w k)/(Abs[Q[z]] i))(e^(-k (ρ^2) Im[(1/(2Q[z]))]  )) (e^(i  ( k (z+ (ρ^2) Re[(1/(2 Q[z]))])- Arg[Q[z]])  ))  (2.1.2a)
where:                         Q[z]=z-i  (Subscript[ζ, 0]); (Subscript[ζ, 0]) =((k (w^(2 )))/2) = ((π (w^(2 )))/λ)                             
                                     Re[(1/(z-i  (Subscript[ζ, 0])))] = (z/((z^2)+  ((Subscript[ζ, 0])^2) )) ; Im[(1/(z-i  (Subscript[ζ, 0])))] = ((Subscript[ζ, 0])/((z^2)+  ((Subscript[ζ, 0])^2)))                 (2.1.2b)
and:                                       Arg[Q[z]]=-ArcTan[((Subscript[ζ, 0])/z)] ;                                  (2.1.2c)

What are the meaning of all these formulas? What can we infer?

Source code Demo basic Gaussian density plot

Manipulate[
If[RUN== 1,
ω=2 π ;
Manipulate[
Animate[DensityPlot[Re[EGauss[{z,ρ,t},{λ,w,ω}]],{z,-5,5}
,{ρ,-3,3}, PlotPoints-> {100,40}, PlotRange-> All
,ColorFunction->"Rainbow"
,AspectRatio->0.6
,ImageSize->{500,300}
]
,{ {t,0,Button["time",BaseStyle->{16, Bold,FontFamily-> "Calibri"}
,Background->Yellow]},0,1,0.01}
,DefaultDuration->10, AnimationRunning->False]
,{{λ,1,Button["wavelength λ",BaseStyle->{16, Bold,White,FontFamily-> "Calibri"},Background->Hue[0.036,1,0.888,1]]},0.5,1.5}
,{{w,1.,Button["Beamwaist w",BaseStyle->{16, Bold,White,FontFamily-> "Calibri"},Background->Hue[0.344,1.,0.634,1]]},0.5,3}
,SynchronousUpdating->False ]

,TextCell["Not running"
,Black,20,Italic,FontFamily->"Calibri",Background->GrayLevel[0.5]]]
,Row[{Button["Toggle box to start, and turn off when not running",BaseStyle->{16,White, Bold,FontFamily-> "Calibri"}
,Background->Black],Control[{{RUN,0,"->"},{0,1}}]
}]
]

Demo basic Gaussian density plot

Source code Demo basic Gaussian 3D plot

Manipulate[
If[RUN== 1,
ω=2 π ;
Manipulate[
Animate[Plot3D[Re[EGauss[{z,ρ,t},{λ,w,ω}]],{z,-5,5}
,{ρ,-3,3}, PlotPoints-> {100,40}, PlotRange-> All
,Mesh-> False, BoxRatios->{5,3,2}
]
,{ {t,0,Button["time",BaseStyle->{16, Bold,FontFamily-> "Calibri"}
,Background->Yellow]},0,1,0.01}
,DefaultDuration->10, AnimationRunning->False]
,{{λ,1,Button["wavelength λ",BaseStyle->{16, Bold,White,FontFamily-> "Calibri"},Background->Hue[0.036,1,0.888,1]]},0.5,1.5}
,{{w,1.,Button["Beamwaist w",BaseStyle->{16, Bold,White,FontFamily-> "Calibri"},Background->Hue[0.344,1.,0.634,1]]},0.5,3}
,SynchronousUpdating->False ]

,TextCell["Not running"
,Black,20,Italic,FontFamily->"Calibri",Background->GrayLevel[0.5]]]
,Row[{Button["Toggle box to start, and turn off when not running",BaseStyle->{16,White, Bold,FontFamily-> "Calibri"}
,Background->Black],Control[{{RUN,0,"->"},{0,1}}]
}]
]

Demo Demo basic Gaussian 3D plot

2.2 Beam intensity

2.2.1 Transverse profile

For beam intensity:
               I ≺ ((|E|)^2)=((Abs[(( w k)/(Abs[Q[z]] i))(e^(-k (ρ^2) Im[(1/(2Q[z]))]  )) (e^(i  ( k (z+ (ρ^2) Re[(1/(2 Q[z]))])- Arg[Q[z]])  ))])^2)      (2.2.1)
                         I ≺ (( (w^2) (k^2))/((Abs[Q[z]])^(2 )))(e^(-2 k (ρ^2) Im[(1/(2 Q[z]))]  ))= (( (w^2) (k^2))/((Abs[Q[z]])^(2 )))(e^(- ((2 (ρ^2))/((W[z])^2))    ))          (2.2.2)

It is a Gaussian profile:                       (( (w^2) (k^2))/((Abs[Q[z]])^(2 )))(e^(- ((2 (ρ^2))/((W[z])^2))    ))
with a beam intensity radius:  
                      ((W[z])^2) =(1/(k Im[(1/(2 Q[z]))]  ))==>    Im[(1/( Q[z]))] =(2/(k ((W[z])^2))) = (λ/(π ((W[z])^2)))        (2.2.3) 
Substitute: 
                                  ((W[z] )^2)=2 (((z^2)+  ((Subscript[ζ, 0])^2))/( k (Subscript[ζ, 0]))) = ((2 (Subscript[ζ, 0]))/( (Subscript[k,  ])))(1+((z^2)/((Subscript[ζ, 0])^2)))                             (2.2.4)
And with                                              (Subscript[ζ, 0]) =((k (w^(2 )))/2)
We obtain:                      ((W[z] )^2)= (w^2) (1+((z^2)/((Subscript[ζ, 0])^2))) = ((W[0])^2)(1+((z^2)/((Subscript[ζ, 0])^2)))              (2.2.5)

The transverse profile of a Gaussian beam is a Gaussian profile everywhere, with the radius changing as W[0](Sqrt[1+((z^2)/((Subscript[ζ, 0])^2))]). We also refer to beam “spot size” as 2 times the beam radius (beam diameter).

Gaussian beam does not originate from a point: there is no such thing as a point source of light, but there is a place where the Gaussian beam has the smallest profile. The smallest beam radius is called beam waist: W[0]. (see note above also that shows that the point source of a Gaussian beam is a complex number)

HW: plot beam intensity

2.2.2 Beam divergence (focusing)

Obviously, at large z, the beam radius diverges as:
                                                         W[0](Sqrt[1+((z^2)/((Subscript[ζ, 0])^2))]) ~((w z)/(Subscript[ζ, 0]))                  (2.2.6)
Thus, one can define the divergence angle (Subscript[θ, 0]) with:
                                             Tan[(Subscript[θ, 0])] ~ Sin[(Subscript[θ, 0])]= (((w z)/(Subscript[ζ, 0]))/z)=(w/(Subscript[ζ, 0]))                        (2.2.7a)
When we deal with small angle:  Tan[(Subscript[θ, 0])] ~ Sin[(Subscript[θ, 0])]~(Subscript[θ, 0])         
Also, with (Subscript[ζ, 0])-> ((π (w^(2 )))/λ)                                                    
                                              Tan[(Subscript[θ, 0])] ~ Sin[(Subscript[θ, 0])]= (λ/(π w))= (2/π) (( λ)/( 2 w))= (2/π) (( λ)/( D))      (2.2.7b).
where D is the “diameter” of the beam. As it turns out, this formula with a slight modification is almost generic for any beam:
       Tan[Divergence angle] =  Tan[(Subscript[θ, 0])]= m (( λ)/( D))     (2.2.7c)
where m is some number that is a characteristic of a specific beam profile, and which is (2/π) for the case of Gaussian. It is 1.22 for example for a top hat beam. (there is no hard definition for the divergence angle either, so the number m is not something fundamental).

Source code

Manipulate[
If[RUN== 1,
Manipulate[k=2 π/λ  ;  (Subscript[ζ, 0])=0.5 (k^2) (w^2) ;W2=(w^2) (1+((z^2)/((Subscript[ζ, 0])^2))) ; 
{Show[Plot[{-w (Sqrt[1+((v^2)/((Subscript[ζ, 0])^2))]),w (Sqrt[1+((v^2)/((Subscript[ζ, 0])^2))])},{v,-50,50},AspectRatio->0.25,PlotRange->{-5,5},Filling->{1->{{2}, Hue[0.42,1,1,0.7]}}]
,Graphics[{Red,Thickness[0.005],Line[{{z,-5},{z,5}}]}] , ImageSize-> {500,175}, PlotRange-> {{-50,50},{-5,5}}]
,DensityPlot[(( (w^2) (k^2))/((z^2)+((Subscript[ζ, 0])^2)))(E^(- ((2 ((x^2)+(y^2)))/W2)    )),{x,-5,5},{y,-5,5}, Mesh-> False, ColorFunction->GrayLevel,PlotRange-> {0,(( (w^2) (k^2))/((z^2)+((Subscript[ζ, 0])^2)))}, PlotPoints->{50,50}]}
,{z,-50,50}
,{{λ,0.8,Button["wavelength λ",BaseStyle->{16, Bold,White,FontFamily-> "Calibri"},Background->Hue[0.036,1,0.888,1]]},0.3,1.5}
,{{w,1,Button["Beamwaist w",BaseStyle->{16, Bold,White,FontFamily-> "Calibri"},Background->Hue[0.344,1.,0.634,1]]},0.3,4}
]
,TextCell["Not running"
,Black,20,Italic,FontFamily->"Calibri",Background->GrayLevel[0.5]]]
,Row[{Button["Toggle box to start, and turn off when not running",BaseStyle->{16,White, Bold,FontFamily-> "Calibri"}
,Background->Black],Control[{{RUN,0,"->"},{0,1}}]
}]
]

Demo beam divergence 1

For small angle:                                              θ≈ (2/π) (( λ)/( D))
Often, we talk about beam full width:        2 θ≈ (4/π) (( λ)/( D ))        (2.2.8)

An example for non-Gaussian beam: for top hat beam, instead of (4/π)~ 1.27, it is 2.44.
Let’s think:  small waist -> large divergence
                    large waist -> small divergence
As seen above:     waist x divergence = constant.
This gives rise to the concept of etendue and related to the definition of brightness.

A similar formula can be used for focusing (see more in section 4).

2.2.3 Longitudinal profile

Again, for the beam intensity above:
               I ≺ ((|E|)^2)=((Abs[(( w k)/(Abs[Q[z]] i))(e^(-k (ρ^2) Im[(1/(2Q[z]))]  )) (e^(i  ( k (z+ (ρ^2) Re[(1/(2 Q[z]))])- Arg[Q[z]])  ))])^2)      (2.2.1)
                         I ≺ (( (w^2) (k^2))/((Abs[Q[z]])^(2 )))(e^(-2 k (ρ^2) Im[(1/(2 Q[z]))]  ))= (( (w^2) (k^2))/((Abs[Q[z]])^(2 )))(e^(- ((2 (ρ^2))/((W[z])^2))    ))          (2.2.2) 
But now, we are interested in variation along z.

At  ρ=0  the center of the beam
                 (( (w^2) (k^2))/((Abs[Q[z]])^(2 )))= (( (w^2) (k^2))/((z^2)+  ((Subscript[ζ, 0])^(2 )))) =(( (w^2) (k^2))/( ((Subscript[ζ, 0])^(2 ))))(( 1)/((1+((z^2)/((Subscript[ζ, 0])^2))) ))= I[0] (( 1)/((1+((z^2)/((Subscript[ζ, 0])^2))) ))                (2.2.9)
As z-> ∞, the intensity drops as (( 1)/(z^2)) as expected for any spherical wave

As the transverse profile of a Gaussian beam gets larger (or smaller), the intensity drops (or increases) as (( 1)/((1+((z^2)/((Subscript[ζ, 0])^2))) )) (think Lorentzian) by virtue of power conservation.

Source code

Manipulate[
If[RUN== 1,
Manipulate[k=2 π /λ ;  (Subscript[ζ, 0])=0.5 (k^2) (w^2) ;W2=(w^2) (1+((z^2)/((Subscript[ζ, 0])^2))) ; 
Column[{Show[Plot[{-w (Sqrt[1+((v^2)/((Subscript[ζ, 0])^2))]),w (Sqrt[1+((v^2)/((Subscript[ζ, 0])^2))])},{v,-50,50},AspectRatio->0.25,PlotRange->{-5,5},Filling->{1->{{2}, Hue[0.42,1,1,0.7]}}]
,Graphics[{Red,Thickness[0.005],Line[{{z,-5},{z,5}}]}] , ImageSize-> {500,175}, PlotRange-> {{-50,50},{-5,5}}]
,Grid[{{DensityPlot[(( (w^2) (k^2))/((z^2)+((Subscript[ζ, 0])^2)))(E^(- ((2 ((x^2)+(y^2)))/W2)    )),{x,-5,5},{y,-5,5}, Mesh-> False, ColorFunction->GrayLevel,PlotRange-> {0,(( (w^2) (k^2))/((Subscript[ζ, 0])^2))}, PlotPoints->{50,50}]
,Plot3D[(( (w^2) (k^2))/((z^2)+((Subscript[ζ, 0])^2)))(E^(- ((2 ((x^2)+(y^2)))/W2)    )),{x,-5,5},{y,-5,5}, Mesh-> False, PlotRange-> {0,(( (w^2) (k^2))/((Subscript[ζ, 0])^2))}, PlotPoints->{50,50}
,BoxRatios->{1,1,1}]}}] }]
,{z,-50,50}
,{{λ,0.8,Button["wavelength λ",BaseStyle->{16, Bold,White,FontFamily-> "Calibri"},Background->Hue[0.036,1,0.888,1]]},0.3,1.5}
,{{w,1,Button["Beamwaist w",BaseStyle->{16, Bold,White,FontFamily-> "Calibri"},Background->Hue[0.344,1.,0.634,1]]},0.3,4}
]
,TextCell["Not running"
,Black,20,Italic,FontFamily->"Calibri",Background->GrayLevel[0.5]]]
,Row[{Button["Toggle box to start, and turn off when not running",BaseStyle->{16,White, Bold,FontFamily-> "Calibri"}
,Background->Black],Control[{{RUN,0,"->"},{0,1}}]
}]
]

Plot[(( 1)/((1+(z^2)) )),{z,-10,10}, PlotRange-> All, AspectRatio->0.2]

Demo beam intensity along the propagation axis

2.2.4 Power conservation

We now can express the beam intensity
                                      I[0] (( 1)/((1+((z^2)/((Subscript[ζ, 0])^2))) ))(e^(- ((2 (ρ^2))/((W[z])^2))  ))              (2.2.10)
where I[0] is the center of the beam intensity at its waist. (the highest intensity or peak intensity of the beam),

Power is obtained by integrating over all area:

(∫_0^∞)(Subscript[Int, 0]) (( 1)/((1+((z^2)/((Subscript[ζ, 0])^2))) ))(E^(- ((2 (ρ^2))/((W[z])^2))  ))2 π ρ \[DifferentialD]ρ

FullSimplify[((π (Subscript[Int, 0]) ((W(z))^2))/(((2 (z^2))/(Subscript[z, 0]^2))+2))   /. ((W(z))^2)->(w^2) (1+(( (z^2))/(Subscript[z, 0]^2)))]

Power is:                                P= ((I[0])/2)π (w^2)  =((I[0])/2) A              (2.2.11a)
where:                                            A=π (w^2)                             (2.2.11b)
is the beam "area".
We may ask why the factor (1/2)? It is because I[0] is the peak intensity at the center of the beam waist:
                                                    I[0] = ((2 P)/(π (w^2) ))                        (2.2.12)
The power is an integration over all intensity of the whole beam, not just the peak. Hence, we can think of (2.2.11a) as:                                   P= (Subscript[I, ave]) A              (2.2.13) 
with average intensity (Subscript[I, ave])=((I[0])/2) . Think like the area of a triangle. The area is 1/2 of the height times the base.

Sometime we write the beam intensity as:
                                    I[ρ,z]=((2 P)/(π (w^2))) (( 1)/((1+((z^2)/((Subscript[ζ, 0])^2))) ))(e^(- ((2 (ρ^2))/((W[z])^2))  ))or   ((2 P)/(π ((W[z])^2))) (e^(- ((2 (ρ^2))/((W[z])^2))    ))        (2.2.14)

2.3 Beam phase

2.3.1 Phase front

Let's look at the beam again:
                       E[z,ρ]=(( w k)/(Abs[Q[z]] i))(e^(-k (ρ^2) Im[(1/(2Q[z]))]  )) (e^(i  ( k (z+ (ρ^2) Re[(1/(2 Q[z]))])- Arg[Q[z]])  ))  (2.1.2a)
where:                         Q[z]=z-i  (Subscript[ζ, 0]); (Subscript[ζ, 0]) =((k (w^(2 )))/2) = ((π (w^(2 )))/λ)                             
                                     Re[(1/(z-i  (Subscript[ζ, 0])))] = (z/((z^2)+  ((Subscript[ζ, 0])^2) )) ; Im[(1/(z-i  (Subscript[ζ, 0])))] = ((Subscript[ζ, 0])/((z^2)+  ((Subscript[ζ, 0])^2)))                 (2.1.2b)
and:                              -Arg[Q[z]]=-Arg[z-i  (Subscript[ζ, 0])]=ArcTan[((Subscript[ζ, 0])/z)] ;                 (2.1.2c)
The phase term is:
                               k (z+ (ρ^2) Re[(1/(2 Q[z]))])- Arg[Q[z]]              (2.3.1)
The phase front is defined as:
                          k (z+ (ρ^2) Re[(1/(2 Q[z]))])- Arg[Q[z]] =constant.    (2.3.2)
This surface is a function of z and ρ.
What are the meanings of each term?  k z is the plane wave term.  k  (ρ^2) Re[(1/(2 Q[z]))] clearly describes the curvature of the wavefront. Arg[Q[z]]=-ArcTan[((Subscript[ζ, 0])/z)] is indepent of ρ, is just a shift that is dependent on z.

If we have a plane wave traveling along a Gaussian wave of the same  (Subscript[k, 1]), do they traveling in step on the axis?

2.3.2 Guoy effect

Not quite. Because for plane wave, it is  (Subscript[k, 1])z ; for Gaussian, it is  (Subscript[k, 1])z + ArcTan[((Subscript[ζ, 0])/z)]  . So as z goes from -∞ to ∞, The Gaussian differs by a total π phase shift as shown in plot below. This is called Guoy's effect.

Plot[If[z<0,0.5+ArcTan[1/z]/π,ArcTan[1/z]/π -0.5],{z,-10,10}, PlotRange->{-0.5,0.5}, PlotStyle->{ Thick,Blue}
,Filling-> {{1-> {Axis,Hue[0.5,0.7,1,0.6]}}}
, Frame-> True, GridLines-> Automatic
,FrameLabel->{"z (unit of Rayleigh range)","Phase (π)"}
,LabelStyle->{18,FontFamily-> "Arial"}]

2.3.3 Phase front curvature

Phase front term:
                               k (z+ (ρ^2) Re[(1/(2 Q[z]))])- Arg[Q[z]]              (2.3.1)
 At large z, where z>>ρ, the variation of the term is ~
                                         z+ (ρ^2) Re[(1/(2 Q[z]))] = constant       (2.3.3)
But this is the equation for a paraboloidal surface: 
                                    z+  ((ρ^2)/(2 R)) = constant                            (2.3.4)
where R is the radius of curvature.                                        
Therefore, we can define:              Re[(1/(Q[z]))] = (1/(R[z] ))               (2.3.5a)
where R[z] is the radius of curvature of the phase front.
From (2.1.2b)                 Re[(1/(Q[z]))]=Re[(1/(z-i  (Subscript[ζ, 0])))] = (z/((z^2)+  ((Subscript[ζ, 0])^2) ))       (2.1.2b)
                                             R[z]=z+((((Subscript[ζ, 0])^2) )/z)                      (2.3.5b)   
Now we find out what the meaning of (1/(Q[z])):                
                                  Re[(1/(Q[z]))] = (1/(R[z]));                            (2.3.5a)
and                           Im[(1/( Q[z]))] =(2/(k ((W[z])^2))) = (λ/(π ((W[z])^2))) ;       (2.3.6)
Thus, we can write:             (1/(Q[z])) = (1/(R[z]))+ i (λ/(π ((W[z])^2) )):       (2.3.7)

Remember the comparizon of paraxial spherical wave vs Gaussian: (1/z) in spherical becomes (1/(Q[z])) in Gaussian. Since (1/z) in spherical is the radius of curvature, it makes sense that the Re part of (1/(Q[z])) is also the radius of curvature. Remember: (λ/(π ((W[z])^2) )) is NOT the Rayleigh range of the beam. We reserve that term only for (λ/(π ((W[0])^2) )), where W[0] is the beam waist.

Source code demo phase front

Manipulate[
If[RUN== 1,
Manipulate[
 k=2 Pi/λ ;  z0=k  (w^2)/2 ;
planewaveplot=Graphics[{Dashed,Hue[0.35,1,0.7]
, Table[Line[{{(i +0.25)λ,-5},{(i +0.25) λ,5}}],{i,-3,Floor[10./λ]}] }];
phasefrontplot=ContourPlot[Mod[k ( z+(( (x^2)  z)/(2((z^2)+ (z0^2)) )))-Arg[z-I  z0]    ,2*Pi]== Pi
	  , {z,-2,10} ,   {x,-5,5  }, ImageSize->{200,200} , PlotPoints->50, ContourStyle-> {Red}];
profileplot=Plot[{-w (Sqrt[1+((z^2)/(z0^2))]),w (Sqrt[1+((z^2)/(z0^2))])},{z,-2,10}, ImageSize->{200,200},
		PlotRange->{-5,5}, PlotStyle-> RGBColor[0.,0.,1] ] ;
Show[ phasefrontplot,profileplot,planewaveplot] 
,
TextCell["This will take time - please wait"
,Black,20,Italic,FontFamily->"Calibri",Background->Hue[0.05,1,1,0.5]]
,{{λ,1,Button["wavelength λ",BaseStyle->{16, Bold,White,FontFamily-> "Calibri"},Background->Hue[0.036,1,0.888,1]]},0.5,1}
,{{w,1,Button["Beamwaist w",BaseStyle->{16, Bold,White,FontFamily-> "Calibri"},Background->Hue[0.344,1.,0.634,1]]},0.5,3}
,SynchronousUpdating->False ]

,TextCell["Not running"
,Black,20,Italic,FontFamily->"Calibri",Background->GrayLevel[0.5]]]
,Row[{Button["Toggle box to start, and turn off when not running",BaseStyle->{16,White, Bold,FontFamily-> "Calibri"}
,Background->Black],Control[{{RUN,0,"->"},{0,1}}]
}]
]

Demo phase front contour (very slow)

The (1/(Q[z])) term, where  Q[z]=z-i  (Subscript[ζ, 0]), tells us the beam radius of curvature in the real part, and the beam transverse profile radius in the imaginary part.

3. Practical techniques to measure Gaussian (or non-Gaussian) beam

3.1 Intensity profiling: beam scan method

How do we measure the profile of a beam?
Why don't we use things like scanning slit?
Scanning of a Gaussian beam

(2/(π (W^(2 ))))(∫_(-∞)^u)(E^(-((2 ((x^2)+(y^2)))/(W^2))))\[DifferentialD]x

(∫_(-∞)^∞) (E^(-((2 (y^2))/(W^2)))) \[DifferentialD]y

(1/(π (W^2)))(Sqrt[(π/2)]) (Sqrt[(W^2)])(Sqrt[(π/2)]) (Sqrt[(W^2)]) (1-erf((Sqrt[2]) (Sqrt[((u^2)/(W^2))])))

Plot[ (1/2) (1-erf((Sqrt[2]) (u/W)))  /. W->1,
		       {u,-2.5,2.5}, GridLines->Automatic, Frame->True, AxesOrigin->{0,0}
,PlotRange->{{-2.5,2.5},{0,1}}]

N[{(1/2) (1-erf((Sqrt[2]) (u/W))) /. u->W/2,(1/2) (1-erf((Sqrt[2]) (u/W))) /. u->-W/2,(1/2) (1-erf((Sqrt[2]) (u/W))) /. u->-W,(1/2) (1-erf((Sqrt[2]) (u/W))) /. u->W}]

FindRoot[(1/2) (1-erf((Sqrt[2]) x))==0.05,{x,0.5}]

1/x /.x->0.822426686269132822`

3.2 Beam divergence (not at thin lens focal plane)

How do we determine the divergence of a beam? Not all laser beams are Gaussian (remember that). But first, what is “divergence”? What is “far field” (FF). Let’s discuss concepts.

The distribution of a beam intensity as a function of angle at infinity is defined as “far field” (FF). We define a function: f[θ,ϕ]. Sometime, it is written: f[Ω], where Ω={θ,ϕ}.
If we look at a Gaussian beam, the results we saw above is that the beam profile is 
(e^(-((2 (ρ^2))/(((Subscript[W, 0])(1+(z/(Subscript[ζ, 0]))))^2)))). We need to convert this into angle {θ,ϕ} as z->∞.  
(e^(-((2 (ρ^2))/(((Subscript[W, 0])(1+(z/(Subscript[ζ, 0]))))^2)))) -> (e^(-(( ((Subscript[ζ, 0])^2) 2 (ρ^2))/((Subscript[W, 0])(z^2))))) =  (e^(-(( 2 ((Subscript[ζ, 0])^(2 )))/(Subscript[W, 0]))  ((ρ^2)/(z^2))))=(e^(-(( 2 ((Subscript[ζ, 0])^(2 )))/(Subscript[W, 0]))((Sin[θ])^2)))  (also use ((Tan[θ])^2))
We define (( ((Subscript[ζ, 0])^(2 )))/(Subscript[W, 0])) = (1/((Sin[(Subscript[θ, 0])])^2)). Thus: (e^(-2 (((Sin[θ])^(2 ))/((Sin[(Subscript[θ, 0])])^2)))): this is the far field profile function
Remember that Sin[(Subscript[θ, 0])] =(( (Subscript[W, 0 ]))/(Subscript[ζ, 0]))=  (( (Subscript[W, 0 ]))/((π ((Subscript[W, 0])^(2 )))/λ)) = (λ/(π (Subscript[W, 0])))  . 
So now we have a practical question: how do we determine a beam far-field?
Go out very far, measure the spatial profile, and divide by distance z? how far? 1 cm? 1 m? 1 km? The key condition is  z >> (Subscript[ζ, 0]): Rayleigh range.
Rayleigh range is the range where the beam diverges very little and the wavefront is nearly planar: it is where the plane wave approximation is good.

Suppose we have a Gaussian beam with 1-cm diameter (that is, beam waist=0.5 cm), what is the beam divergence? What is its Rayleigh range?

We need to know the wavelength also. Let's λ=0.6328 μm; the beam Rayleigh range is:

((π (0.5 cm)^2)/(0.6328 μm)) /. μm->10^(-4) cm

Suppose we want to measure at a location ~ 10 times the Rayleigh range, it'll be:
1241 m or 1.2 km! How big is the beam there?

2 (Subscript[w, 0]) (Sqrt[1+((z^2)/((Subscript[ζ, 0])^2))]) /. (Subscript[w, 0]) -> 0.5 cm /. ((z^2)/((Subscript[ζ, 0])^2))->100

Clearly, it is not practical to do it this way.

To measure the far field of any beam, we can just focus the beam with a lens (the reason will be clearer in Chapter 9). We can also then measure the profile after focusing at a reasonable distance (do not confuse this with measuring the focal plane of a thin lens that we will learn in Chapter 9).

But then, the divergence we get is just the divergence of the focused beam, not the original beam! How do we convert? We can convert for a Gaussian beam using formula in Section 4 below. But first, let's look at a concept: etendue. 

For a beam or an optical system, the product of accepting area and accepting solid angle is called beam etendue. We can have 1-D etendue or 2-D etendue. Earlier, we ask a question about Gaussian. Let's see:
Beam divergence: Tan[(Subscript[θ, 0])] = (λ/(π (Subscript[W, 0]))) or 2 (Subscript[θ, 0])≈ ((2 λ)/(π (Subscript[W, 0])));
Beam waist: 2 (Subscript[W, 0]); Thus, 1-D etendue = 2 (Subscript[W, 0]) ((2 λ)/(π (Subscript[W, 0])))  = ((4 λ)/(π )): constant!
(What is its 2-D etendue?)
So now, if a lens transforms a Gaussian beam into another Gaussian beam (just of different beam waist), then if we measure the etendue of one, we knows the etendue of the other! Apply this to the laboratory exercise.

3.3 Laboratory exercise: (M^2) and number of times of “diffraction-limited”

Ever wonder when you focus a beam like laser beam, what is at the focus spot? NOT all beams will give the same spot size even with the same focusing convergence angle! This is the concept of “diffraction-limited”. We will study more in Chapter 8, Diffraction.

Perform a Gaussian beam measurement in the lab and report your measurement: locate and measure beam waist, measure divergence angle. Focus the beam with a lens and determine the beam (M^2).

4. Gaussian beam transformation through lenses (and paraboloidal/spherical mirrors)

In practice, rarely any beam is used without any transformation through optics.

Except for integrated optics, which use waveguide to conduct a beam, and the beam profile is determined by the waveguide modes, many systems still use free-space optics.

Even at a very macroscopic level... 
The beam profile properties are essential.

4.1 Review of thin lens and spherical or paraboloidal mirror

Thin lens and low-curvature spherical or paraboloidal mirror has an interesting property: they convert input beam into an output beam with an approximately paraboloidal wavefront:   (e^(-i k ((ρ^2)/(2 f)))) for thin lens where f is the lens focal point and  (e^(i k ((ρ^2)/(Subscript[R, M])))) where (Subscript[R, M]) is the mirror radius of curvature. (derivation comes from paraxial approximation). Convention: f>0 for a positive lens, and (Subscript[R, M]) < 0 for concave, (Subscript[R, M]) > 0 for convex.
Derivation for this comes from ray approximation in which the phase at location ρ is proportional to (ρ^2). (Linked to Simple Optical Elements).
Thus, if a lens is located at location (Subscript[z, 1]), and the input wave is U[x,y,(Subscript[z, 1])], then the output wave at the same location is: U[x,y,(Subscript[z, 1])](e^(-i k ((ρ^2)/(2 f)) )). Whatever happens to the wave afterward? at arbitrary z? this is precisely the interesting property of these optics. In particular, we will see that Gaussian beam behaves quite nicely with this type of optics and that is the reason we study it here.

4.2 Gaussian beam after a thin lens. (derivation)

A Gaussian beam remains a Gaussian through a thin lens (or equivalent parabolic mirror) but with a different beam waist location, value, or direction of propagation (the case of mirror).

Let's set z=0 at the waist of the input beam. Then, let’s write the input Gaussian profile at (Subscript[z, 1]) just before the lens. From: 
                E[z,ρ]=(( w k)/(Abs[Q[z]] i))(e^(-k (ρ^2) Im[(1/(2Q[z]))]  )) (e^(i  ( k (z+ (ρ^2) Re[(1/(2 Q[z]))])- Arg[Q[z]])  ))  (2.1.2a)
to:                E[(Subscript[z, 1]),ρ]≺(e^(-k (ρ^2) Im[(1/(2Q[(Subscript[z, 1])]))]  )) (e^(i  ( k ((Subscript[z, 1])+ (ρ^2) Re[(1/(2 Q[(Subscript[z, 1])]))])- Arg[Q[(Subscript[z, 1])]])  ))            (4.2.1)
(We drop the term (( w k)/(Abs[Q[z]] i)) which is just an overall amplitude factor without any effect on profile and phase for simplicity).
Remember that:       Im[(1/(2 Q[(Subscript[z, 1])]))]=(1/(k ((W[(Subscript[z, 1])])^2))) = (λ/(2 π ((W[(Subscript[z, 1])])^2)))                   (4.2.2a)
and:                           Re[(1/(2 Q[(Subscript[z, 1])]))]= (1/(2R[(Subscript[z, 1])]))                              (4.2.2b)
Then, by substitute (4.2.2), (4.2.1) becomes:
                                 (e^(- ((ρ^2)/((W[(Subscript[z, 1])])^2))  )) (e^(i  ( k ((Subscript[z, 1])+  ((ρ^2)/(2R[(Subscript[z, 1])])))- Arg[Q[(Subscript[z, 1])]])  ))                 (4.2.3)
The output beam just after the lens is:  
                (e^(- ((ρ^2)/((W[(Subscript[z, 1])])^2))  )) (e^(i  ( k ((Subscript[z, 1])+  ((ρ^2)/(2R[(Subscript[z, 1])])))- Arg[Q[(Subscript[z, 1])]])  )) (e^(-i k ((ρ^2)/(2 f))))
                   =(e^(- ((ρ^2)/((W[(Subscript[z, 1])])^2))  ))(e^(i k (Subscript[z, 1 ])))(e^(i k ((  (ρ^2))/(2 ))((( 1)/(R[(Subscript[z, 1])]))-(1/f)) ))(e^(-i Arg[Q[(Subscript[z, 1])]]))   (4.2.4)

If we define:                          ((  1)/(R'[(Subscript[z, 1])]))=((  1)/(R[(Subscript[z, 1])]))-(1/(f )),                                 (4.2.5)
(Note: R'[(Subscript[z, 1])] is NOT the derivative of R[z] at (Subscript[z, 1]). )
then we can write:           (e^(- ((ρ^2)/((W[(Subscript[z, 1])])^2))  ))(e^(i k (Subscript[z, 1 ])))(e^(i k ((  (ρ^2))/(2 R'[(Subscript[z, 1])] ))))(e^(-i Arg[Q[(Subscript[z, 1])]]))    (4.2.6)

Let's look at the expression. What does it tell us?
Obviously, it has a Gaussian amplitude profile: (e^(- ((ρ^2)/((W[(Subscript[z, 1])])^2)))); and a paraboloidal phase term: (e^(i k (Subscript[z, 1 ])))(e^(i k ((  (ρ^2))/(2 R'[(Subscript[z, 1])] )))) . This is a typical Gaussian profile! (note: The phase term (e^(-i Arg[Q[(Subscript[z, 1])]])), like the amplitude term (( w k)/(Abs[Q[z]] i)) that we drop, is ONLY a constant phase term without any ρ-dependency hence, just a factor). 

The only difference between this output Gaussian and the input Gaussian is that the phase front has a radius of curvature  R'[(Subscript[z, 1])] instead of R[(Subscript[z, 1])]. They both have the same beam spot radius at (Subscript[z, 1]): W[(Subscript[z, 1])].  Thus, this is also a solution of the Helmholtz equation.
The output of a Gaussian beam through a an on-axis thin lens is a Gaussian beam. BUT this Gaussian beam is NOT the same as the input (inspite of the fact that BOTH have the same beam spot radius at (Subscript[z, 1]): they have difference phase front at (Subscript[z, 1])).
So, to find out the property of this output Gaussian from the lens, we must find a Gaussian beam expression that gives us exactly the beam spot radius: W[(Subscript[z, 1])] and radius of curvature  R'[(Subscript[z, 1])] at location (Subscript[z, 1]).
Let's write the expression again: 
(e^(- ((ρ^2)/((W[(Subscript[z, 1])])^2))  ))(e^(i k (Subscript[z, 1 ])))(e^(i k ((  (ρ^2))/(2 R'[(Subscript[z, 1])] ))))(e^(-i Arg[Q[(Subscript[z, 1])]]))= (e^((1/2)i k (ρ^2) ((( 1)/( R'[(Subscript[z, 1])]))+i  (2/(k ((W[(Subscript[z, 1])])^2)))) ))(e^(i k (Subscript[z, 1 ])))(e^(-i Arg[Q[(Subscript[z, 1])]]))
Now we recall that in general,              (( 1)/( R'[(Subscript[z, 1])]))+i (2/(k ((W[(Subscript[z, 1])])^2)))=(1/(P[(Subscript[z, 1])] ))    (4.2.7)
where P[(Subscript[z, 1])] is the Q-function (virtual point) for some Gaussian beam.

How do we do that?

Basically we must find the beam waist location (Subscript[z, 2]) and beam waist (Subscript[W, 0])'.  We write the general expression for a Gaussian beam:
                   (e^(- ((ρ^2)/((W[z])^2)))) (e^(i  ( k (z+  ((ρ^2)/(2R'[z])))-Arg[z- i (Subscript[ζ, 0])'])  ))                (4.2.8)
This is expression for a beam with waist at z=0. We shift the coordinate to (Subscript[z, 2]):  
            (e^(- ((ρ^2)/((W'[z-(Subscript[z, 2])])^2))  )) (e^(i  ( k (z-(Subscript[z, 2])+  ((ρ^2)/(2R'[z-(Subscript[z, 2])])))-Arg[z-(Subscript[z, 2])- i (Subscript[ζ, 0])'])  ))       (4.2.9)     
At z=(Subscript[z, 1]): 
         (e^( ((ρ^(2 ))/( ((W'[(Subscript[z, 1])-(Subscript[z, 2])])^2)  )))) (e^(i  ( k ((Subscript[z, 1])-(Subscript[z, 2])+  ((ρ^2)/(2R'[(Subscript[z, 1])-(Subscript[z, 2])])))-Arg[(Subscript[z, 1])-(Subscript[z, 2])- i (Subscript[ζ, 0])'])  ))       (4.2.10) 
Now, this beam must match with the output Gaussian beam from the lens at (Subscript[z, 1]), right?
So, we compare terms by terms:
1- Amplitude profile:              (e^(- ((ρ^2)/((W'[(Subscript[z, 1])-(Subscript[z, 2])])^2))  )) = (e^(-((ρ^2)/( ((W[(Subscript[z, 1])])^2)))    ))           (4.2.11) 
which yields:((W'[(Subscript[z, 1])-(Subscript[z, 2])])^2) =((W[(Subscript[z, 1])])^(2  )); or:
                     (((Subscript[W, 0])')^2)(1+(((((Subscript[z, 1])-(Subscript[z, 2])))^2)/(((Subscript[ζ, 0])')^2)))=((Subscript[W, 0])^2)(1+(((Subscript[z, 1])^2)/((Subscript[ζ, 0])^2)))  
                    (((Subscript[W, 0])')^2)(1+(((((Subscript[z, 1])-(Subscript[z, 2])))^2)/(((π^2) (((Subscript[W, 0])')^(4 )))/(λ^2))))=((Subscript[W, 0])^2)(1+(((Subscript[z, 1])^2)/(((π^2) ((Subscript[W, 0])^(4 )))/(λ^2))))   
                    (((Subscript[W, 0])')^2)+(((λ^2) ((((Subscript[z, 1])-(Subscript[z, 2])))^2))/((π^2) (((Subscript[W, 0])')^2)))=((Subscript[W, 0])^2)+(((λ^2) ((Subscript[z, 1])^2))/((π^2) ((Subscript[W, 0])^2) ))
                    (π/λ)(((Subscript[W, 0])')^2)+((λ ((((Subscript[z, 1])-(Subscript[z, 2])))^2))/(π (((Subscript[W, 0])')^2)))=(π/λ)((Subscript[W, 0])^2)+((λ ((Subscript[z, 1])^2))/(π ((Subscript[W, 0])^2) ))
                      (Subscript[ζ, 0])'+(( ((((Subscript[z, 1])-(Subscript[z, 2])))^2))/((Subscript[ζ, 0])'))=(Subscript[ζ, 0])+(( ((Subscript[z, 1])^2))/(Subscript[ζ, 0 ]))                     (4.2.12)        
2- Phase profile: 
       (e^(i  ( k ((Subscript[z, 1])-(Subscript[z, 2])+  ((ρ^2)/(2R'[(Subscript[z, 1])-(Subscript[z, 2])])))-Arg[(Subscript[z, 1])-(Subscript[z, 2])- i (Subscript[ζ, 0])'])  ))≺ (e^(i k (Subscript[z, 1 ])))(e^(i k ((  (ρ^2))/(2 ))(((  1)/(R[(Subscript[z, 1])]))-(1/f)) ))(e^(-i Arg[(Subscript[z, 1])-i (Subscript[ζ, 0])] ))    (4.2.13)   
Again, separate the ρ-dependent term from the others:  
            (e^(i   k  ((ρ^2)/(2R'[(Subscript[z, 1])-(Subscript[z, 2])] ))))= (e^(i k ((  (ρ^2))/(2 ))(((  1)/(R[(Subscript[z, 1])]))-(1/f)) ))=>((  1)/(R'[(Subscript[z, 1])-(Subscript[z, 2])]))=((  1)/(R[(Subscript[z, 1])]))-(1/f)         (4.2.14)   
which is precisely the wavefront relation we mention earlier.
The other term is NOT important here since it is just a constant phase shift (no ρ or z-dependence). 
Recalling that ((  1)/(R[(Subscript[z, 1])]))= ((Subscript[z, 1])/(((Subscript[z, 1])^2)+  ((Subscript[ζ, 0])^2))), we have the equation:
                                             ((  (Subscript[z, 1])-(Subscript[z, 2]))/(((((Subscript[z, 1])-(Subscript[z, 2])))^2)+  (((Subscript[ζ, 0])')^2)))=((Subscript[z, 1])/(((Subscript[z, 1])^2)+  ((Subscript[ζ, 0])^2)))-(1/(f  ))        (4.2.15)

Solving the equations

Let's rewrite the Eqs. again:
                                          (Subscript[ζ, 0])'+(( ((((Subscript[z, 1])-(Subscript[z, 2])))^2))/((Subscript[ζ, 0])'))=(Subscript[ζ, 0])+(( ((Subscript[z, 1])^2))/(Subscript[ζ, 0]))            (4.2.12)
                                  ((  (Subscript[z, 1])-(Subscript[z, 2]))/(((((Subscript[z, 1])-(Subscript[z, 2])))^2)+  (((Subscript[ζ, 0])')^2)))=((Subscript[z, 1])/(((Subscript[z, 1])^2)+  ((Subscript[ζ, 0])^2)))-(1/(f  ))               (4.2.15)
Changing (Subscript[z, 1])-(Subscript[z, 2]) -> -d
                                          (Subscript[ζ, 0])'+(( (d^2))/((Subscript[ζ, 0])'))=(Subscript[ζ, 0])+(( ((Subscript[z, 1])^2))/(Subscript[ζ, 0]))                  (4.2.16a)
                                          ((  d)/((d^2)+  (((Subscript[ζ, 0])')^2)))=(1/f)-((Subscript[z, 1])/(((Subscript[z, 1])^2)+  ((Subscript[ζ, 0])^2)  ))              (4.2.16b)
So we have 2 eqs. with 2 unknowns.

FullSimplify[Solve[{ ζ0p+(( (d^2))/(ζ0p))==ζ0 +(( ((Subscript[z, 1])^2))/(ζ0))
				,((  d)/((d^2)+   ((ζ0p)^2)))==(1/f) -((Subscript[z, 1])/(((Subscript[z, 1])^2)+  ((ζ0)^2)))},{d,ζ0p}]]

4.3 Result discussion

4.3.1 Relations between the beams

The first result:            (Subscript[ζ, 0])' = (Subscript[ζ, 0]) ((f^(2 ))/(((((Subscript[z, 1])-f))^2)+((Subscript[ζ, 0])^2))) or                   (4.3.1)
                           
                           (Subscript[W, 0])' = (Subscript[W, 0])(Sqrt[((f^(2 ))/(((((Subscript[z, 1])-f))^2)+((Subscript[ζ, 0])^2)))]) = (Subscript[W, 0]) M         (4.3.2)
where we define: M=(Sqrt[((f^(2 ))/(((((Subscript[z, 1])-f))^2)+((Subscript[ζ, 0])^2)))]) =((|f| )/(((Subscript[z, 1])-f)))(Sqrt[(1/(1+(((Subscript[ζ, 0])^2)/((((Subscript[z, 1])-f))^2))))])     (4.3.3)
Next:
                 d->((((Subscript[z, 1])-f) (f^2))/((f^2)-2 (Subscript[z, 1]) f+((ζ0)^2)+(Subscript[z, 1]^2)))+f =>  d -f =((Subscript[z, 1])-f) ((f^2)/(((((Subscript[z, 1])-f))^2)+((ζ0)^2)))      
                                                  d  =f+((Subscript[z, 1])-f) (M^2)           (4.3.4)
The quantity M is interesting: it is the ratio between the 2 beam waists. Thus we define M to be the "magnification" factor (as if we use the lens to magnify or focus the beam). We see these relations are somewhat similar to ray-optics of lens that we are familiar. 
In fact we have other relations involving M:
Rayleigh range: (Subscript[ζ, 0])' = (Subscript[ζ, 0]) ((f^(2 ))/(((((Subscript[z, 1])-f))^2)+((Subscript[ζ, 0])^2))) = (M^2) (Subscript[ζ, 0])             (4.3.5)
Beam divergence:   θ' =(λ/(π (Subscript[W, 0])')) =  (1/M) (λ/(π (Subscript[W, 0])))= (1/M)θ .      (4.3.6)

Write a code calculate output beam waist and location as a function of input beam parameters

Transforming: Given a beam, and you wish to transform it into another beam, can you find a lens to do it (HW)

4.3.2 Meaning of the "magnification" (or reduction) term

Let's look at the term: 
                                                M=(Sqrt[((f^(2 ))/(((((Subscript[z, 1])-f))^2)+((Subscript[ζ, 0])^2)))])              (4.3.3)
If the beam is point-like: it has a small beam waist and short Rayleigh range, then:
  M=(Sqrt[((f^(2 ))/(((((Subscript[z, 1])-f))^2)+((Subscript[ζ, 0])^2)))]) = ((|f )/(((Subscript[z, 1])-f)))|(Sqrt[(1/(1+(((Subscript[ζ, 0])^2)/((((Subscript[z, 1])-f))^2))))]) ≈ |((f )/(((Subscript[z, 1])-f)))|   (4.3.7)
If (Subscript[|z, 1])-f| >>(Subscript[ζ, 0 ]) so that (((Subscript[ζ, 0])^2)/((((Subscript[z, 1])-f))^2)) ≈ 0.
But this is the typical ray-optics lens relation:

Image magnification:   (Subscript[W, 0])' = (Subscript[W, 0]) M  
In fact, we look at other relation:
                                                     d -f =((Subscript[z, 1])-f) (M^2)
or:                                                d -f =((Subscript[z, 1])-f) ((f^2)/((((Subscript[z, 1])-f))^2))
                                                    d -f =((f^2)/(((Subscript[z, 1])-f)))
                                                     (d -f) ((Subscript[z, 1])-f)=(f^2)
                                                     d (Subscript[z, 1]) -f(d+(Subscript[z, 1]))+(f^2)=(f^2)
                                                     (((d+(Subscript[z, 1])))/(d (Subscript[z, 1]))) = (1/f)=> (1/(Subscript[z, 1]))+(1/d)= (1/f)    (4.3.8)
This is the lens imaging relation in geometrical (ray) optics.
Depth of focus:                                   2(Subscript[ζ, 0])' = (M^2) 2(Subscript[ζ, 0 ])            (4.3.9)
All of these relations are compatible with ray optics.
When is this concept not relevant? obviously, when the beam waist is large: the beam is more "plane-wave" like rather than spherical wave emanating from a spot: now we know when we can use simple geometrical optics and when we should not when dealing with Gaussian beam: the Rayleigh ranges << f or the waist is very small.

4.3.3 Gaussian beam focusing or mode matching

Now we have the other extreme: our beam is pencil like, which means it has a large Rayleigh range:  (Subscript[ζ, 0])>> f. Also, the beam waist should also be close to the lens, in other words,  (Subscript[z, 1]) is NOT >> f. What will happen?

Let's look at beam waist: 
(((Subscript[W, 0])' )^2)= ((Subscript[W, 0])^2)((f^(2 ))/(((((Subscript[z, 1])-f))^2)+((Subscript[ζ, 0])^2))) =(((Subscript[W, 0 ])^2)/((Subscript[ζ, 0])^2)) ((f^(2 ))/((((((Subscript[z, 1])-f))^2)/((Subscript[ζ, 0])^2))+1))=(((λ^2) (f^(2 )))/((π^2)((Subscript[W, 0])^2)))(1/((((((Subscript[z, 1])-f))^2)/((Subscript[ζ, 0])^2))+1))     
Since we assume that (Subscript[ζ, 0]) >> |(Subscript[z, 1])-f| , we can write: 
                                                (((Subscript[W, 0])')^2) ≈ (1/((π^2)((Subscript[W, 0])^2)))(λ^2) (f^(2  ))
or the famous Gaussian beam transform formula:

π (Subscript[W, 0]) (Subscript[W, 0])' = λ f                       (4.3.10) 
Notice how the beam waists are inverse of each other: input a large beam, we have a very small focused spot. Input a small spot, we have a large focus spot.

4.3.3 Focal-focal conjugate

Now suppose we want the "symmetric case". Can we make a beam such that the input is similar to the output except for a displacement of the waist?
What if we just put a Gaussian waist at the focal plane of a lens? That is:  (Subscript[z, 1])=f or (Subscript[z, 1])-f=0 ?
From the above:                        d  =f+((Subscript[z, 1])-f) (M^2) = f          
This means that the output beam waist is at the other focal plane. Sure, this is simple enough to understand: if we put a beam waist at one focal plane, the lens just bend the phase such that the output is at the other focal plane. But what happens to beam waist size? For that, we need:
                                                  M=(Sqrt[((f^(2 ))/(((((Subscript[z, 1])-f))^2)+((Subscript[ζ, 0])^2)))]) = (f/(Subscript[ζ, 0]))         
If we want the two beams also to be nearly equal in waist, we need  f ~(Subscript[ζ, 0]), in other words, the Rayleigh range of the beam has to be nearly equal to the lens f. So we can't arbitrary make beams like that for a given lens. Either we choose the lens f or we must have beam with proper waist!

Demo of Gaussian beam through a lens

Demo Homework lens selection for Gaussian beam transform

Consider focusing a beam:

Consider focusing. From the formula:
                                                     π (Subscript[W, 0]) (Subscript[W, 0])' = λ f                       (4.3.10) 
                                                      (Subscript[W, 0])' =(( λ f)/(π (Subscript[W, 0]) ))   
If we want a very small focused spot, we want (Subscript[W, 0]) as large as possible. But the input beam can’t be larger that the lens! otherwise, the power will not be collected. Hence, we can say that (Subscript[W, 0]) is limited by the radius of the lens r. Hence:                               (Subscript[W, 0])' ~(( λ f)/(π r )) 
But the ratio (f/r)~ (1/(Sin[θ])), which is also called nunerical aperture (NA) of the lens. We see that it doesn’t matter large lens or small lens, the only thing that matter is the NA:

The above 2 lenses give exactly the same spot size when focusing:
                                                         2(Subscript[W, 0])' ~2(( λ)/(π NA))            (4.3.11) 
That’s why NA is the only rating that matter in lenses when it comes to resolution: how fine the details we can see or how small a spot we can focus.

4.3.4 Off-axis beam: (advanced HW only)

What if the beam is off the axis?  
(e^(-((ρ^2)/( ((W[(Subscript[z, 1])])^2)))  )) (e^(i  ( k((Subscript[z, 1])+  ((ρ^2)/(2R)))- Arg[Q[(Subscript[z, 1])]])  )) (e^(-i k (((((ρ Cos[θ]+D))^2)+(ρ^(2 ))((Sin[θ])^2))/(2 f)) ))
(((ρ Cos[θ]+D))^2)+(ρ^(2 ))((Sin[θ])^2)=(ρ^(2 )) + (D^2)+2 ρ D Cos[θ]  
(e^(-((ρ^2)/( ((W[(Subscript[z, 1])])^2)))  )) (e^(i  ( k((Subscript[z, 1])+  ((ρ^2)/(2R)))- Arg[Q[(Subscript[z, 1])]])  )) (e^(-i k ((ρ^(2 ))/(2 f)))) (e^(-i k ((D^2)/(2 f))))(e^(-i k ((ρ  D Cos[θ])/f)))
(e^(-((ρ^2)/( ((W[(Subscript[z, 1])])^2)))  ))(e^(i k (Subscript[z, 1 ])))(e^(i k ((  (ρ^2))/(2 ))(((  1)/R)-(1/f)) ))(e^(-i (Arg[Q[(Subscript[z, 1])]]+k ((D^2)/(2 f))))) (e^(-i k ((ρ  D Cos[θ])/f)))
Meaning of the various terms:
(e^(-((ρ^2)/( ((W[(Subscript[z, 1])])^2)))  )) As before, this is the Gaussian amplitude profile
(e^(i k (Subscript[z, 1 ])))(e^(i k ((  (ρ^2))/(2 ))(((  1)/R)-(1/f)) )) This is the plane wave combined with a wavefront with curvature R':
                               (1/(R'))=((  1)/R)-(1/f)
(e^(-i (Arg[Q[(Subscript[z, 1])]]+k ((D^2)/(2 f)))))  There is a NET constant phase shift (e^(-i k ((D^2)/(2 f)))): this is the "beam
                              spatially averaged phase shift" from the lens.
Finally, the interesting term: (e^(-i k ((ρ  D Cos[θ])/f))) : beam position-dependent phase shift. We can just write: (e^(-i k y((D )/f))). This shows that the phase of the beam has been distorted in the x-y plane: along the y axis, the beam experiences a linear phase shift.
Remember: this is only the beam profile at location (Subscript[z, 1]). This is NOT necessarily the solution at any z!

4.3.5 Linear decomposition approach

We can try the linear composition approach.
First: Decompose the beam profile into x and y, forget all about constant phase
(e^(-((((x^2)+(y^2) ))/( ((W[(Subscript[z, 1])])^2)))))(e^(i k ((  ((x^2)+(y^2) ))/(2 (Subscript[R, out]) )))) (e^(-i k y ((D )/f))) 
(∫_□^□)(e^(-((((x^2)+(y^2) ))/( ((W[(Subscript[z, 1])])^2)))))(e^(i k ((  ((x^2)+(y^2) ))/(2 (Subscript[R, out]) )))) (e^(-i k y ((D )/f)))(e^(-i (Subscript[q, 1]) x))(e^(-i (Subscript[q, 2]) y))dxdy 
Let's do some simplification: 
(e^(-((x^2)/( ((W[(Subscript[z, 1])])^2)))))(e^(i k ((  (x^2))/(2 (Subscript[R, out]) )))) = (e^((1/2)i k (x^2)((1/( (Subscript[R, out])))+i(2/( k ((W[(Subscript[z, 1])])^2))))))=(e^(i (k/(2 P[(Subscript[z, 1])])) (x^2))) where we define  (1/P)=((1/( (Subscript[R, out])))+i(2/( k ((W[(Subscript[z, 1])])^2))))

((E^(I (k/(2 P[(Subscript[z, 1])])) (y^2)))(E^(-I k y ((D )/f))))(∫_(-∞)^∞)(E^(I (k/(2 P[(Subscript[z, 1])])) (x^2))) (E^(-I (Subscript[q, 1]) x))\[DifferentialD]x

We get:
((e^(i (k/(2 P[(Subscript[z, 1])])) (y^2)))(e^(-i k y ((D )/f))))  (e^(-((i P((Subscript[z, 1])) (Subscript[q, 1]^2))/(2 k)))) (Sqrt[2 π]) (Sqrt[((i P((Subscript[z, 1])))/k)])
Next:

(∫_(-∞)^∞)(E^(I (k/(2 P[(Subscript[z, 1])])) (y^2)))(E^(-I k y ((D )/f))) (E^(-I (Subscript[q, 2]) y))\[DifferentialD]y

(e^(-((i P((Subscript[z, 1])) (((D k+f (Subscript[q, 2])))^2))/(2 (f^2) k)))) (Sqrt[2 π]) (Sqrt[((i P((Subscript[z, 1])))/k)])=(e^(-((i P((Subscript[z, 1])) )/(2  k))((((Subscript[q, 2])+((D k)/f)))^2))) (Sqrt[2 π]) (Sqrt[((i P((Subscript[z, 1])))/k)]) 
The whole term is:  
(e^(-((i P((Subscript[z, 1])) (Subscript[q, 1]^2))/(2 k)))) (Sqrt[2 π]) (Sqrt[((i P((Subscript[z, 1])))/k)])(e^(-((i P((Subscript[z, 1])) )/(2  k))((((Subscript[q, 2])+((D k)/f)))^2))) (Sqrt[2 π]) (Sqrt[((i P((Subscript[z, 1])))/k)]) 
which is:   
(e^(-((i P((Subscript[z, 1])) (Subscript[q, 1]^2))/(2 k)))) (e^(-((i P((Subscript[z, 1])) )/(2  k))((((Subscript[q, 2])+((D k)/f)))^2)))2 π((i P((Subscript[z, 1])))/k)

To reconstruct the beam linearly:  
((i P((Subscript[z, 1])))/k)(∫_□^□)(e^(-((i P((Subscript[z, 1])) (Subscript[q, 1]^2))/(2 k)))) (e^(-((i P((Subscript[z, 1])) )/(2  k))((((Subscript[q, 2])+((D k)/f)))^2)))(e^(i (Subscript[q, 1])x))(e^(i (Subscript[q, 2])y))(e^(i (Sqrt[(k^2) - ((Subscript[q, 1])^2) -((Subscript[q, 2])^2)])z))d(Subscript[q, 1])d(Subscript[q, 2])
(Sqrt[(k^2) - ((Subscript[q, 1])^2) -((Subscript[q, 2])^2)]) ~ k-(((Subscript[q, 1])^2)/(2 k))-(((Subscript[q, 2])^2)/(2 k))
((i P((Subscript[z, 1])))/k)(e^(i k z))(∫_□^□)(e^(-i ((Subscript[q, 1]^2)/(2 k))( P((Subscript[z, 1]))+z) )) (e^(-((i )/(2  k))(P((Subscript[z, 1])) ((((Subscript[q, 2])+((D k)/f)))^2)+((Subscript[q, 2])^2)z)))(e^(i (Subscript[q, 1])x))(e^(i (Subscript[q, 2])y))d(Subscript[q, 1])d(Subscript[q, 2])

(∫_(-∞)^∞)(E^(-I ((Subscript[q, 1]^2)/(2 k))( P((Subscript[z, 1]))+z) ))(E^(I (Subscript[q, 1])x))\[DifferentialD](Subscript[q, 1])

FullSimplify[(∫_(-∞)^∞)(E^(-((I )/(2  k))(P((Subscript[z, 1])) ((((Subscript[q, 2])+ k Tan[α]))^2)+((Subscript[q, 2])^2)z)))(E^(I (Subscript[q, 2])y))\[DifferentialD](Subscript[q, 2])]

(E^((FullSimplify[(((I P((Subscript[z, 1])) tan(α)-I y))^2)])/(FullSimplify[4 (((I z)/(2 k))+((I P((Subscript[z, 1])))/(2 k)))])))

(e^((i k (((y-P((Subscript[z, 1])) tan(α)))^2))/(2 (z+P((Subscript[z, 1])))))) (Sqrt[2 π]) (Sqrt[-((i k)/(P((Subscript[z, 1]))+z))]) (e^(-(1/2)i k P((Subscript[z, 1])) (tan^2)(α) )) 
The solution: 
((i P((Subscript[z, 1])))/k)(e^(i k z)) (e^((i k (x^2))/(2 (z+P((Subscript[z, 1])))))) (Sqrt[2 π]) (Sqrt[-((i k)/(z+P((Subscript[z, 1]))))])(e^((i k (((y-P((Subscript[z, 1])) tan(α)))^2))/(2 (z+P((Subscript[z, 1])))))) (Sqrt[2 π]) (Sqrt[-((i k)/(P((Subscript[z, 1]))+z))])  (e^(-(1/2)i k P((Subscript[z, 1])) (tan^2)(α))) 
simplify:    
((2 π P((Subscript[z, 1])) )/(z+P((Subscript[z, 1])))) (e^(i k z)) (e^((i k (x^2))/(2 (z+P((Subscript[z, 1]))))))  (e^((i k (((y-P((Subscript[z, 1])) tan(α)))^2))/(2 (z+P((Subscript[z, 1]))))))    (e^(-(1/2)i k P((Subscript[z, 1])) (tan^2)(α)))  
We can also make a transformation of the y axis. We shift by D it so that it is on the lens axis. y-> y-D
If we are interested only in beam intensity: 
((i k (x^2))/(2 (z+P((Subscript[z, 1])))))=((i k (x^2))/2)(1/(z+P[(Subscript[z, 1])] )) and (( 1)/( R'))+i (2/(k ((W[(Subscript[z, 1])])^2)))=(1/(P[(Subscript[z, 1])]))
We find the transformation such that: z+P[(Subscript[z, 1])]= Q[z-(Subscript[z, 2])]

Expand[((( 1)/( R'))+I (2/(k ((W[(Subscript[z, 1])])^2))))(1+(( z)/( R'))-I(2/(k ((W[(Subscript[z, 1])])^2))))]

BeamOut[beamin_, lensParm_,D_, x_,y_,z_]:=Module[
		               {λ,W0, z0, Wz1,Rz1,R21,f,z1,k,P1, Tα },
		                {f,z1}=lensParm; 
		                 {λ, W0 }=beamin ;
		              k=2*Pi/λ ; Tα=D/f  ;
		                  z0=k (W0^2)/2 ; Wz1=W0(Sqrt[1+((((z1/z0)))^2)]) ;
	                   	Rz1=(( z1)/(((z1^2)+ (z0^2)) ));
		                R21=Rz1-1/f  ;
		                P1=1/(R21+I(2/( k (Wz1^2)))) ; 
		             
		      ((2 π P1 )/(z+P1)) (E^(I k z)) (E^((I k (x^2))/(2 (z+P1))))  (E^((I k (((y-D-P1 Tα ))^2))/(2 (z+P1))))   (E^(-(1/2) I k P1 ((Tα)^2))) ] ;

(* BeamIntOut[beamin_, lensParm_,D_, x_,y_,z_]:=Module[
		               {λ,W0, z0, Wz1,Rz1,R21,f,z1,k,P1, Tα },
		                {f,z1}=lensParm; 
		                 {λ, W0 }=beamin ;
		              k=2*Pi/λ ; Tα=D/f  ;
		                  z0=k (W0^2)/2 ; Wz1=W0(Sqrt[1+((((z1/z0)))^2)]) ;
	                   	Rz1=(( z1)/(((z1^2)+ (z0^2)) ));
		                R21=Rz1-1/f  ;
		                P1=((R21-I(2/( k (Wz1^2))))/((R21^2)+(4/( (k^2) (Wz1^4))))) ;
(*  z+((R21-I(2/( k (Wz1^2))))/((R21^2)+(4/( (k^2) (Wz1^4)))))= ((z((R21^2)+(4/( (k^2) (Wz1^4))))+R21-I(2/( k (Wz1^2))))/((R21^2)+(4/( (k^2) (Wz1^4)))))
			          
		             
		     Abs[ ((2 π P1 )/(z+P1))]  (E^((I k (x^2))/(2 (z+P1))))  (E^((I k (((y-P1 Tα ))^2))/(2 (z+P1))))   *) ] *)

lens={1,0}*10^4  ;  (* unit micron  *)
beamin={1.5,0.3*10^4}    (* unit micron  *) ;
Plot3D[Re[BeamOut[beamin, lens,0,0,y,z]],
		         {z,9980,10020},{y,-10,10}
, PlotPoints->{120,100}, BoxRatios->{2,1,0.5}
,Mesh->False, PlotRange->All, ViewPoint->{0,0,10}]

lens={1,0}*10^4  ;  (* unit micron  *)
beamin={1.5,0.3*10^4}    (* unit micron  *) ;
Plot3D[Re[BeamOut[beamin, lens,3000,0,y,z]],
		         {z,9980,10020},{y,-10,10}
, PlotPoints->{120,100}, BoxRatios->{2,1,0.5}
,Mesh->False, PlotRange->All, ViewPoint->{0,0,10}]

lens={1,0}*10^4  ;  (* unit micron  *)
beamin={1.5,0.3*10^4}    (* unit micron  *) ;
Plot3D[Re[BeamOut[beamin, lens,10000,0,y,z]],
		         {z,9980,10020},{y,-10,10}
, PlotPoints->{120,100}, BoxRatios->{2,1,0.5}
,Mesh->False, PlotRange->All, ViewPoint->{0,0,10}]

lens={1,1}*10^4  ;  (* unit micron  *)
beamin={1.5,1.5}    (* unit micron  *) ;
Plot3D[Re[BeamOut[beamin, lens,0,0,y,z]],
		         {z,9980,10020},{y,-10,10}
, PlotPoints->{120,100}, BoxRatios->{2,1,0.5}
,Mesh->False, PlotRange->All, ViewPoint->{0,0,10}]

lens={1,2}*10^4  ;  (* unit micron  *)
beamin={1.5,2.5}    (* unit micron  *) ;
Plot3D[Re[BeamOut[beamin, lens,0,0,y,z]],
		         {z,19980,20020},{y,-10,10}
, PlotPoints->{120,100}, BoxRatios->{2,1,0.5}
,Mesh->False, PlotRange->All, ViewPoint->{0,0,10}]

lens={1,2}*10^4  ;  (* unit micron  *)
beamin={1.5,2.5}    (* unit micron  *) ;
Plot3D[Re[BeamOut[beamin, lens,3000,0,y,z]],
		         {z,19980,20020},{y,-3010,-2990}
, PlotPoints->{120,100}, BoxRatios->{2,1,0.5}
,Mesh->False, PlotRange->All, ViewPoint->{0,0,10}]

lens={2.5,4.9}*10^4  ;  (* unit micron  *)
beamin={4.8,12}    (* unit micron  *) ;
Plot3D[Abs[BeamOut[beamin, lens,0,0,y,z]]^2,
		         {z,50000,50500},{y,-200,200}
, PlotPoints->{100,50}, BoxRatios->{5,4,4}
,Mesh->False, PlotRange->All, ViewPoint->{0,0,10}]

When you go into the lab and try to focus a beam to as a small spot as possible, what do you do? Project: use whatever lenses you can find in the lab, see who can obtain the smallest focused spot.

Dr. Bujin Guo is using a pinhole and try to focus a beam through. Obtain from him necessary parameters to calculate how much the throughput is (% of power that gets through the pin hole).

Look at the diode beam output. Why do you think the beam diverges in one dimension so much more than the other? Write a theory for x-y asymmetric beam.

Diode beam focus: how do you couple a beam from a diode into a rectangular waveguide? (which orientation to use?)

Where is the beam focused? 
 d -f =((Subscript[z, 1])-f) ((f^2)/(((((Subscript[z, 1])-f))^2)+((ζ0)^2))) = ((Subscript[z, 1])-f)((f^2)/((ζ0)^2)) (1/((((((Subscript[z, 1])-f))^2)/((ζ0)^2))+1)) ≈ ((Subscript[z, 1])-f)((f^2)/((ζ0)^2))
               d=   f (1 + ((((Subscript[z, 1])-f)f)/((ζ0)^2))) ≈f
very near the focal point.

If we want to collimate and expand the beam to shoot far away, where do we put the beam?

use the relation:  d -f =((Subscript[z, 1])-f) ((f^2)/(((((Subscript[z, 1])-f))^2)+((ζ0)^2) )). If the input beam waist is right at the focal point: (Subscript[z, 1])=f, then:    d -f =0: the output beam waist is at the other focal point.

4.3.6 Numerical aperture of a lens

If we want a lens of focus or collimate a beam, how do we order it? Let's look at the external cavity:

link to Appl Optics

Does the lens diameter matter? Why do we want "aspheric" lens? what is the relation between distance (Subscript[L, 1]) and (Subscript[L, 2]) and focal length f ? Where should the location of the lens should be if we are given a fixed distance (Subscript[L, 1])+(Subscript[L, 2]).

What can the lateral slit do for us?

4.3.7 Summary

Now we know why Gaussian beam is so interesting: it remains a Gaussian (although transformed) through lenses or curve mirrors. All the above relations with lens also apply with spherical (or paraboloidal mirror).

- With appropriate conditions, we can find a set of Gaussian solutions that remain invariant in a system of optics: laser cavity or lens guide (obsolete applications).

- We'll see that Gaussian turns out to be just one solution among many cavity modes; not surprisingly, it is the simplest solution, or the lowest order solution.

- More importantly: many concepts above are NOT exclusive to Gaussian beam. Even non-Gaussian beams have similar behavior: this is because they all are governed by diffraction. What we learn is the diffraction of paraxial light beam. Gaussian is a special case.

In the vertical dimension, a diode beam coming out of its facet has a flat wavefront and an amplitude profile like this:

Plot the wave profile of the beam as it leaves the diode facet along the x and z dimension. Use reasonable assumption of everything. (you must figure it out yourself).

Assignment: do all the excersizes in section 4.3 above and this:

John measured a beam with a Gaussian intensity profile (e^(-((2(ρ^2))/(w^2)))) and obtained w= 2.5 μm. The wavelength was λ=1 μm. He then measured the beam divergence which has a full width half max (FWHM) of 35 degree. Is this posible?

5. Advanced demonstration: Large NA beam

More accurate calculation needed for very large NA beam. Polarization needs be included. Illustration: