Why a paraxial source cannot focus tightly

The Gaussian beam is a solution of the paraxial wave equation, not of the full Helmholtz equation. It is derived by writing the field as a slowly varying envelope on a carrier, \(E = u(x,y,z)\,e^{ikz}\), and then discarding \(\partial^2 u/\partial z^2\). In the plane-wave picture that approximation is the expansion

\[k_z = \sqrt{k^2 - k_\perp^2} \;\approx\; k - \frac{k_\perp^2}{2k},\]

which is accurate only while \(k_\perp \ll k\), i.e. while every plane wave in the beam travels within a small angle of the axis.

In a Gaussian beam’s geometry, waist and far-field divergence are locked together, \(\theta = \lambda/\pi w_0\), so the beam’s angular content is fixed by the requested waist:

\[w_0 = \frac{\lambda}{\pi\,\mathrm{NA}} \quad\Longleftrightarrow\quad \theta = \mathrm{NA}.\]

For instance, at NA 0.9 the paraxial formula assigns the beam a divergence of 0.9 radian, about 52 degrees, while the true cone half-angle for NA 0.9 is \(\arcsin(0.9) \approx 64\) degrees. Two things then go wrong at once:

1. The propagation phase is wrong. Writing \(x = (k_\perp/k)^2\), the exact and paraxial axial wavenumbers are \(k_z/k = \sqrt{1-x}\) and \(1 - x/2\), and \(\sqrt{1-x} < 1 - x/2\) for every \(x > 0\). The paraxial expansion therefore overestimates \(k_z\), by more and more as \(k_\perp\) grows: it has each off-axis plane wave accumulating axial phase faster than it really does. At NA 0.9 the outermost component sits at \(k_\perp/k = 0.9\), where the true \(k_z/k = \sqrt{0.19} \approx 0.436\) against a paraxial \(0.595\) — 37% too large, which over the 2 um from source plane to focus is a phase error of about 4 radians for that component. The ansatz sets the wavefront curvature on the source plane so that every component would arrive in phase at the focus if the paraxial \(k_z\) held; since the outer components in fact accumulate less axial phase than that, they arrive mistimed and stop adding coherently. The paraxial form also assigns propagating character to spatial frequencies with \(k_\perp > k\), which in reality are evanescent and cannot reach the focus at all.

2. The polarization is wrong. The paraxial Gaussian is purely transverse. But \(\nabla\cdot\mathbf{E} = 0\) forces a longitudinal component as soon as the transverse field varies on the scale of a wavelength, \(E_z \sim (i/k)\,\partial E_x/\partial x\). A real tightly focused beam always has \(E_z\); the paraxial ansatz has none.

The result is that, beyond a certain NA, the field a paraxial source launches simply does not converge any further, no matter how small a waist_radius it is given.

td.ThinLensBeam makes no paraxial assumption, instead following the fully vectorial method of the following papers:

1. M. Mansuripur, “Distribution of light at and near the focus of high-numerical-aperture objectives,” J. Opt. Soc. Am. A 3, 2086–2093 (1986).

2. M. Mansuripur, “Certain computational aspects of vector diffraction problems,” J. Opt. Soc. Am. A 6, 786–805 (1989).

3. M. Mansuripur, “Distribution of light at and near the focus of high-numerical-aperture objectives: erratum, Certain computational aspects of vector diffraction problems: erratum,” J. Opt. Soc. Am. A 10, 382–383 (1993).

What to expect from the result

Two consequences of the above change how the numbers should be read.

The thin-lens focus is an Airy pattern, not a Gaussian. With fill_lens=True the pupil is uniformly illuminated, and a uniformly illuminated circular pupil focuses to an Airy core with rings. That is a different function from a Gaussian mode, so it has a different width for the same diffraction limit: fitting a Gaussian to an Airy core gives \(w\,\mathrm{NA}/\lambda \approx 0.421\) where a true Gaussian mode gives \(1/\pi \approx 0.318\).

The thin-lens focus is elliptical at high NA. The polarization rotation that supplies \(E_z\) also makes the spot wider along the polarization direction than across it. At NA 0.9 the two transverse axes differ by about 35%, so a single “waist” number is no longer the whole story; the sweep below fits both axes rather than averaging them together.

