Gaussian to Flat-top

Zemax Knowledge Base gives a design method to convert Gaussian irradiance profile to flat top profile (see reference [1]).
This post (A) gives another analytic form and (B) examines the 1-D situation.

A.
If the input Gaussian profile has the 1/e² radius of W, the output flat-top has the half-width of K, then the Zemax KB article gives:
        S = K [1-exp(-2X²/W²)]^1/2                                              (1)
This is a ray mapping formula to get the output coordinate value S for every input coordinate X.

For Zemax, this formula can be written in another form. Let R_EP be the radius of the entrance pupil, then

X W

=

X/REP W/REP

= XNA1/2

(2)

where X_N is the normalized coordinate and A is the system apodization factor. Eq(1) can therefore be written as
        S = K [1-exp(-2X_N²A)]^1/2                                              (3)
This formula is better suited for Zemax because both X_N and A are the Zemax direct parameters. Using this formula, the macro file for generating merit function can be simplified.

B. How about 1-D?
Very often people need a 1-D flat-top, e.g., flat-top on X-axis and Gaussian on Y-axis. In this case, the irradiance integration is different.

For input Gaussian profile, the contained power in a 1-D variable slit is:

Pi(X) = Ii ∫

X -X

exp(-2x2/Wx2)dx

∫

+∞ -∞

exp(-2y2/Wy2)dy

(4)

where I_i is the peak irradiance, W_x and W_y are the 1/e² beam radii. For output profile, the contained power in the 1-D variable slit is:

Po(S) = Io ∫

S -S

dξ

∫

+∞ -∞

exp(-2η2/Wη2)dη

(when S<=Wξ)

(5)

Po(S) = Io2Wξ

∫

+∞ -∞

exp(-2η2/Wη2)dη

(when S>Wξ)

(5)

where I_o is the peak irradiance, W_ξ is the half-width of the flat-top and W_η is the 1/e² half-width of the Gaussian profile.

First, let P_i(X=∞) = P_o(S=∞), the conservation of total power gives:

Io Ii

=

π1/2WxWy 21/22WξWη

(6)

To get the ray-mapping relation, let P_i(X) = P_o(S) and plug in Eq.(6):
        S = W_ξ erf (2^1/2X/W_x) = W_ξ erf (2^1/2X_NA^1/2).                                              (7)
where erf is error function. So for 2D and 1D flat-top generation, the ray-mapping formula is different.

To write a macro for generating merit function for 1-D flat-top, one can use
       GETSYSTEMDATA 1
       Apodization = VEC1(2)
to read the apodization factor.

To caculate erf function, Zemax does not have it built-in. One can use an approximation from wikipedia [2].

References:
[1] http://kb-en.radiantzemax.com/Knowledgebase/How-to-Design-a-Gaussian-to-Top-Hat-Beam-Shaper
[2] http://en.wikipedia.org/wiki/Error_function

Gaussian Beam Size

A.
For circular fundamental mode Gaussian beam, irradiance (intensity) profile is
     I(r) = exp(-2r²)                                              (1)
Here r is scaled by w (1/e² radius), i.e., when r=1, the irradiance drops to 1/e²=13.5% of the peak value.

1. D86 diameter
Assuming we measure power with a variable aperture whose diameter is ρ, the contained power is an integration of the irradiance:

P(ρ) = ∫

ρ 0

exp(-2r2)rdr

∫

2π 0

dθ =

π 2

[1-exp(-2ρ2)]

(2)

The total power is when ρ=∞: P_tot = P(ρ=∞) = π/2. Then
    P(ρ)/P_tot = 1-exp(-2ρ²)                                                                              (3)
The D86 diameter is when 86.5% of the total power is contained. From the above equation, D86 diameter is equal to the 1/e² diameter because P(ρ=1)/P_tot = 1-1/e² = 0.865.

2. D63 diameter
I see nowhere uses D63 diameter except one: the IEC Laser Safety Standard. The D63 diameter is when 63.2% of the total power is contained in a variable aperture. From Eq.(3), P(ρ=1/$\sqrt{2}$)/P_tot = 0.632. So D63 diameter is 1/$\sqrt{2}$=0.707 times the D86 diameter. From Eq.(1), I(1/$\sqrt{2}$)=1/e. So D63 diameter is the width at 1/e intensity points.

