Light Propagation (part 1)
Spatial Behavior

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

Segment 4

8. (Advanced topics 6358/7358) Introduction to the diffraction concept and far-field Fraunhofer approximation

8.1 Linear superposition and diffraction formulism

The general aspect of Gaussian beam properties we discussed above is NOT specific to just Gaussian beam. In fact, it comes from diffraction behavior, a topic we will study in depth later. We use Gaussian beam as a study object because it is mathematically convenient, and is also important in terms of practicality. But we can begin to gain some insight by consider the problem more generally.

In the above, we consider this expression: 
                         E[x,z;t]=(e^(-i ω t))(∫_(-∞)^∞)c[(Subscript[k, x])] (e^(i (x (Subscript[k, x]) + z (Sqrt[((Subscript[k, 0])^2)ϵ- ((Subscript[k, x])^2)])))) d(Subscript[k, x])      (7.2.1b)
which is about constructing a solution with linear superposition of plane wave (e^(i (x (Subscript[k, x]) + z (Sqrt[((Subscript[k, 0])^2)ϵ- ((Subscript[k, x])^2)])))) weighed by the function c[(Subscript[k, x])], which is called the spatial spectral function of the beam. We see that if we know the field at a particular plane, e. g. z=0, then 
                    E[x,z=0;t=0]=(∫_(-∞)^∞)c[k] (e^(i x k)) dk                     (8.1.1)
where we drop the subscript x for convenience.

Then:
                         c[k] =(1/(2 π))(∫_(-∞)^∞)(Subscript[E, 0])[x] (e^(-i x k)) dx                          (8.1.2)
as the inverse Fourier transform of (Subscript[E, 0])[x]≡E[x,z=0;t=0] .
We can substitute (8.1.2) into (7.2.1b) above and eliminate c[k]:
                 E[x,z;t]=(e^(-i ω t))(1/(2 π))(∫_(-∞)^∞)(∫_(-∞)^∞)(Subscript[E, 0])[ξ] (e^(-i ξ k)) dξ (e^(i (x k + z (Sqrt[((Subscript[k, 0])^2)ϵ- (k^2)])))) dk      (8.1.3)
Note that we cannot use the same variable x because it is a dummy integration variable and hence we replace it with ξ to distinguish it from the true variable x.

Now we can switch the integration order:
                         E[x,z;t]=(e^(-i ω t))(∫_(-∞)^∞) dξ (Subscript[E, 0])[ξ]  (1/(2 π))(∫_(-∞)^∞)  (e^(-i ξ k)) (e^(i (x k + z (Sqrt[((Subscript[k, 0])^2)ϵ- (k^2)])))) dk      (8.1.4)
The integral                               (1/(2 π))(∫_(-∞)^∞)  (e^(-i ξ k)) (e^(i (x k + z (Sqrt[((Subscript[k, 0])^2)ϵ- (k^2)])))) dk                     (8.1.5)
is a function that links the field in one plane (Subscript[E, 0])[ξ]  to the rest of field further away: this is the essence of the diffraction integral concept. But this integral is very difficult to calculate hence, we need approximation like we did above.

8.2 Paraxial approximation

For paraxial approximation, we let:
                   (Sqrt[((Subscript[k, 0])^2)ϵ-(k^2)]) =(Sqrt[((Subscript[k, n])^2)-(k^2)]) =(Subscript[k, n])(Sqrt[1-((k^2)/((Subscript[k, n])^2))]) ≈(Subscript[k, n]) (1-((k^2)/(2 ((Subscript[k, n])^2)))) = (Subscript[k, n]) - ((k^2)/(2 (Subscript[k, n])))    (8.2.1)
Let’s take a detailed look at the paraxial approximation concept.

8.2.1 Paraxial phase front in r-space,

The above approximation is made in k-space. But for illustration, we will look at the paraxial approximation in r-space, which is {x,y,z}. Of course the 2 spaces are different, but the math concept is equivalent. For  r-space, let’s look at the phase front. This is the phase front of a spherical wave, which (obviously) is spherical!

y=0 ;
ContourPlot[Arg[(E^(I 2 π (Sqrt[x^2+y^2+z^2]))) ]==0,{z,-10,10}, {x,-10,10}]