# Standard python imports
import matplotlib.pyplot as plt
import numpy as np

# Tidy3D import
import tidy3d as td
from tidy3d import web

# Import base object
from tidy3d.components.beam import BeamProfile

td.config.logging.level = "ERROR"

wvl0 = 0.5
freq0 = td.C_0 / wvl0

beam_NA = np.linspace(0.05, 0.9, 20)
expected_waist = wvl0 / np.pi / beam_NA

# 10 x 10 um source plane, injected ``focus_distance`` upstream of z = 0 so the requested waist
# lands on the focal-plane monitor. All three FDTD sources radiate on this same plane, which for a
# source is a real aperture -- but see the note below: the same aperture does not cost the three
# sources the same amount of light.
source_size = 20 * wvl0
focus_distance = 2.0

# FDTD domain and run time. 16 um transverse keeps the low-NA beams (w0 = 3.2 um at NA 0.05)
# clear of the boundary: the worst case over the sweep leaves 3e-6 of the peak intensity at the
# domain edge, so truncation is negligible.
sim_size = (16, 16, 15)
run_time = 60 / freq0

The building of the Quasi-Gaussian beam profile is done with the BeamProfile object. Details can be found in this notebook.

def QuasiGaussian(
    size,
    source_time,
    center=(0, 0, 0),
    direction="+",
    pol_angle=0,
    angle_phi=0,
    angle_theta=0,
    w0=1,
    waist_distance=0,
):
    """Quasi-Gaussian beam, injected as a ``td.CustomFieldSource``.

    The analytic profile carries the leading non-paraxial correction. Sampling it on the injection
    plane and handing the tangential field to a custom source is what lets FDTD propagate it.
    """
    freqs = [source_time.freq0]

    class QuasiGaussian_obj(BeamProfile):
        """Component for constructing quasi-Gaussian beam data. The normal direction is
        implicitly defined by the ``size`` parameter.
        """

        def scalar_field(self, points, background_n):
            """Scalar field for quasi-Gaussian beam.
            Scalar field corresponding to the analytic beam in coordinate system such that the
            propagation direction is z and the ``E``-field is entirely ``x``-polarized. The field is
            computed on an unstructured array ``points`` of shape ``(3, ...)``.
            """
            x, y, z = points
            z = z - waist_distance
            r = np.sqrt(x**2 + y**2)

            k = 2 * np.pi * np.array(freqs) / td.C_0 * background_n
            zR = np.real(w0**2 * k / 2)
            wz = w0 * np.sqrt(1 + z**2 / zR**2)
            inv_Rz = z / (z**2 + zR**2)
            Fpp = np.sqrt(1 + r**2 * inv_Rz**2)
            phi0 = np.arctan2(z, zR)

            phase = -((r / Fpp) ** 2) / wz**2 - 1j * k * z + 1j * phi0
            if np.all(inv_Rz) != 0:
                phase -= 1j * k * (Fpp - 1) / inv_Rz

            E = w0 / wz * 1 / Fpp**2 * np.exp(phase)

            return E

    profile = QuasiGaussian_obj(
        center=center,
        size=size,
        freqs=freqs,
        angle_theta=angle_theta,
        angle_phi=angle_phi,
        pol_angle=pol_angle,
        direction=direction,
    )
    field_data = profile.field_data
    field_dataset = td.FieldDataset(
        Ex=field_data.Ex,
        Ey=field_data.Ey,
        Ez=field_data.Ez,
        Hx=field_data.Hx,
        Hy=field_data.Hy,
        Hz=field_data.Hz,
    )
    return td.CustomFieldSource(
        center=center, size=size, source_time=source_time, field_dataset=field_dataset
    )

Why We Simulate All Three

A beam profile object evaluates a closed-form expression at whatever points you hand it. For a paraxial Gaussian that expression is an ansatz: it asserts a waist of waist_radius at the focal plane. Evaluating it there returns what the formula claims, which is not the same thing as what the source delivers — and, per the discussion above, the two part company exactly where the paraxial approximation stops holding.