3. Knife edge width
Now we measure power with a traveling knife edge along the x-axis, the transmitted power is an integration of the irradiance:

P(a) = ∫

a -∞

exp(-2x2)dx

∫

+∞ -∞

exp(-2y2)dy =

π 4

[1+erf($\sqrt{2}$ a)]

(4)

where erf is error function. The total power P_tot=P(∞)=π/2. So
    P(a)/P_tot = 1/2 [1+erf($\sqrt{2}$ a)]                                                         (5)
From Eq.(5), P(-0.5)/P_tot = 0.159 and P(0.5)/P_tot = 0.841. This means that the knife edge traveling distance between the 15.9% and 84.1% power levels is equal to the 1/e² radius.
Alternatively, if one measures 10% - 90% power levels, the distance is equal to 1.281*1/e² radius.

4. Full width or half width at half maximum
Let I(r) = exp(-2r²) = 0.5, we get ln(0.5) = -2r², or r = $\sqrt{\mathrm{(ln2)/2}}$ = 0.5887. So when 1/e² radius is 1, the half width at half maximum is 0.5887.

5. D4 sigma (2nd moment) diameter
This is the ISO standard (ISO 11146-1). For the fundamental Gaussian, this definition and the traditional 1/e² definition are identical. This definition heavily weights the tails or outer wings of the intensity profile; so for non-ideal beams having side-lobes, the 2nd moment diameter can be substantially larger than their central lobe diameter. I will not go through the math detail here.

B.
For asymmetric 1-D Gaussian beam (Gaussian form at x, arbitrary form at y), irradiance (intensity) profile is
     I(x,y) = exp(-2x²) I(y)                                              (6)
Here x is scaled by w_x (1/e² half-width). I use "width" to replace "diameter" (and "half-width" to replace "radius") for non-circular beam.

1. D86 width
Instead of a circular aperture, a variable 1-D rectangular aperture is used here. Assuming a is the half-width of the aperture, the transmitted power is

P(a) = ∫

a -a

exp(-2x2)dx

∫

+∞ -∞

I(y)dy =

$\sqrt{\mathrm{\pi /2}}$

erf($\sqrt{2}$ a)

∫

+∞ -∞

I(y)dy

(7)

Divided by total power P_tot=P(∞),
   P(a)/P_tot = erf($\sqrt{2}$ a)                                                        (8)
Solve for erf($\sqrt{2}$ a)=0.865, we get a=0.747. So for 1-D Gaussian beam, D86 width is 0.747 times the 1/e² width.

2. D63 width
Using Eq.(8), solve for erf($\sqrt{2}$ a)=0.632, we get a=0.450. So for 1-D Gaussian beam, D63 width is 0.45 times the 1/e² width. Also for 1-D Gaussian beam, D63 width is 0.45/0.707=0.636 times smaller than that of the 2-D circular Gaussian beam.

3. Knife edge width
From Eq.(4), since the knife edge scan is a 1-D scan, the knife edge width is the same for 1-D and 2-D Gaussian.

4. FWHM
FWHM is an "intensity point relative to peak" type of width, it has nothing to do with contained power, so it is the same for 1-D and 2-D beams, i.e., it is still 0.5887 times the 1/e² width.

5. Second-moment width
2nd-moment width is inherently defined for separate x and y, so it is the same for 1-D and 2-D beams.

Summary:
Use this table to convert between different Gaussian beam sizes:

	Circular 2D Gaussian	1D Gaussian
1/e2 half-width	1	1
D86 half-width	1	0.747
1/e half-width	0.707	0.707
D63 half-width	0.707	0.45
Knife edge width 15.9%-84.1% clip	1	1
Knife edge width 10%-90% clip	1.281	1.281
HWHM	0.5887	0.5887

Note: for ideal fundamental Gaussian form only.

References:
[1] User Manual, ModeMaster PC - M² beam propagation analyzer. Coherent inc.
[2] A.E.Siegman, "How to Measure Laser Beam Quality". PDF copy available online.

Laser Beam Collimation

This post discusses one question: how to achieve the best collimation of a laser beam?