Now if we have very large z (or any other main propagation direction), and we look in the z x plane, we can compare 2 wavefronts: spherical and paraxial for a small viewing angle in the propagation direction:

y=0. ;
ContourPlot[{Arg[(E^(I 2 π (Sqrt[x^2+y^2+z^2]))) ]==0, Arg[(E^(I 2 π  (z+((x^2)/(2 z))))) ]==0}
,{z,40,60}, {x,-10,10},ContourStyle-> {{Thick, Red},{Dashed, Thick, Blue}}]

We see that the paraxial approximation works very well along path nearly parallel to the axis: paraxial. Of course, we don't need contourplot to even plot it out. We know:
                     (z+((x^2)/(2 z)))== m

z= .; x=. ;
Solve[(z+((x^2)/(2 z)))== m ,{x}]

m=50; 
Plot[{Table[{- (Sqrt[(m^2) -(z^2)]),  (Sqrt[(m^2) -(z^2)])} ,{m,40,60,2}],Table[{ -(Sqrt[2]) (Sqrt[m z-(z^2)]),(Sqrt[2]) (Sqrt[m z-(z^2)])} ,{m,40,60,2}]}, {z,30,60}, PlotRange-> {{0,60},{-40,40}},PlotStyle-> {Red,Blue},AspectRatio-> 8./6]

We will see that even for angle where the concept of nearly “parallel” is a stretch, e. g. 30 degree, paraxial approx is still reasonable. At larger angle where the deviation is obvious, it is actually still fine if appropriate correction terms are added. Paraxial is not applicable only for large angle >45

Continued

Here in k-space, the same paraxial approximation can be applied on the condition that the beam does not have much (Subscript[k, x])-component compared with (Subscript[k, z]). This means the beam has small divergence. However, we will still use the approximation with the understanding that the main beam in the narrow center is valid and usually what we care. (the beam portion with large divergence is usually very low in power).
With  the paraxial approximation:
          (1/(2 π))(∫_(-∞)^∞)  (e^(-i ξ k)) (e^(i (x k + z (Sqrt[((Subscript[k, 0])^2)ϵ- (k^2)])))) dk ~(1/(2 π)) (e^(i (Subscript[k, n]) z))(∫_(-∞)^∞) (e^(-i k ξ)) (e^(i (k x  -z ((k^2)/(2 (Subscript[k, n]))) )))dk            (8.2.2)
which is can be integrated in closed form:

k=. ;
(∫_(-∞)^∞) (E^(-I k ξ)) (E^(I (k x  -z ((k^2)/(2 q)) )))\[DifferentialD]k

Hence:      (1/(2 π))(∫_(-∞)^∞) (e^(-i k ξ)) (e^(i (k x  -z ((k^2)/(2 (Subscript[k, n]))) )))dk = (Sqrt[((Subscript[k, n])/(2 π i z))])(e^((i (Subscript[k, n]) (((x-ξ))^2))/(2 z)))            (8.2.3)
Substitute  (8.2.3) in (8.1.3):
                    E[x,z;t]=(e^(i ((Subscript[k, n]) z -ω t)))(Sqrt[((Subscript[k, n])/(2 π i z))])(∫_(-∞)^∞)dξ (Subscript[E, 0])[ξ] (e^((i (Subscript[k, n]) (((x-ξ))^2))/(2 z)))      (8.2.4)

We will see that for a finite aperture, this important expression  (8.2.4) is equivalent with the Fresnel diffraction integral.
We can expand the exponent term and obtain
                    E[x,z;t]=(e^(i ((Subscript[k, n]) z -ω t)))(Sqrt[((Subscript[k, n])/(2 π i z))])(∫_(-∞)^∞)dξ (Subscript[E, 0])[ξ] (e^((i (Subscript[k, n])  ((ξ^2)+(x^2)-2 ξ x))/(2 z)))                    
                   E[x,z;t]=(e^(i ((Subscript[k, n]) z -ω t)))(e^(i (Subscript[k, n])((x^2)/(2z)) ))(Sqrt[((Subscript[k, n])/(2 π i z))])(∫_(-∞)^∞)dξ (Subscript[E, 0])[ξ] (e^(-i (Subscript[k, n]) ξ (x/z) ))(e^(i (Subscript[k, n])  ((ξ^2)/(2z)) ))       (8.2.5)   