So the measurement is made the way an experiment would make it. Each source is injected into an empty FDTD domain on a plane focus_distance upstream of its nominal focus, and a field monitor at \(z=0\) records the spot that arrives. FDTD solves the full Maxwell equations, so nothing about the propagation is paraxial: the source hands the solver a tangential field on its injection plane, and every spatial frequency in it then travels with its true \(k_z\), with \(k_\perp > k\) components evanescent rather than propagating. Where the paraxial approximation is good this reproduces \(\lambda/\pi\mathrm{NA}\) and confirms the source works as advertised. Where it is not, the simulated field shows what actually turns up.

Three simulations are built per NA, identical except for the source:

• td.GaussianBeam, the paraxial reference,
• the quasi-Gaussian profile above, wrapped in a td.CustomFieldSource,
• td.ThinLensBeam, the vectorial thin lens.

Putting the thin lens through the same solver as the other two is what makes the comparison apples-to-apples. All three share the injection plane, the domain, the grid, the run time and the monitor, and are fit over the same footprint; the only difference between the three runs is the field on the injection plane. The domain is homogeneous vacuum with no structures, and its transverse extent is large enough that the beam is negligible at the boundary for every NA in the sweep, so the fitted spot is set by the source and the physics rather than by the simulation aperture.

def create_beam_sims(w0, na):
    """Three otherwise-identical vacuum simulations of a beam focused to ``w0`` / ``na``.

    Returns ``(QG_sim, G_sim, TL_sim)``: the quasi-Gaussian custom source, ``td.GaussianBeam`` and
    ``td.ThinLensBeam``. Domain, grid, run time and monitors are shared, so the fitted waists are
    directly comparable.
    """
    QG_size = (source_size, source_size, 0)
    QG_source_time = td.GaussianPulse(freq0=freq0, fwidth=freq0 / 10)
    # Short propagation distance from source to monitor keeps the beam radius
    # small relative to the domain even at high NA (fast-diverging, small w0).
    QG_waist_distance = -focus_distance
    QG_center = (0, 0, QG_waist_distance)

    QG_source = QuasiGaussian(
        center=QG_center,
        size=QG_size,
        source_time=QG_source_time,
        w0=w0,
        waist_distance=QG_waist_distance,
    )

    G_source = td.GaussianBeam(
        center=QG_center,
        size=QG_size,
        source_time=QG_source_time,
        direction="+",
        waist_radius=w0,
        waist_distance=QG_waist_distance,
    )

    # The thin lens is specified by NA rather than by a waist, and fills its pupil uniformly,
    # so it aims at the Airy limit rather than the Gaussian-mode limit.
    TL_source = td.ThinLensBeam(
        center=QG_center,
        size=QG_size,
        source_time=QG_source_time,
        direction="+",
        numerical_aperture=na,
        waist_distance=QG_waist_distance,
        fill_lens=True,
        num_plane_waves=101,
    )

    field = td.FieldMonitor(size=(td.inf, td.inf, 0), name="field", freqs=[freq0])

    QG_sim = td.Simulation(
        size=sim_size,
        sources=[QG_source],
        structures=[],
        monitors=[field],
        run_time=run_time,
    )

    G_sim = QG_sim.updated_copy(sources=[G_source])
    TL_sim = QG_sim.updated_copy(sources=[TL_source])

    return QG_sim, G_sim, TL_sim

Run the NA Sweep Batch

Three simulations per NA, submitted as a single batch. Tasks are keyed by sweep index — "q{i}" for the quasi-Gaussian source, "{i}" for td.GaussianBeam and "t{i}" for td.ThinLensBeam — so the returned BatchData indexes straight back against beam_NA.

That is 60 vacuum simulations on a 16 x 16 x 15 um domain at the default auto grid (dl = 0.05 um, 10 steps per wavelength), so it costs real FlexCredits — call gauss_benchmark_batch.estimate_cost() first if that matters.

gauss_benchmark_sims = {}

for i, (w, na) in enumerate(zip(expected_waist, beam_NA)):
    qg, g, tl = create_beam_sims(float(w), float(na))
    gauss_benchmark_sims[f"q{i}"] = qg  # quasi-Gaussian custom source
    gauss_benchmark_sims[str(i)] = g  # td.GaussianBeam
    gauss_benchmark_sims[f"t{i}"] = tl  # td.ThinLensBeam

