Airy beams are a class of non-diffracting, self-accelerating optical fields whose propagation behaviors—such as intensity profile evolution, phase dynamics, and auto-focusing—are of great interest in optics. This article uses the Angular Spectrum Method (ASM) to simulate the propagation of one‑dimensional (1D) and two‑dimensional (2D) Airy beams and analyzes the output characteristics at various distances using MATLAB implementations.

Fundamental Properties of Airy Beams

The key features of Airy beams are their non-diffracting nature and self-accelerating trajectory. The mathematical model of an Airy beam is:

\[ U(x,y,0) = \mathrm{Ai}\!\left(\frac{x}{x_0}\right) \exp\!\left(\frac{a x}{x_0}\right) \cdot \mathrm{Ai}\!\left(\frac{y}{y_0}\right) \exp\!\left(\frac{a y}{y_0}\right) \]

where \(\mathrm{Ai}(\cdot)\) is the Airy function, \(x_0, y_0\) are transverse scaling factors controlling the beam width, and \(a\) is the apodization parameter (typically \(0.01-0.1\)) that truncates the infinite energy tail of an ideal Airy beam. The exponential factor \(\exp(a x/x_0)\) suppresses the side lobes and makes the beam carry finite energy.

• Non-diffracting: In theory, the main lobe width remains constant during propagation.
• Self-accelerating: The beam’s peak follows a parabolic trajectory with a deflection angle \(\theta \approx \sqrt{2a}\).
• Self-healing: Even if part of the beam is obscured, the original distribution can re-emerge downstream.

The Angular Spectrum Method (ASM)

ASM is an exact numerical technique for simulating field propagation. It decomposes the initial field into a superposition of plane waves (the angular spectrum), propagates each plane wave by multiplying with a phase factor, and then reconstructs the field via inverse Fourier transform.

One-dimensional casee: The propagation of a 1D Airy beam along \(z\) is described by:

\[ U(x,z) = \mathcal{F}^{-1}\left\{ \mathcal{F}\{U(x,0)\} \cdot \exp\!\left( i\,z\,\sqrt{k_0^2 - k_x^2} \right) \right\} \]

where \(\mathcal{F}\) denotes the Fourier transform, \(k_0 = 2\pi/\lambda\) is the free‑space wavenumber, \(k_x\) is the spatial frequency in the \(x\) direction, and the phase factor \(\exp(i z \sqrt{k_0^2 - k_x^2})\) accounts for propagation over distance \(z\).

Two-dimensional extension: Simply replace the 1D transform with a 2D Fourier transform and use \(k_x, k_y\) coordinates. The method remains the same.

Advantages of ASM:

• Exact: No paraxial or Fresnel approximation is required; valid from near field to far field.
• Flexible: Can handle arbitrary initial field distributions (Airy, Gaussian, vortex, etc.).
• Fast: Uses FFT, complexity \(O(N \log N)\).

1D Airy Beam Propagation

Parameter Setup

• Wavelength: \(\lambda = 532\,\text{nm}\) (green light).
• Scaling factor: \(x_0 = 0.1\,\text{mm}\).
• Apodization: \(a = 0.05\).
• Propagation distances: \(z = 0, 1, 2, 4\,\text{m}\).
• Spatial grid: \(N = 512\) points, total width \(L = 1\,\text{mm}\), step \(\Delta x = L/N \approx 2\,\mu\text{m}\).

MATLAB Code

% 1D Airy beam propagation
lambda = 532e-9;
k0 = 2 * pi / lambda;
scale = 0.1e-3;      % x0
decay = 0.05;        % a
z_vals = [0, 1, 2, 4];
N = 512;
L = 1e-3;
dx = L / N;
x = (-N/2:N/2-1) * dx;

% Initial field
U0 = airy(x / scale) .* exp(decay * x / scale);
U0 = U0 / max(abs(U0));

figure('Position', [100 100 800 600]);
for idx = 1:length(z_vals)
    z = z_vals(idx);
    % Fourier domain
    F = fftshift(fft(U0));
    % Spatial frequencies
    df = 1 / L;
    kx = (-N/2:N/2-1) * df;
    % Propagation factor
    kz = sqrt(k0^2 - kx.^2);
    prop = exp(1i * kz * z);
    % Back to space
    Uz = ifftshift(F .* prop);
    Uz = Uz / max(abs(Uz));
    
    subplot(2, 2, idx);
    plot(x * 1e3, abs(Uz).^2, 'b-', 'LineWidth', 1.5);
    xlabel('x (mm)'); ylabel('Normalized intensity');
    title(sprintf('z = %d m', z));
    grid on; axis tight;
end
sgtitle('1D Airy Beam Propagation');

