Intensity distribution of Gaussian beam Propagated in space through a Lens at twice the focal Length

[SudipKCTUM]

Hi there!

I wanted to express my gratitude for developing the incredibly useful LightPipe module. Currently, I'm using it to simulate the propagation of a Gaussian beam through space via a lens over a distance of 10 km, with a focal length of 5000 using the Fresnel propagator command.

Specifically, I'm interested in simulating the beam's propagation through a normal lens (thin lens) at a distance twice the focal length. My expectation was that the intensity pattern of the Gaussian beam would remain unchanged at twice the focal length, similar to its initial state (Near field), provided I chose the right grid size to match the beam's width. However, I observed a spread in the intensity distribution of the Gaussian beam with larger propagation distances, which was unexpected.

I'd greatly appreciate your assistance in adjusting my code to ensure that the intensity pattern remains unchanged when the beam reaches twice the focal length through the lens. This adjustment is crucial as it will facilitate further analysis of the intensity distribution in subsequent repetitions, such as when the beam bounces off a mirror.

Any advice or guidance you can provide to resolve this matter would be immensely helpful. Thank you for considering my request.

Sincerely, Sudip

Code 1

from LightPipes import *
import matplotlib.pyplot as plt
import numpy as np
from scipy.interpolate import interp1d

def find_fwhm(x, y):
    half_max = max(y) / 2.0
    sign_change_points = np.where(np.diff(np.sign(y - half_max)))[0]
    if len(sign_change_points) >= 2:
        left_root = np.interp(half_max, [y[sign_change_points[0]], y[sign_change_points[0] + 1]], [x[sign_change_points[0]], x[sign_change_points[0] + 1]])
        right_root = np.interp(half_max, [y[sign_change_points[1]], y[sign_change_points[1] + 1]], [x[sign_change_points[1]], x[sign_change_points[1] + 1]])
        return right_root - left_root, (left_root, right_root)
    else:
        return None, None

# Initial settings
GridSize = 77* 1e-2  # 75 cm in meters
GridDimension = 2000
lambda_ = 1064 * 1e-9  # 1064 nm in meters

Field = Begin(GridSize, lambda_, GridDimension)

w = 5 * 1e-2  # Width of the Gaussian Beam in meters
Field = GaussAperture(w, 0, 0, 1, Field)

I0 = Intensity(1, Field)
phi = Phase(Field)

# Straight-forward lens
fn = 5000 #  m
zn = 10000 # m

F = Lens(fn, Field)  # 50 meters focal length
F = Fresnel(Field, zn)  # 100 meters propagation
I1 = Intensity(1, F)
phi1 = Phase(F)

# Method with spherical coordinates
f = 5000  # 50 meters focal length
f1 =10000  # 100 meters propagation
f2 = f1 * f / (f1 - f)
frac = f / f1
newsize = frac * GridSize

F0 = Lens(f, Field)
F1 = LensFresnel(f2, f, F0)
F1 = Convert(F1)
phi2 = Phase(F1)
I2 = Intensity(1, F1)

# Plotting
fig, axes = plt.subplots(2, 3, figsize=(15, 12))

x = [(-GridSize / 2 + i * GridSize / GridDimension) / 1e-3 for i in range(GridDimension)]
y = [(-newsize / 2 + i * newsize / GridDimension) / 1e-3 for i in range(GridDimension)]

# Intensity line plots
ax1, ax2, ax3 = axes[0]
ax1.plot(x, I0[int(GridDimension / 2)])
ax1.set_title('Near Field')
ax1.set_xlabel('x [mm]')
ax1.set_ylabel('Intensity [a.u.]')

ax2.plot(x, I1[int(GridDimension / 2)])
ax2.set_title('Normal lens')
ax2.set_xlabel('x [mm]')
ax2.set_ylabel('Intensity [a.u.]')

ax3.plot(y, I2[int(GridDimension / 2)])
ax3.set_title('Spherical Coordinate')
ax3.set_xlabel('x [mm]')
ax3.set_ylabel('Intensity [a.u.]')

# Add suptitle
# Add suptitle with dynamic values
plt.suptitle(f'Propagation of a Gaussian beam in free space through lens (f = {fn:.2f} m, z = {zn:.3f} m, Grid Size ={GridSize:.3f} m)')

# Calculate and plot FWHM
for ax, data, x_vals in zip([ax1, ax2, ax3], [I0, I1, I2], [x, x, y]):
    intensity_profile = data[int(GridDimension / 2)]
    fwhm, roots = find_fwhm(x_vals, intensity_profile)
    if fwhm:
        ax.axvline(roots[0], color='r', linestyle='--')
        ax.axvline(roots[-1], color='r', linestyle='--')
        ax.text(0.05, 0.9, f'FWHM: {fwhm:.2f} mm', transform=ax.transAxes, color='red')