gauss_benchmark_batch = web.Batch(simulations=gauss_benchmark_sims)

gauss_benchmark_data = gauss_benchmark_batch.run()

20:31:16 UTC Started working on Batch containing 60 tasks.                      

20:31:17 UTC Maximum FlexCredit cost: 1.500 for the whole batch.                

             Use 'Batch.real_cost()' to get the billed FlexCredit cost after    
             completion.                                                        

20:31:25 UTC Batch complete.                                                    

def focal_intensity(task_name):
    """Focal-plane intensity from one task, as a 2-D ``DataArray`` over ``(x, y)``.

    This is the FDTD counterpart of "the delivered spot": whatever the source launched, after the
    solver has carried it ``focus_distance`` downstream to z = 0.
    """
    return gauss_benchmark_data[task_name]["field"].intensity.isel(f=0).squeeze(drop=True)

Fit Beam Waist from Intensity Profiles

The following section defines a Gaussian fitting helper and uses x/y slices through peak intensity to estimate the waist. Every curve below is fit with this same helper over the same footprint, so the three methods are directly comparable.

def gaussian(coords, amp, center, wz):
    wz = float(wz)
    g_values = amp * np.exp(-2*(center - coords)**2/wz**2)
    return g_values

from scipy import optimize
from scipy.optimize import OptimizeWarning
import warnings

warnings.simplefilter('ignore', OptimizeWarning)

def estimate_waist(intensity_data, expected_waist_radius, plot=False):
    indices = intensity_data.argmax(dim=['x', 'y'])

    data_x = np.squeeze(intensity_data.isel(y=indices['y']).values)
    data_y = np.squeeze(intensity_data.isel(x=indices['x']).values)
    x, y = intensity_data.x.values, intensity_data.y.values
    expected_amp = float(np.max(intensity_data))

    x0 = float(np.atleast_1d(x[indices['x']])[0])
    y0 = float(np.atleast_1d(y[indices['y']])[0])
    p0_x = [expected_amp, x0, expected_waist_radius]
    p0_y = [expected_amp, y0, expected_waist_radius]

    popt_x, pcov_x = optimize.curve_fit(gaussian, x, data_x, p0=p0_x)
    popt_y, pcov_y = optimize.curve_fit(gaussian, y, data_y, p0=p0_y)

    if plot:
        fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(10, 4))

        ax1.plot(x, data_x, label="Data")
        ax1.plot(x, gaussian(x, popt_x[0], popt_x[1], popt_x[2]), label="Fit")
        ax1.plot([popt_x[1], popt_x[1] + popt_x[2] / 2], [popt_x[0] / 2, popt_x[0] / 2])
        ax1.legend()
        ax1.set_title("X data")

        ax2.plot(y, data_y, label="Data")
        ax2.plot(y, gaussian(y, popt_y[0], popt_y[1], popt_y[2]), label="Fit")
        ax2.plot([popt_y[1], popt_y[1] + popt_y[2] / 2], [popt_y[0] / 2, popt_y[0] / 2])
        ax2.legend()
        ax2.set_title("Y data")

        plt.tight_layout()
        plt.show()

    ave_width = popt_x[2] / 2 + popt_y[2] / 2
    return ave_width

def estimate_waist_xy(intensity_data, expected_waist_radius):
    """Same fit as ``estimate_waist``, but returns the two transverse axes separately.

    At high NA the thin-lens focal spot is elliptical (broader along the polarization direction), so
    averaging the axes hides a real vectorial effect.
    """
    indices = intensity_data.argmax(dim=["x", "y"])

    data_x = np.squeeze(intensity_data.isel(y=indices["y"]).values)
    data_y = np.squeeze(intensity_data.isel(x=indices["x"]).values)
    x, y = intensity_data.x.values, intensity_data.y.values
    expected_amp = float(np.max(intensity_data))

    x0 = float(np.atleast_1d(x[indices["x"]])[0])
    y0 = float(np.atleast_1d(y[indices["y"]])[0])

    popt_x, _ = optimize.curve_fit(gaussian, x, data_x, p0=[expected_amp, x0, expected_waist_radius])
    popt_y, _ = optimize.curve_fit(gaussian, y, data_y, p0=[expected_amp, y0, expected_waist_radius])

    return abs(popt_x[2]), abs(popt_y[2])