We can obtain an alternative approximation for (8.2.5) for very large z>>x:
                   E[x,z;t]=((e^(i ((Subscript[k, n]) r -ω t)))/(Sqrt[2 π  i (Subscript[k, n]) r]))(Subscript[k, n])(∫_(-∞)^∞)dξ (Subscript[E, 0])[ξ] (e^(-i (Subscript[k, n]) ξ (x/r) ))(e^(i (Subscript[k, n])  ((ξ^2)/(2r)) ))       (8.2.6)

8.3 Illustration

Consider we have a square beam:

Plot[UnitBox[x],{x,-1,1}, Filling-> {1-> {Axis, Blue}}]

Such a beam is not realistic at the discontinuous edge, but we will consider that as an approximation for a sharp cut off beam, such as a plane wave with an hard aperture.
Applying (8.2.4):                (∫_(-∞)^∞)dξ (Subscript[E, 0])[ξ] (e^((i (Subscript[k, n]) (((x-ξ))^2))/(2 z))) =(∫_(-w)^w)dξ  (e^((i (Subscript[k, n]) (((x-ξ))^2))/(2 z)))    (8.3.1)
where w is the beam half-width,

k=. ; (Sqrt[(k/(2 π I z))])(∫_(-w)^w) (E^((I k  (((x-ξ))^2))/(2 z)))\[DifferentialD]ξ

FullSimplify[((1/2)-(I/2))  (Sqrt[-((I )/2)])  ]

hence:  
                  E[x,z;t]= (e^(i ((Subscript[k, n]) z -ω t)))(i/2) (Erfi[(((e^(i π/4)) (Sqrt[k]) (x+w))/(Sqrt[2 z]))]-Erfi[(((e^(i π/4)) (Sqrt[k]) (x-w))/(Sqrt[2 z]))] )     (8.3.2)

Manipulate[k=2 π ;  
Plot3D[Re[(I/2)(E^(I k z )) (Erfi[((((1/2)+(I/2)) (Sqrt[k]) (x-w))/(Sqrt[z]))]-Erfi[((((1/2)+(I/2)) (Sqrt[k]) (x+w))/(Sqrt[z]))]) ] ,{x,-5,5},{z,0,5}, PlotRange-> All,BoxRatios->{2,1,1}, Mesh-> False,PlotPoints->{80,40}],{w,0.5,5}]

Manipulate[k=2 π ;  
DensityPlot[Re[(I/2)(E^(I k z )) (Erfi[((((1/2)+(I/2)) (Sqrt[k]) (x+w))/(Sqrt[z]))]-Erfi[((((1/2)+(I/2)) (Sqrt[k]) (x-w))/(Sqrt[z]))]) ] ,{z,0,5}, {x,-5,5},PlotRange-> All, Mesh-> False,PlotPoints->{80,40},ColorFunction->GrayLevel, AspectRatio-> 2],{w,0.5,5}]

compare with diffraction of water around an edge

Note that the solution at large angle is NOT accurate because it violates the assumption of paraxial approximation. The part with small divergence angle is more accurate. We will see that the farfield divergence vs. beam waist behaves the same as that of Gaussian beam.

8.4 Fraunhofer diffraction and farfield profile.

In many cases, we are not interested the field near the beam waist or aperture as shown above, but care more about the beam divergence behavior at a large distance away, which referred to as “far field” (FF). In this case, we don’t care to describe the farfield in terms of cartesian coordinates x and y, but rather in terms of viewing angle θ, i. e. polar coordinate. It turns out that we can also make further approximation as well to make the integration a lot easier, yet the FF approximation is still very good, and in fact, much better than at near field.

Let’s look at the quadratic term in the exponent (e^(i (Subscript[k, n])  ((ξ^2)/(2z)) )) of (8.2.5):
                                                      (Subscript[k, n])  ((ξ^2)/(2z))= ((π (ξ^2))/((Subscript[λ, n])z))= ((π  ξ)/(Subscript[λ, n])) (ξ/z) =π ((ξ/(Subscript[λ, n])) )/((z/ξ) )    (8.4.1)