# Intensity images
ax1, ax2, ax3 = axes[1]
ax1.imshow(I0, cmap='rainbow')
ax1.axis('off')
ax1.set_title('Near Field')

ax2.imshow(I1, cmap='rainbow')
ax2.axis('off')
ax2.set_title('Normal Lens')

ax3.imshow(I2, cmap='rainbow')
ax3.axis('off')
ax3.set_title('Spherical Coordinate')

plt.tight_layout()
plt.show()

Code 2

from LightPipes import *
import matplotlib.pyplot as plt

from mpl_toolkits.mplot3d import Axes3D
import numpy as np

# Initial settings
GridSize = 77* 1e-2  # 75 cm in meters
GridDimension = 1000
lambda_ = 1064 * 1e-9  # 1064 nm in meters

# Initialize the field
Field = Begin(GridSize, lambda_, GridDimension)

# Gaussian beam parameters
w = 5 * cm  # Width of the Gaussian Beam
Field = GaussAperture(w, 0, 0, 1, Field)

# Propagate through lens to focal point 5000 m
fn = 5000 *m
zn = 10000 *m
Field = Lens(fn * m, Field)

# #Method with spherical coordinates
# f = 5000*m
# f1 = 10000*m
# f2 = f1*f/(f1-f)
# frac = f/f1
# newsize = frac*GridSize

# F0 = Lens(5000*m,Field)
# F1 = LensFresnel(f2,f, F0)
# F1 = Convert(F1)
# phi2 = Phase(F1)
# I2 = Intensity(1,F1)

# Define propagation distance and steps
start_distance = 0 * m
end_distance = zn * m
num_steps = 30
step_size = (end_distance - start_distance) / num_steps

# Initialize array to store intensities
intensity_profiles = []

# Propagate the field in small steps and store the intensity
for i in range(num_steps):
    distance = start_distance + i * step_size
    Field = Fresnel(Field, step_size)
    intensity_profiles.append(Intensity(1, Field))

# Convert intensity_profiles to a numpy array for easier plotting
intensity_profiles = np.array(intensity_profiles)

# Create a meshgrid for plotting
x = np.linspace(-GridSize/2/mm, GridSize/2/mm, GridDimension)
z = np.linspace(start_distance/m, end_distance/m, num_steps)
X, Z = np.meshgrid(x, z)

# Plotting the intensity distribution
plt.figure(figsize=(10, 8))
plt.imshow(intensity_profiles[:, int(GridDimension/2), :], extent=[x.min(), x.max(), z.min(), z.max()],
           aspect='auto', cmap='rainbow')
plt.colorbar(label='Intensity [a.u.]')
plt.xlabel('x [mm]')
plt.ylabel('Propagation Distance [m]')
# plt.title('Intensity Distribution from 0 m to 10000 m, Gaussian aperture')
plt.suptitle(f'Propagation of a Gaussian beam in free space through lens (f = {fn:.2f} m, z = {zn:.3f} m, Grid Size ={GridSize:.3f} m)')
plt.show()
plt.show()

[ldoyle]

Dear Sudip,

thanks for the kind words. I have not tested your code in detail, but I instead investigated a different aspect of your problem: While the code may also have some problems, I believe it is correct that you do not get exactly the same beam diameter back.

This relates to a discussion I had with fellow physicists a couple of times in the past: due to diffraction, even a "collimated" Gaussian beam will diverge. The strength of this effect is indicated/calculated by the Rayleigh range $z_R$. Now here comes the interesting bit: If the waist of the Gaussian beam lies on the plane of the lens, the beam coming onto the lens "from the left", cannot be collimated. Instead, by placing the waist exactly on the lens, we imply the condition that an even larger Gaussian beam (coming from "the left") converges to its waist $w_0$ on the lens. This becomes important once the focal length is not much smaller than the Rayleigh range (i.e. $z_R >> f$ is not fulfilled) as we will see below.

I therefore usually approach the problem from a different perspective: from ray optics and the "thin-lens equation" we know that at $z=2f$ an image is formed of an object at distance $z=-2f$. This image has a magnification of $1$, so indeed the image size at $2f$ behind the lens must always match the beam diameter at $2f$ before the lens. Thus, I would recommend placing the waist not on the lens, but a distance $2f$ before the lens.

For your beam parameters ($\lambda=1064nm$, $w_0=50mm$) and (purposely chosen) focal length of e.g. $f=5m$, this yields the following plot (code and calculations below):