The Two Diffraction Limits

The paraxial sources and the thin lens aim at two different diffraction limits, so each needs its own reference line before anything can be compared.

For a Gaussian mode the coefficient is \(1/\pi\) by definition, since \(w = \lambda/\pi\mathrm{NA}\). For a uniformly filled circular pupil the focus is the Airy pattern \([2J_1(v)/v]^2\) with \(v = k\,\mathrm{NA}\,r\), and the cell below obtains its coefficient by fitting the same Gaussian model used everywhere else to that profile. Working in the normalized coordinate \(v\) makes the answer independent of NA, so it is a single number that serves the whole sweep.

from scipy import special

v = np.linspace(-40, 40, 80001)
v_safe = np.where(v == 0, 1e-12, v)
airy_profile = (2 * special.j1(np.abs(v_safe)) / np.abs(v_safe)) ** 2

popt_airy, _ = optimize.curve_fit(
    gaussian, v, airy_profile, p0=[1.0, 0.0, 2.6], bounds=([0, -1, 1e-3], [np.inf, 1, 50])
)

# w = V / (k NA)  =>  w NA / lambda = V / (2 pi)
AIRY_COEFF = abs(popt_airy[2]) / (2 * np.pi)
GAUSS_COEFF = 1 / np.pi

print(f"Gaussian-mode target : w NA / lambda = {GAUSS_COEFF:.4f}")
print(f"Filled-pupil target  : w NA / lambda = {AIRY_COEFF:.4f}")
print(f"ratio                : {AIRY_COEFF / GAUSS_COEFF:.3f}x")

Gaussian-mode target : w NA / lambda = 0.3183
Filled-pupil target  : w NA / lambda = 0.4207
ratio                : 1.322x

Sweep NA

Every point on both panels below is now read back from an FDTD focal-plane monitor, and all three are fit by estimate_waist, which collapses the two transverse axes into a single averaged waist. Same solver, same domain, same aperture, same estimator, same footprint — the curves differ only because the injected fields differ.

The thin lens is additionally fit axis-by-axis with estimate_waist_xy, because its focus is genuinely elliptical at high NA and the averaged number hides that. The average of those two axes is identical to what estimate_waist returns, so the grey band in the right panel brackets the same marker it is drawn around.

g_beam_estimated_waist = []
qg_beam_estimated_waist = []
tl_beam_estimated_waist = []
tl_waist_x = []
tl_waist_y = []

for i, e in enumerate(expected_waist):
    g_beam_estimated_waist.append(estimate_waist(focal_intensity(str(i)), e))
    qg_beam_estimated_waist.append(estimate_waist(focal_intensity(f"q{i}"), e))

    tl_intensity = focal_intensity(f"t{i}")
    tl_beam_estimated_waist.append(estimate_waist(tl_intensity, e))
    # both axes kept for the thin lens, so the high-NA ellipticity stays visible
    wx, wy = estimate_waist_xy(tl_intensity, e)
    tl_waist_x.append(wx)
    tl_waist_y.append(wy)

g_beam_estimated_waist = np.array(g_beam_estimated_waist)
qg_beam_estimated_waist = np.array(qg_beam_estimated_waist)
tl_beam_estimated_waist = np.array(tl_beam_estimated_waist)
tl_waist_x = np.array(tl_waist_x)
tl_waist_y = np.array(tl_waist_y)

Conclusion

A paraxial Gaussian source cannot deliver waists at high NA. Both paraxial sources follow their \(1/\pi\) target up to NA \(\approx\) 0.3 and then leave it, and in absolute terms the delivered waist simply stops shrinking: td.GaussianBeam bottoms out near 0.38 um and the quasi-Gaussian near 0.27 um, however small a waist_radius is requested. At NA 0.9 those are 2.2x and 1.5x the requested waist.