To collimate a laser diode or fiber's output using a lens, a common belief is to make the output beam's waist location as distant as possible form the lens.  For example, see Saleh and Teich [1].  Well, this is only approximately correct.

The best collimation is achieved when the output beam has either:
(1) the smallest divergence angle; or
(2) the largest waist size; or
(3) the longest Rayleigh range.
The above three conditions are equivalent, i.e., achieving any single one means the other two are achieved simultaneously.

To reach the above conditions, the input waist needs be at exactly the object focal plane of the collimating lens. Then the output beam waist will be at exactly the conjugate (image) focal plane [2]:

Figure 1: Geometry of the perfect collimation.

I plot the output waist location (relative to the image focal plane) and the output waist diameter as a function of the input waist location (relative to the object focal plane). This plot shows that (a) output waist size is at maximum and (b) waist location is at the image focal plane when input beam waist is at the object focal plane. (Assuming input beam waist diameter = 5 um, lens focal length = 10 mm, wavelength = 1 um.)

Figure 2: Output beam properties vs. input beam waist location. Notice my custom sign rule: if "relative to" is at left of the reference point, the sign is plus. For example, positive waist location means it is at the left side of the focal plane.

This figure also indicates that the output beam's waist location cannot be placed at infinity.  The farthest distance one can make is when the input beam's waist is located at Z_R + f from the lens (Z_R is input beam's Rayleigh range and f is focal length). This is the common belief that the best collimation should be reached; however, the output waist size is not the largest.



References:
[1] Saleh and Teich, Fundamentals of Photonics, 2nd Ed. Page 90.
[2] Sidney A. Self, "Focusing of spherical Gaussian beams," Applied Optics 22, 658 (1983).

Diffraction theory summary

This is a summary of the scalar diffraction theory after reading Ref. 1 and 2.

0. Scalar diffraction theory neglects the vectorial nature of the e&m wave. Two conditions need be met for scalar theory to produce accurate result: (1) the diffracting aperture larger than wavelength; and (2) the observing field not too close to the aperture. (Ref.1, p.35)

1. Scalar diffraction formula with least approximation:

In figure above, Σ is the diffration aperture, P₁ is a point in it and U₀(P₁) is the complex amplitude waiting to be diffracted, P is the point we want to solve. The complex amplitude at P is:

U(P) =

1 jλ

∫∫ U0(P1)

ejkr r

cos(n, r) ds

(1)

The approximation needed here is r >> λ . This is the Rayleigh-Summerfeld solution.

In frequency domain,  let G_z(f_x, f_y) and G_z=0(f_x, f_y) be the Fourier transform of U(P) and U₀(P₁) respectively; they are the amplitudes of each fundamental plane-waves (decomposed by Fourier transform).  By solving Helmholtz wave equation, Eq.(1) in frequency domain becomes (Ref.2, p.83):
     G_z(f_x, f_y) = G₀(f_x, f_y) exp[jkz(1-α²+β²)^1⁄₂]                     (2)
where
    f_x = α/λ,  f_y = β/λ                                                                  (2a)
are the spatial frequencies on 2-D planes normal to the z-axis; they depend on these plane-wave's propagation direction and α, β are the directional cosines relative to x and y axes respectively.  Eq.(1) analyzed in frequency domain is called "Angular Spectrum" approach (Ref.1, p.55).

2. Fresnel diffraction:
If further approximation is made:
    cos(n,r) ≈ 1                                                                                (3)
then

U(x, y) =

ejkz jλz

∫∫ U0(x1, y1) exp[jk

(x-x1)2+(y-y1)2 2z

] dx1dy1

(4)

2.1
Eq.(4) can be written as:

U(x, y) =

ejkz jλz

∫∫ U0(x1, y1) h(x-x1, y-y1) dx1dy1

(5)

=

ejkz jλz

U0(x, y) ∗ h(x, y)

(5a)

where

h(x, y) = exp(jk

x2+y2 2z

)

(5b)

Eq.(5)-(5b) have the following physical meaning:
• Eq.(5b): the impulse response h is the complex amplitude of the spherical wave from a unit point source at (x, y) = (0, 0).
• Eq.(5a): the diffracted field is simply a convolution between the initial field U₀ and the impulse response h.
• Fresnel diffraction is a space-invariant linear system.

2.2
Eq.(4) can also be written as:

U(x, y) =

ejkz jλz

exp(jk

x2+y2 2z

) F[U0(x1, y1) exp(jk

x12+y12 2z

)]

(6)

Eq.(6) means:
• The diffracted field is a Fourier transform of U₀ multiplied by a diverging spherical wave.

3. Fraunhofer diffraction:
If an even more stringent approximation is made:

z >>

(x12+y12)max 2λ

(7)

then

U(x, y) =

ejkz jλz

exp(jk

x2+y2 2z

) F[U0(x1, y1)]

(8)

So Fraunhofer diffraction is simple: the diffracted field is simply a Fourier transform of the initial field.

Eq.(7) is more often expressed by "Fresnel number":

F =

R2 λz

(7a)

where R is aperture's radius. Fraunhofer diffraction is when F << 1. 

In summary,
     • When diffraction distance >> λ , scalar theory can be used. In frequency domain, the angular spectrum approach is very physically intuitive.
     • When diffraction angle is small (Eq.(3)), Fresnel diffraction is used.
     • When Eq. (3) and (7) are both satisfied, Fraunhofer diffraction can be used.

References:
[1] Goodman, "Introduction to Fourier Optics", 2nd Ed.
[2] 黄婉云，《傅里叶光学教程》，1985年第1版.

Near Field and Far Field

1. In diffraction optics,

Fresnel number:

F =

R2 λz

(1)

where R is aperture's radius.

Far field:
     • F << 1.
     • It is the Fraunhofer diffraction area.
     • The diffraction pattern remains its shape; only its size gets bigger at farther distance.  

Near field:
     • F ≥ 1.
     • It is the Fresnel diffraction area.
     • The diffraction pattern changes shape as a function of distance.

 There is also an "extreme near field" region where the field observed behind the aperture is simply a geometrical projection of the aperture. Geometrical optics zone is when F >> 1.

2. For laser beam,

Rayleigh range:

z0 =

πW02 λ

=

π w02 λ M2

(2)

where W₀ is 1/e² waist radius and w₀ = W₀/M is the 1/e² waist radius of the embedded fundamental Gaussian beam. (i.e., a Gaussian beam and its embedded Gaussian beam have the same Rayleigh range).

Near field is near the laser beam waist, within its Rayleigh range.
Far field is outside of the Rayleigh range.

References:
[1] http://en.wikipedia.org/wiki/Fresnel_number
[2] 黄婉云，《傅里叶光学教程》，1985年第1版. Page 95.

Zemax POP macro: normalized plot and calculation of beam size

Zemax Physical Optics Propagation (POP) generates cross-sectional beam profile.  Two functions are lacking that are useful to me:
1. An irradiance-normalized plot. (The absolute "watts-per-sq-millimeter" value does not matter to me almost all the time).
2. Zemax gives "second-moment" beam diameter but again, almost all the time I just need a simple 1/e^2 diameter regardless of perfect Gaussian beam or not. (The second-moment beam diameter definition heavily weighs the tails or outer wings of the intensity profile. Therefore beams with side-lobes will have second-moment width substantially larger than their central lobe widths. [1])

So I wrote a macro file to do these two tasks:

!-----------------------------------------------------------
!FIRST, SAVE SETTINGS IN POP WINDOW. 
!The "Show As" in "Display" tab must be set as "Cross X (or Y)".
!----------------------------------------------------------- 
!Input POP parameters saved in POP window:
!-----------------------------------------------------------
sampling = 512
zoom = 8

AA$ = $TEMPFILENAME()
GETTEXTFILE AA$, POP    # Import POP text data file
OPEN AA$

!-----------------------------------------------------------
! Create a loop to read POP data
!-----------------------------------------------------------
invalid_row = 13        # The first 13 rows are text, not data
For i,1,invalid_row,1
    READ a,b
Next i

Grid_size = sampling/zoom        # POP grid size (=Sampling/zoom)
DECLARE xx, DOUBLE, 1, Grid_size+1                   
DECLARE Irrad, DOUBLE, 1, Grid_size+1
For i,1,Grid_size+1,1
    READ a,b
    xx(i) = a
    Irrad(i) = b
    print xx(i), ", ", Irrad(i)
Next i
CLOSE

!-----------------------------------------------------------
! Find peak
!-----------------------------------------------------------
Max_i = 1
For j,2,Grid_size+1,1
!    PRINT "j=",j,"xx=",xx(j),"Irrad=",Irrad(j)
    IF ( Irrad(j)>Irrad(j-1) ) & ( Irrad(j)>Irrad(Max_i) ) THEN Max_i = j
Next j

!-----------------------------------------------------------
! Normalize lineshape
!-----------------------------------------------------------
DECLARE Irrad_norm, DOUBLE, 1, Grid_size+1
For j,1,Grid_size+1,1
    Irrad_norm(j) = Irrad(j)/Irrad(Max_i)
    PRINT "j=",j," xx=",xx(j)," Irrad_norm=",Irrad_norm(j)
Next j
!-----------------------------------------------------------
! Find 1/e^2 radius at extreme left (linear interpolation)
!-----------------------------------------------------------
Radius_Li = 1
For j,1,Max_i,1
    IF ( Irrad_norm(j)<=0.135 ) & ( Irrad_norm(j+1)>=0.135 ) THEN GOTO 1    #once the condition is met, jump out loop.
Next j
LABEL 1
Radius_Li = j
R_L1 = xx(Radius_Li)
R_L2 = xx(Radius_Li+1)
Irrad_L1 = Irrad_norm(Radius_Li)
Irrad_L2 = Irrad_norm(Radius_Li+1)
R_L = (0.135-Irrad_L1)*(R_L1-R_L2)/(Irrad_L1-Irrad_L2) + R_L1
PRINT "At ", R_L, " Irradiance = 0.135"
!-----------------------------------------------------------
! Find 1/e^2 radius at extreme right (linear interpolation)
!-----------------------------------------------------------
Radius_Ri = 1
For j,Max_i,Grid_size,1
    IF ( Irrad_norm(j)>=0.135 ) & ( Irrad_norm(j+1)<=0.135 ) THEN Radius_Ri = j+1
Next j
R_R1 = xx(Radius_Ri)
R_R2 = xx(Radius_Ri-1)
Irrad_R1 = Irrad_norm(Radius_Ri)
Irrad_R2 = Irrad_norm(Radius_Ri-1)
R_R = (0.135-Irrad_R1)*(R_R1-R_R2)/(Irrad_R1-Irrad_R2) + R_R1
PRINT "At ", R_R, " Irradiance = 0.135"
!-----------------------------------------------------------
! Calculate 1/e^2 beam diameter
!-----------------------------------------------------------
dia = R_R - R_L
PRINT "Beam Diameter = ", dia
!-----------------------------------------------------------
! Plot normalized lineshape
!-----------------------------------------------------------
PLOT NEW
PLOT DATA, xx,Irrad_norm,Grid_size+1,0,0,0
PLOT BANNER, "Normalized POP beam profile"
PLOT RANGEX, xx(1), xx(Grid_size+1)
PLOT RANGEY, 0, 1
PLOT COMM1, $DATE()
PLOT COMM2, $FILEPATH()
PLOT COMM3, "1/e^2 Gaussian Diameter = ", $STR(dia)
PLOT GO

Save the above text as a "xxx.zpl" file in Zemax's ZPL folder (the folder path can be viewed/changed in File -> Preferences). To run this macro, press F9, choose this file and click "Execute".

Here is an example. The POP window generated by Zemax shows an aberrated beam profile:


My macro generates a plot shown below.  First, it is normalized.  Second, the 1/e^2 beam diameter is given.



Compare three beam sizes:
1. POP pilot beam radius (paraxial 1/e^2) = 0.0070
2. POP "second moment" beam radius = 0.0071
3. Simple 1/e^2 radius calculated by my POP macro = 0.0024
The three sizes are different. Most often I just need know the 3rd value.

Next: It would be nice to also generate a Gaussian-fit curve.  I will perhaps write a MATLab-linked code to do this task.

Reference
[1] A.E.Siegman, "How to measure laser beam quality". PDF available on web by googling it.