The following visualisation shows Hermite-Gaussian (HG) mode solutions to the paraxial wave equation. The field amplitude at any plane along the beam axis is proportional to \[U(x, y) \propto H_n H_m \exp(-i (n + m) \arctan(z/z_R)) \] where \(H_n\) are the Hermite polynomials, \(z\) is the distance from the beam focus, \(z_R\) is the Rayleigh range and \(n\) and \(m\) are the mode order parameters.

The left panel shows the beam in 3-D space. The top panel shows the amplitude/intensity profile of the beam and the bottom panel shows the phase of the electric field.

Code Samples

Matlab

function hgbeam

    n = 3;                  % Beam order
    m = 1;                  % Beam order
    w0 = 2.0;               % Beam waist
    k = 2*pi/532.0e-9;      % Wavenumber of light

    zR = k*w0^2/2;      % Calculate the Rayleigh range

    % Setup the cartesian grid for the plot at plane z
    z = 0.0;
    [xx, yy] = meshgrid(linspace(-5, 5), linspace(-5, 5));

    U00 = 1/(1 + 1i*z/zR) .* exp(-(xx.^2 + yy.^2)/w0^2./(1 + 1i*z/zR));
    Hn = hermiteH(n, xx);
    Hm = hermiteH(m, yy);

    U = U00.*Hn.*Hm.*exp(-1i*(n + m).*atan(z/zR));

    figure;
    subplot(1, 2, 1);
    imagesc(abs(U).^2);
    title('Intensity');
    subplot(1, 2, 2);
    imagesc(angle(U));
    title('Phase');

end

% hermiteH function was introduced in R2017a (but we are sing R2014a)
function ret = hermiteH(n, x)
    if n == 0
        ret = 1;
    elseif n == 1
        ret = 2*x;
    else
        ret = 2*x.*hermiteH(n-1, x) - 2*(n - 1)*hermiteH(n - 2, x);
    end
end

Python

import numpy as np
import matplotlib.pyplot as plt

n = 3;                  # Beam order
m = 1;                  # Beam order
w0 = 2.0;               # Beam waist
k = 2*np.pi/532.0e-9;   # Wavenumber of light

zR = k*w0**2.0/2;       # Calculate the Rayleigh range

# Setup the cartesian grid for the plot at plane z
z = 0.0;
[xx, yy] = np.meshgrid(np.linspace(-5, 5), np.linspace(-5, 5));

def hermiteH(n, x):
    if n == 0:
        return 1;
    elif n == 1:
        return 2*x;
    else:
        return 2*x*hermiteH(n-1, x) - 2*(n - 1)*hermiteH(n - 2, x);

U00 = 1.0/(1 + 1j*z/zR)*np.exp(-(xx**2 + yy**2)/w0**2/(1 + 1j*z/zR));
Hn = hermiteH(n, xx);
Hm = hermiteH(m, yy);

U = U00*Hn*Hm*np.exp(-1j*(n + m)*np.arctan(z/zR));

plt.figure()
plt.title('Intensity')
plt.pcolor(abs(U)**2);
plt.axis('equal')

plt.figure()
plt.title('Phase')
plt.pcolor(np.angle(U)**2);
plt.axis('equal')

plt.show()

Hermite-Gaussian Beams

Code Samples