The thin lens keeps focusing. td.ThinLensBeam tracks its filled-pupil target across the entire sweep, staying within roughly 10% of \(0.421\,\lambda/\mathrm{NA}\) where the paraxial sources drift by 50–120%. In absolute terms it is still shrinking at NA 0.9, where both paraxial sources have long since stalled.

How to Read the Plots

The left panel is in absolute units. The right panel is the same data normalized as \(w\,\mathrm{NA}/\lambda\): a horizontal line there means the spot is still shrinking in proportion to \(1/\mathrm{NA}\), while a rising curve means focusing has stalled.

There are two dashed targets, one per family of source, because the two are aiming at different diffraction limits:

• Gaussian-mode target, \(\lambda/\pi\mathrm{NA}\), or \(1/\pi = 0.318\) normalized. This is what td.GaussianBeam and the quasi-Gaussian source are asking for, and where they should sit if the paraxial description holds.
• Filled-pupil target, \(0.421\,\lambda/\mathrm{NA}\). This is the Airy core of a uniformly filled circular pupil, which is what td.ThinLensBeam should hit. It lies above the Gaussian-mode line because an Airy core and a Gaussian mode have different widths at the same diffraction limit — not because it is a worse focus.

The grey band around the thin-lens points is the x–y spread. Its two edges are the waist fitted along \(x\), the polarization direction, and along \(y\), across it. At low NA the two are equal and the band is invisible. As NA grows the polarization rotation widens the focus along \(x\) and tightens it along \(y\), so the band opens out; at NA 0.9 the two axes differ by about 35%. The black markers are the average of the two edges, which is what a single-waist fit would report.

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 5.6))

# ---- absolute units ----
ax1.plot(beam_NA, expected_waist, "--", color="steelblue", lw=1.6,
         label="Gaussian-mode target  λ/πNA")
ax1.plot(beam_NA, AIRY_COEFF * wvl0 / beam_NA, "--", color="dimgray", lw=1.6,
         label=f"Filled-pupil target  {AIRY_COEFF:.3f}λ/NA")
ax1.scatter(beam_NA, g_beam_estimated_waist, marker="s", color="green", s=42,
            label="Gaussian source (FDTD)", zorder=3)
ax1.scatter(beam_NA, qg_beam_estimated_waist, marker="H", color="red", s=42,
            label="Quasi-Gaussian source (FDTD)", zorder=3)
ax1.scatter(beam_NA, tl_beam_estimated_waist, marker="o", color="black", s=22,
            label="Thin lens (FDTD)", zorder=3)
ax1.set_xlabel("Beam NA")
ax1.set_ylabel("Delivered waist (µm)")
ax1.set_title("Delivered focal waist")
ax1.set_ylim(0, 3.6)
ax1.grid(alpha=0.3)
ax1.legend(fontsize=8)

# ---- normalized: flat means "still focusing as 1/NA" ----
ax2.axhline(GAUSS_COEFF, ls="--", color="steelblue", lw=1.6, label="Gaussian-mode target  1/π")
ax2.axhline(AIRY_COEFF, ls="--", color="dimgray", lw=1.6,
            label=f"Filled-pupil target  {AIRY_COEFF:.3f}")
ax2.fill_between(beam_NA,
                 tl_waist_y * beam_NA / wvl0,
                 tl_waist_x * beam_NA / wvl0,
                 color="black", alpha=0.15, label="Thin lens, x–y spread (vectorial)")
ax2.scatter(beam_NA, g_beam_estimated_waist * beam_NA / wvl0, marker="s", color="green", s=42,
            label="Gaussian source (FDTD)")
ax2.scatter(beam_NA, qg_beam_estimated_waist * beam_NA / wvl0, marker="H", color="red", s=42,
            label="Quasi-Gaussian source (FDTD)")
ax2.scatter(beam_NA, tl_beam_estimated_waist * beam_NA / wvl0, marker="o", color="black", s=22,
            label="Thin lens (FDTD)")
ax2.set_xlabel("Beam NA")
ax2.set_ylabel("w · NA / λ   (flat = focusing as 1/NA)")
ax2.set_title("Normalized: does it keep focusing?")
ax2.grid(alpha=0.3)
ax2.legend(fontsize=8)

plt.tight_layout()
plt.show()