In (8.4.1), ξ is the integration variable over which (Subscript[E, 0])[ξ] is appreciable. This approach is applicable only for (Subscript[E, 0])[ξ] with finite spatial extent, i. e. we can set a rough value (Subscript[ξ, M]) for which (Subscript[E, 0])[ξ]~ 0 for |ξ| > (Subscript[ξ, M]). For Gaussian for example, we can set (Subscript[ξ, M]) ~ 2 w. This is NOT an exact definition. Then, ({{((|ξ|)/(Subscript[λ, n])) <, ((Subscript[ξ, M])/(Subscript[λ, n]))}}) is roughly the size of the beam (or aperture) relative to the medium wavelength (Subscript[λ, n]). Likewise, (z/ξ) is the ratio of the location where we observe the beam relative to the beam size. 
Now, we see that if (z/(Subscript[ξ, M])) >> ((Subscript[ξ, M])/(Subscript[λ, n])), then (Subscript[k, n])  ((ξ^2)/(2z))<<1 over the significant region of integration in (∫_(-∞)^∞)dξ (Subscript[E, 0])[ξ] (e^(-i (Subscript[k, n]) ξ (x/z) ))(e^(i (Subscript[k, n])  ((ξ^2)/(2z)) )). 
What it means is that we can drop (e^(i (Subscript[k, n])  ((ξ^2)/(2z)) )) and our life gets easier (integration-wise) from (8.2.6):
                       E[x,z;t]=((e^(i ((Subscript[k, n]) r -ω t)))/(Sqrt[2 π i (Subscript[k, n]) r]))(Subscript[k, n])(∫_(-∞)^∞)dξ (Subscript[E, 0])[ξ] (e^(-i (Subscript[k, n]) ξ (x/r) ))(e^(i (Subscript[k, n])  ((ξ^2)/(2r)) )) 
                     ~((e^(i ((Subscript[k, n]) r -ω t)))/(Sqrt[2 π i (Subscript[k, n]) r]))(Subscript[k, n])(∫_(-∞)^∞)dξ (Subscript[E, 0])[ξ] (e^(-i (Subscript[k, n]) ξ (x/r) ))      (8.4.2)  
Now, we notice that:                             (x/r)=Sin[θ]                                         (8.4.3)
and we obtain the Fraunhofer diffraction expression:
                       E[θ,z;t]=((e^(i ((Subscript[k, n]) r -ω t)))/(Sqrt[2 π i (Subscript[k, n]) r]))(Subscript[k, n])(∫_(-∞)^∞)dξ (Subscript[E, 0])[ξ] (e^(-i (Subscript[k, n]) ξ Sin[θ] ))               (8.4.4)
what is this integral?                  (∫_(-∞)^∞)dξ (Subscript[E, 0])[ξ] (e^(-i ξ ((Subscript[k, n])  Sin[θ]) ))               (8.4.5)
It is the Fourier transform of (Subscript[E, 0])[ξ], the field at the aperture or in the plane of the field original measurement.

The far field (FF) beam profile (in the Fraunhofer FF condition) is approximately the Fourier transform of the nearfield beam profile. The Fraunhofer FF condition requires that the observation distance from the source plane is much larger than the source Rayleigh’s length:            z>>((π ((Subscript[ξ, M])^2))/(Subscript[λ, n])).

Note: if z<((π ((Subscript[ξ, M])^2))/(Subscript[λ, n])) or not much larger, it is said to be in the Fresnel zone to indicate that the quadratic term of Fresnel diffraction is neccessary, as opposed to be in the Fraunhofer (farfield) zone of the above.

8.5 Fraunhofer diffraction illustration

Let’s consider a uniform field with an aperture to give us a top-hat shape field profile as in 8.3 above:

Plot[UnitBox[x],{x,-1,1}, Filling-> {1-> {Axis, Blue}}]

We have seen that the field looks pretty complex at near field in 8.3. What about far field?
              (∫_(-∞)^∞)dξ (Subscript[E, 0])[ξ] (e^(-i ξ ((Subscript[k, n])  Sin[θ]) ))=(∫_(-w)^w)  (e^(-i ξ ((Subscript[k, n])  Sin[θ]) ))dξ       (8.5.1)