For your $w_0=50mm$ beam (58.8mm FWHM), the Rayleigh range is $7.3km$ and the beam therefore "quasi-collimated" on the scale of meters. On the plot, the incoming beam looks parallel, the focus is diffraction limited at $\approx 34\mu m$ and the beam expands to exactly 50mm after a distance $2f$.

Now let's try with your desired focal length $f=5000m$. This is on the same order as the Rayleigh range, so the picture changes drastically:

By placing the beam waist at $z=-2f$ it is evident that diffraction is not irrelevant left of the lens. Since the beam diverges, the smallest spot behind the lens is not at a distance $f$ and probably not of diffraction limited size for the given beam diameter, instead it comes at $\approx 6.6km$ (this is not analytic, but the next calculated step, so maybe slightly off). Nevertheless, at a distance $2f$ behind the lens, the original beam diameter is perfectly preserved.

Now let's try placing the waist very close to the lens ($1m$) which is the setup you are simulating in your code:

Surprisingly, the smallest waist is closer than $f$ at $\approx 3.4km$. The beam waist after 10km is 84mm (FWHM intensity $\approx 99mm$). (Note: I calculated the FWHM in intensity from the waist in E-field, so there is a $\sqrt{2}$ that sneaks in, I hope I got it right.)

With all that being said, I think you should consider what effect the Rayleigh range has for the setup parameters you chose and where a suitable distance to place the beam waist before the lens would be. Now some comments on the coding side (again, I have not run your code, just looking at it and testing my own code):

• In your first example, you might have a typo rendering the lens ineffective (F = Lens(fn, Field); F = Fresnel(Field, zn) should probably be ... F = Fresnel(F, zn))
• In your second method, I will have to double check the lens coefficients, maybe you are applying a too strong lens (F0 = Lens(f, Field) maybe should be F0 = Lens(f1, Field))
• Given my prior discussions with colleagues about this Rayleigh range problem, I instead decided to create a new code which uses the analytic formulas for tracing a Gaussian beam through a thin lens system (based on ABCD matrix formalism and "complex beam parameter")
• While coding this, I noticed there may be a typo in LightPipes leading to errors when I coded the example with LightPipes.ABCD
• I therefore re-implemented the necessary formulas directly in the code file (below). This also made it much faster, since I am not operating on any 2D grid for this example.
• While typing my answer, I noticed the discrepancy between your near field beam diameter (83mm) and my nearfield beam diameter (59mm) which is exactly $\sqrt{2}$. The reason is most likely the definition of LightPipes.GaussAperture which states in the docs :param w: 1/e intensity width. This is not the typical textbook definition, I will look into that.

Best regards, Lenny

[ldoyle]

And here is the code:

"""

Use matrix formalism to trace Gaussian beam through a thin lens.

The confusion for me arises when the Rayleigh range of the incoming beam is
not negligible (i.e. not z_R >> f). Then the smallest waist behind the lens
will not be in the focal spot at distance $f$ behind the lens.

* based on ABCD command reference example
* it seems there is a sign error in LightPipes Gauss functions, so here
  reimplement direct calculation (also faster for the problem at hand
  since not dealing with 2D fields but just complex beam parameter)
* values close to problem of Github issue #84

The ABCD-matrix formalism can be used to trace Gaussian beams with the
"complex beam parameter" (see e.g. Wikipedia or optics books)

Created on Thu May 30 23:16:54 2024

@author: Leonard.Doyle
"""

import numpy as np
import matplotlib.pyplot as plt

#import LightPipes as lp # not used for analytic calculation
from LightPipes.units import * # m, mm, um and so on

def M_lens(f):
    """matrix for a thin lens with focal length $f$"""
    return [[1.0,       0.0],
            [-1.0/f,    1.0]]

def M_propagate(z):
    """matrix for free space propagation by distance $z$"""
    return [[1.0,       z],
            [0.0,     1.0]]

def w0_to_zR(w0, lam):
    return np.pi * w0**2 / lam

def zR_to_w0(zR, lam):
    return np.sqrt(zR * lam / np.pi)

def w_to_intensity_FWHM(w):
    """Calculate intensity d_FWHM given field waist w(z). Important:
        assuming w(z) = w0 * sqrt(1+(z/zR)²) defined in E-Field, but
        assuming d_FWHM in intensity (square of E-field)! Otherwise,
        changes by factor of sqrt(2)"""
    return np.sqrt(2*np.log(2)) * w

def w0_to_q(w0, lam):
    """Calculate complex beam parameter $q$ at minimum waist w0 for
    given wavelength $lam$"""
    zR = w0_to_zR(w0, lam)
    q = 1j*zR
    return q

def q_to_w0(q, lam):
    """Return the minimum beam waist w0 of this beam (not the waist at
    current location $z$!)"""
    zR = np.imag(q)
    return zR_to_w0(zR, lam)

