Laser Beam Laguerre-Gaussian TEM modes¶
Author: Calvin Z He
The electric field plotted is of the form ($z$ and phase offsets ignored)$^1$: $$E_{pl}(r,\phi,z=0) = E_0\left(\frac{2r^2}{w_0^2}\right)^{|l|/2} L_p^{|l|}\left(\frac{2r^2}{w_0^2}\right)\exp\left(-\frac{r^2}{w_0^2}\right)\cos(l\phi)$$
The generalized Laguerre Polynomial in sum notation$^2$: $$L_p^l(\rho) = \sum_{m=0}^p (-1)^m \frac{(p+l)!}{(p-m)!(l+m)!m!}\rho^m$$
The generalized Laguerre Polynomial can also be generated using scipy's genlaguerre function.
References:
import numpy as np
from scipy.special import factorial, genlaguerre
from matplotlib import pyplot as plt
from matplotlib import rcParams
rcParams['text.usetex'] = True
x = np.linspace(-100,100,101)
y = np.linspace(-100,100,101)
X,Y = np.meshgrid(x,y)
R = np.sqrt(X**2+Y**2)
Phi = np.atan2(Y,X)
pmax = 6
lmax = 3
w0 = 10
for p in range(pmax):
for l in range(lmax):
Clp = np.sqrt(2*factorial(p)/(np.pi*factorial(p+np.abs(l))))
#Lpl = 0
#for m in range(p+1):
# Lpl += (-1)**m*factorial(p+l)/(factorial(p-m)*factorial(l+m)*factorial(m))*(2*R**2/w0**2)**m
Lpl = genlaguerre(p,l)(2*R**2/w0**2)
Epl = np.cos(Phi*l)*np.exp(-R**2/w0**2)*(np.sqrt(2)*R/w0)**np.abs(l)*Clp*Lpl/w0
if l==0:
Epl_row = [Epl]
else:
Epl_row = np.append(Epl_row, [Epl],axis=0)
if p==0:
Epl_dat = [Epl_row]
else:
Epl_dat = np.append(Epl_dat, [Epl_row],axis=0)
print(np.shape(Epl_dat))
(6, 3, 101, 101)
%matplotlib inline
plotextent = [min(x), max(x), min(y), max(y)]
fig, ax = plt.subplots(figsize=(12,18),ncols=lmax,nrows=pmax)
for p in range(pmax):
for l in range(lmax):
axc = ax[p,l]
dat = Epl_dat[p,l]
im = axc.imshow(dat,cmap='seismic',vmax = 0.08,vmin=-0.08,norm='linear',extent=plotextent)
plt.colorbar(im, ax=axc)
axc.set_xlabel('$x$')
axc.set_ylabel('$y$')
axc.set_title(f'$(p,l)=({p},{l})$')
plt.tight_layout()
plt.show()
D4σ and M-squared¶
The D4σ is a conventional way of quantifying an arbitrary beam's size. Is is equal to four times the variance $\sigma$ calculated in the $x$ and $y$ directions via: $$\sigma_x = \sqrt{\frac{\iint x^2 I(x,y)dxdy}{\iint I(x,y)dxdy}}$$ $$\sigma_y = \sqrt{\frac{\iint y^2 I(x,y)dxdy}{\iint I(x,y)dxdy}}$$ So that $$D4\sigma_x = 4\sigma_x = 2w_x$$ $$D4\sigma_y = 4\sigma_y = 2w_y$$ where $w_x$ and $w_y$ are the beam radii in the $x$ and $y$ directions respectively.
Since $D4\sigma_x$ and $D4\sigma_y$ is not necessarily the same in general, I choose to define an averaged quantity $$D4\sigma_{avg} = (D4\sigma_x+D4\sigma_y)/2$$ in order to quantify the "overall" width of the profile, given that it is reasonably square.
The $M^2$ value of a beam quantifies deviation of the profile from the diffraction limited beam. It is implied from the equation for the divergence half angle: $$\theta_D = M^2 \frac{\lambda}{\pi w}$$ that compared to the diffraction limited case: $$\theta_D = \frac{\lambda}{\pi w_0}$$ that the non-diffraction limited beam radius $w$ is related to the diffraction limited beam radius $w_0$ via: $$w=M^2w_0$$ This means the $M^2$ value for each TEM mode can be "measured" based on the ratio of their widths to that of the TEM00 mode: $$M^2 = \frac{D4\sigma_{avg}(p,l)}{D4\sigma_{avg}(0,0)}$$
One contradiction I've ran into, however, is that the $M^2$ value of the Laguerre-Gaussian modes is given by$^1$: $M^2 = 2p+|l|+1$, however the $D4\sigma$ measurements show a weaker linear dependence on $p$ and $|l|$.
for p in range(pmax):
for l in range(lmax):
dat = Epl_dat[p,l]
D4sigma_x = np.sqrt(np.sum(X**2*dat**2)/np.sum(dat**2))*4
D4sigma_y = np.sqrt(np.sum(Y**2*dat**2)/np.sum(dat**2))*4
D4sigma_avg = np.mean([D4sigma_x,D4sigma_y])
if l==0:
D4sigma_x_row = [D4sigma_x]
D4sigma_y_row = [D4sigma_y]
D4sigma_avg_row = [D4sigma_avg]
else:
D4sigma_x_row = np.append(D4sigma_x_row,[D4sigma_x],axis=0)
D4sigma_y_row = np.append(D4sigma_y_row,[D4sigma_y],axis=0)
D4sigma_avg_row = np.append(D4sigma_avg_row,[D4sigma_avg],axis=0)
if p==0:
D4sigma_x_list = [D4sigma_x_row]
D4sigma_y_list = [D4sigma_y_row]
D4sigma_avg_list = [D4sigma_avg_row]
else:
D4sigma_x_list = np.append(D4sigma_x_list, [D4sigma_x_row],axis=0)
D4sigma_y_list = np.append(D4sigma_y_list, [D4sigma_y_row],axis=0)
D4sigma_avg_list = np.append(D4sigma_avg_list, [D4sigma_avg_row],axis=0)
M2_measured = D4sigma_avg_list/w0/2
for l in range(lmax):
plt.plot(M2_measured[:,l],'.',label=f'Measured: $l={l}$')
for l in range(lmax):
plt.plot(2*np.array(range(pmax)) + l +1,'o',label=f'Theoretical: $l={l}$',alpha=0.5)
plt.xlabel('$p$')
plt.ylabel('$M^2$')
plt.legend()
plt.show()
The normalization coefficient¶
According to the wikipedia section on Laguerre-Gaussian profiles$^3$, the normalization coefficient necessary for the TEM mode to integrate to unity over all space is: $$C_p^l = \sqrt{\frac{2p!}{\pi(p+|l|)!}}$$ However, when using this coefficient, I find that the integral is unity only for $l=0$ modes, whereas it is $0.5$ for $l\neq0$. Perhaps the positive and negative values for a given $|l|$ are supposed to be added together to achieve unity?
dx = x[1]-x[0]
dy = y[1]-y[0]
for p in range(pmax):
for l in range(lmax):
dat = Epl_dat[p,l]
print(f'(p,l)=({p},{l}): {np.sum(dat**2*dx*dy):.2f}')
(p,l)=(0,0): 1.00 (p,l)=(0,1): 0.50 (p,l)=(0,2): 0.50 (p,l)=(1,0): 1.00 (p,l)=(1,1): 0.50 (p,l)=(1,2): 0.50 (p,l)=(2,0): 1.00 (p,l)=(2,1): 0.50 (p,l)=(2,2): 0.50 (p,l)=(3,0): 1.00 (p,l)=(3,1): 0.50 (p,l)=(3,2): 0.50 (p,l)=(4,0): 1.00 (p,l)=(4,1): 0.50 (p,l)=(4,2): 0.50 (p,l)=(5,0): 1.00 (p,l)=(5,1): 0.50 (p,l)=(5,2): 0.50