Potential inconsistency in Zernike masking and beam waist calculation

[zzpwahaha]

Hi, I'm currently working at commit 76bb028cc9b5a621aa7a439940f041875c65b82d (Jan 3, 2025).

1. Zernike Masking in zernike_sum

In the file holography.toolbox.phase.py, the function zernike_sum() contains a potentially confusing or incorrect aperture masking operation around line 889 with the code:

mask = np.square(x_grid * x_scale) + np.square(y_grid * y_scale) <= 1

The x_scale and y_scale parameters used here are generated by the zernike_aperture() function. These values correspond to the reciprocal of the input beam amplitude diameter (in appropriate units), assuming the SLM has already undergone wavefront calibration and the input beam amplitude was characterized in the process.

This implies that the beam, at its 1/e² intensity contour, spans a circular region with diameter = 1, radius = 0.5, and centered at the origin (0, 0).

However, in zernike_sum(), the aperture mask appears to compute:

mask = np.square(x * x_scale) + np.square(y * y_scale) <= 1

This checks whether each pixel is within a unit radius circle (i.e., compares to 1 = (diameter)²) rather than 0.5² = 0.25 = (radius)², which seems inconsistent with the assumption that the beam radius is 0.5. Should this threshold be:

mask = np.square(x * x_scale) + np.square(y * y_scale) <= 0.25

instead?

2. Beam Waist from Source Amplitude in SLM

The beam waist used above is returned by self.source["amplitude_radius"] in the SLM class in hardware.slms.slm.py and is computed with SLM.fit_source_amplitude around line 921:

elif method == "moments":
    # Do moments in power-space, not amplitude.
    center = analysis.image_positions(np.square(amp))
    std = np.sqrt(2 * analysis.image_variances(np.square(amp), centers=center)[:2])

    center = np.squeeze(center)

However, for a Gaussian beam of the form:

I(x) ∝ exp(-2 * x² / ω₀²)

the variance of intensity is ω₀² / 4, so the beam waist ω₀ should be:

ω₀ = sqrt(4 * variance)

Thus, the current factor of 2 seems incorrect and may underestimate the beam waist.

Best,
Zhenpu

[ichristen]

Hi Zhenpu,

Thank you for your comprehensive report here!

Regarding 1:

mask = np.square(x * x_scale) + np.square(y * y_scale) <= 1

is correct. The intent of x_scale and y_scale is to define the edge of the unit disk and thus define the extent of canonical Zernike polynomials.

Currently, this edge is set (intentionally though somewhat arbitrarily) by slm.get_source_zernike_scaling() at the 1/e^2 amplitude radius (1/e^4 intensity radius), which is twice the "amplitude_radius" 1/e amplitude (1/e^2 intensity) beamradius parameter. This should probably be more documented.

I found that setting the mask to be smaller caused issues in practice, as there is a significant fraction of the power and amplitude outside the 1/e amplitude (1/e^2 intensity) beamradius. If the mask is used, this fraction of power is ignored. If the mask is not used, this fraction experiences r^O terms for high Zernike order O, which explode for r > 1.

If you find a contradictory reference for "Zernike-Gaussian beams" which defines things otherwise, we would consider changing the current behavior. Whatever the case, it should be documented more.

Regarding 2:

I'm fairly sure this is a bug. There shouldn't be the factor of 2 inside the sqrt. Whatever the case, "moments" and "fit" should return the same value for an analytic Gaussian beam that doesn't extend too far outside the SLM bounds: that's the test to verify that things are working. I'm testing this now.

Best,
Ian

[zzpwahaha]

Hi Ian,

Thank you for your prompt reply and for clarifying the definition of the 1/e⁴ intensity radius—it really helps.

I agree that rigorously defining a “unit disk” for an unbounded Gaussian beam is tricky. For my purposes, we’re simply using Zernike polynomials to correct the wavefront errors in our optics. As long as the correction converges effectively, and the polynomials remain largely orthogonal without necessitating prohibitively higher orders, the precise definition of the unit disk isn't critical for our purposes.

Thanks again for your insight!

Best regards,
Zhenpu

[ichristen]

Testing 2:

I tested 2 with the following code:

slm = SimulatedSLM((1600, 1200), pitch_um=(8,8))
slm.set_source_analytic(sim=True)

for method in ["fit", "moments"]:
    slm.fit_source_amplitude(method=method, force=True)
    print(method, slm.source["amplitude_radius"])

and got:

fit 2400.0
moments 2400.0

so I think things are working as intended. I think that the factor of 2 is to account for the fact that the moment calculation is done in power-space, not amplitude space. I don't precisely remember why I did that, but it's probably to suppress error from wavefront calibration in the edges of the SLM's field of view, fitting mainly the population you care about in the center.

However, I think there's still a bug here, as the standard "normal distribution" definition of a Gaussian used in gaussian2d

I(x) ∝ exp(- x² / 2 ω²)

Does not conform with the standard beamradius definition of a Gaussian beam expected by "amplitude_radius"

I(x) ∝ exp(-2 x² / ω²)

and I think this is currently incorrect (or at least documented incorrectly) in slmsuite.

[ichristen]

I'll give a bit of thought to how we want to define this. I think I found a few more areas where this is incorrect.