Results and Discussion

• \(z = 0\) m: The intensity profile shows the typical Airy pattern: a narrow main lobe and rapidly decaying side lobes on one side.
• \(z = 1\) m (near field): The main lobe width remains nearly unchanged, confirming the non‑diffracting nature.
• \(z = 2\) m (mid field): The main lobe begins to move sideways—the self‑acceleration effect becomes visible.
• \(z = 4\) m (far field): The deflection angle increases; side lobes become weaker as energy concentrates into the main lobe.

2D Airy Beam Propagation

Parameter Setup

Same wavelength, scaling \(x_0 = y_0 = 0.1\) mm, apodization \(a = 0.05\), propagation distances \(0, 1, 2, 4\) m. A 2D grid of \(512 \times 512\) points over \(1 \times 1\) mm² is used.

MATLAB Code

% 2D Airy beam propagation
lambda = 532e-9;
k0 = 2 * pi / lambda;
scale = 0.1e-3;
decay = 0.05;
z_vals = [0, 1, 2, 4];
N = 512;
L = 1e-3;
dx = L / N;
x = (-N/2:N/2-1) * dx;
y = (-N/2:N/2-1) * dx;
[X, Y] = meshgrid(x, y);

% Initial field (product of two 1D Airy functions)
U0 = airy(X / scale) .* exp(decay * X / scale) ...
   .* airy(Y / scale) .* exp(decay * Y / scale);
U0 = U0 / max(abs(U0));

figure('Position', [100 100 800 600]);
for idx = 1:length(z_vals)
    z = z_vals(idx);
    F = fftshift(fft2(U0));
    df = 1 / L;
    kx = (-N/2:N/2-1) * df;
    ky = (-N/2:N/2-1) * df;
    [Kx, Ky] = meshgrid(kx, ky);
    Kz = sqrt(k0^2 - Kx.^2 - Ky.^2);
    prop = exp(1i * Kz * z);
    Uz = ifftshift(F .* prop);
    Uz = Uz / max(abs(Uz));
    
    subplot(2, 2, idx);
    imagesc(x * 1e3, y * 1e3, abs(Uz).^2);
    axis square; colorbar; colormap('jet');
    xlabel('x (mm)'); ylabel('y (mm)');
    title(sprintf('z = %d m', z));
end
sgtitle('2D Airy Beam Propagation');

Results and Discussion

• \(z = 0\) m: The intensity is a 2D Airy pattern—a bright central lobe surrounded by fainter side lobes.
• \(z = 1\) m (near field): The central lobe stays mostly circular; side lobes begin to spread.
• \(z = 2\) m (mid field): The main lobe shifts along the diagonal direction (self‑acceleration), with an approximate deflection angle \(\theta \approx 2a\).
• \(z = 4\) m (far field): The deflection becomes more pronounced; side lobes almost vanish as energy consolidates into the main lobe.

Practical Considerations and Optimizations

Choice of the Apodization Parameter \(a\)

• Issue: A ideal Airy beam carries infinite energy (\(\int_{-\infty}^{\infty} |\mathrm{Ai}(x)|^2 dx\) diverges). The parameter \(a\) trunncates the side lobes to make the beam physically realizable.
• Guideline: \(a\) is typically chosen between 0.01 and 0.1. Larger values suppress lobes more but slightly broaden the main lobe.

Sampling and Aliasing

• Issue: The Nyquist criterion requires \(\Delta x < \lambda / (2 \sin\theta_{\max})\), where \(\theta_{\max}\) is the maximum diffraction angle. For an Airy beam with \(a=0.05\), \(\theta_{\max} \approx 30^\circ\), so \(\Delta x < 532\,\text{nm}\). Our grid step \(\Delta x \approx 2\,\mu\text{m}\) satisfies this.
• Optimization: In practice, one can increase \(N\) or reduce \(L\) if higher spatial frequencies are present.

Computational Efficiency

• Use MATLAB’s built‑in fft/fft2 for fast computation.
• For multiple propagation distances, parfor can parallelize the loop (requires Parallel Computing Toolbox).

Summary

Using the Angular Spectrum Method, we have simulated the propagation of both 1D and 2D Airy beams over different distances. The results verify the key properties:

• The main lobe of a 1D Airy beam maintains its width (non‑diffraction) while its centroid follows a parabolic path (self‑acceleration).
• The 2D Airy beam evolves into a moving spot that shifts diagonally, with side lobes fading as energy coalesces into the central lobe.

These features make Airy beams promising for applications such as free‑space optical communications, super‑resolution imaging, and optical trapping. Future work may involve phase modulation (e.g., using cotangent or cosine masks) to further tailor the propagation behavior.