k=. ; x= .; z= .; ξ= . ; a= .;
(∫_(-w)^w)  (E^(-I k ξ Sin[θ]))\[DifferentialD]ξ

Hence:         E[θ,r;t]=((e^(i ((Subscript[k, n]) r -ω t)))/(Sqrt[2 π i (Subscript[k, n]) r]))2w (Subscript[k, n])((Sin[(Subscript[k, n]) w Sin[θ] ])/( (Subscript[k, n]) w Sin[θ]))               (8.5.2)
What does it look like?
The part ((e^(i ((Subscript[k, n]) r -ω t)))/(Sqrt[2 π i (Subscript[k, n]) r])) is a 2D circular wave as seen below:

Manipulate[
DensityPlot[Re[((E^(I 2 π ((Sqrt[(x^2)+(z^2)]) -t)))/((((x^2)+(z^2)))^(1/4)))],{z,100,110},{x,-20,20}, ColorFunction->GrayLevel, AspectRatio->4, PlotPoints->{40,80}, Mesh-> False],{t,0,1}]

The part ((Sin[(Subscript[k, n]) w Sin[θ] ])/((Subscript[k, n]) w Sin[θ])) should be plotted in polar coordinate:

Manipulate[
PolarPlot[Abs[((Sin[2 π w Sin[θ] ])/(2 π w Sin[θ])) ],{θ,-π/2,π/2} ,PlotStyle-> {Thick,Red}, PlotRange-> {{-0.1,1},{-0.3,0.3}}],{w,0.5,5}]

We see both part:

Manipulate[
DensityPlot[((Sin[2 π w (x/(Sqrt[(x^2)+(z^2)])) ])/(w 2 π Sin[(x/(Sqrt[(x^2)+(z^2)]))])) Re[((E^(I 2 π ((Sqrt[(x^2)+(z^2)]) -t)))/((((x^2)+(z^2)))^(1/4)))],{z,250,270},{x,-50,50}, ColorFunction->GrayLevel, AspectRatio->2.5, PlotPoints->{80,60}, Mesh-> False, PlotRange-> {-(1./(Sqrt[250])),(1./(Sqrt[250]))}],{{w,1.5},0.5,4},{t,0,1} ]

Manipulate[
ParametricPlot3D[{r  Cos[θ],r  Sin[θ] ,((Sin[2 π w Sin[θ] ])/(w 2 π Sin[θ])) Re[((E^(I 2 π  (r -t)))/(Sqrt[r]))]},{r,250.,260},{θ,-0.2*π/2,0.2*π/2},  BoxRatios->{20,150,20}, PlotPoints->{100,80}, Mesh-> False, PlotRange-> {{240,260},{-75,75},{-0.1,0.1}}],{{w,2.},1,4},{t,0,1} ]

Manipulate[
ParametricPlot3D[{r  Cos[θ],r  Sin[θ] ,(((((Sin[2 π w Sin[θ] ])/((Sqrt[r]) w 2 π Sin[θ]))))^2)},{r,250.,260},{θ,-0.5*π/2,0.5*π/2},  BoxRatios->{20,150,20}, PlotPoints->{20,70},Mesh-> False, PlotRange-> {{240,260},{-75,75},{0,0.005}}],{{w,2.},1,4} ]

Manipulate[z=250. ;
Off[DensityPlot::exclul] ;
Show[ DensityPlot[If[x== 0.,1. ,(((((Sin[2 π w (x/(Sqrt[(x^2)+(z^2)])) ])/(w  2 π Sin[(x/(Sqrt[(x^2)+(z^2)]))])) ))^2) ] (1/(Sqrt[(x^2)+(z^2)])),{x,-75,75},{y,-20,20} , ColorFunction->GrayLevel, AspectRatio->0.4, PlotPoints->{40,2}, Mesh-> False, PlotRange-> All]
,Graphics[{Hue[0.15,1,1,0.3],Rectangle[{-w,-20},{w,20}] }] ]
,{{w,1.5},0.5,4} ]

HW: Divergence from a diode laser facet

8.6 Advanced assignment: Simulate Young’s double-slit experiment.