def q_to_z(q):
    """Return the current location $z$ with respect to the Gaussian
    waist."""
    z = np.real(q)
    return z

def q_to_w(q, lam):
    """Return the current beam waist $w$ of this beam at the
    current location $z$ (not at minimum waist!)."""
    temp = -(1/q).imag
    w = np.sqrt(lam/(np.pi*temp))
    return w

def q_to_R(q):
    """Return the current radius of curvature $R$ of this beam at the
    current location $z$. At waist, this is $inf$."""
    inv_R = (1/q).real
    R = 1/inv_R if inv_R else np.inf # avoid div by 0 error
    return R

def q_after_M(q, M):
    """Apply matrix M to q and return q' resulting complex beam parameter.
    M must be a tuple or list of shape [[A,B],[C,D]]"""
    ((A,B),(C,D)) = M
    q_out = (A*q + B)/(C*q + D)
    return q_out

lam = 1064*nm
#siz = 1*m # not used in analytical calculation
#N = 1024

w0 = 50*mm
f = 5000*m
# propagate $z1$ to the lens, and $z2$ behind the lens
z1 = 1*m
z2 = 10000*m

print(f'Rayleigh range of this beam: {w0_to_zR(w0, lam):.2f}m')

z0 = -z1
z_before_lens = np.linspace(z0, 0, 101) # odd number to make spacing even
z_after_lens = np.linspace(0, z2, 101)

q0 = w0_to_q(w0, lam)
q1 = np.zeros_like(z_before_lens, dtype=complex)
for i, z in enumerate(z_before_lens):
    # propagate from waist w0 at -z1 to 0 in steps
    q1[i] = q_after_M(q0, M_propagate(z-z0))
q_before_lens = q1[-1] # last step is on lens
q_after_lens = q_after_M(q_before_lens, M_lens(f))
q2 = np.zeros_like(z_after_lens, dtype=complex)
for i, z in enumerate(z_after_lens):
    # start back at lens each time, calculate after distance z2
    q2[i] = q_after_M(q_after_lens, M_propagate(z))

plt.figure(figsize=(10,5))
plt.plot(z_before_lens, q_to_w(q1, lam)/mm), plt.plot(z_after_lens, q_to_w(q2, lam)/mm)
plt.xlabel('Distance in [m], lens at z=0')
plt.ylabel('Waist at $z$ [mm]')

plt.annotate(f'({z0:.0f}|{w0/mm:.1f})', (z0, w0/mm), ha='center', va='bottom')
w_z1 = q_to_w(q1[-1], lam)
plt.annotate(f'({0:.0f}|{w_z1/mm:.1f})', (0, w_z1/mm), ha='center', va='bottom')

i_min = np.argmin(q_to_w(q2, lam))
z_min = z_after_lens[i_min]
w_min = q_to_w(q2[i_min], lam)
plt.annotate(f'w0 ({z_min:.0f}|{w_min/mm:.1e})', (z_min, w_min/mm), ha='center', va='top')

w_z2 = q_to_w(q2[-1], lam)
plt.annotate(f'({z_after_lens[-1]:.0f}|{w_z2/mm:.1f})', (z_after_lens[-1], w_z2/mm), ha='center', va='bottom')

plt.title(f'Waist / FWHM at z={z_after_lens[-1]/m:.2f}m: {w_z2/mm:.2f}mm / {w_to_intensity_FWHM(w_z2)/mm:.2f}mm')

plt.tight_layout() # try to make borders nicer

[SudipKCTUM]

Dear Lenny, Thank you so much for the detailed feedback, especially regarding the divergence of the Gaussian beam because of diffraction when the Rayleigh range of the beam is comparable to the focal length of the thin lens. Your code on tracing the Gaussian beam using ABCD matrix formalism and complex beam parameters made it remarkably clear how this phenomenon works.

Using your recommendation of positioning the waist of the beam not on the lens but 2f before the lens I am able to get the same FWHM after 10 km, which I calculated from the waist in E-field this time to mitigate the discrepancy involving √2. This helped me realize my mistake in my simulation before, where the waist was very close to the lens (1m).

Thanks again for your support.

Sincerely, Sudip

[ldoyle]

Dear Sudip,

thanks for the feedback, you're welcome. One more comment: In your second code, inside the loop you use Field = Fresnel(Field, step_size) iteratively. I recommend starting from the initial field each time, e.g. F = Fresnel(F0, i*step_size). I did see a minor difference in the output between both methods, and I would assume the one which always starts at the same input to be more accurate since it has less chance to collect numerical errors.

More generally, however, you should always assume some numerical imprecision and therefore compare different propagators or change e.g. the step size and see if that changes the outcome, as well as cross-check the results with your expectation.