# Wavefront Aberration Wrought by Stochastic Surface Roughness and Mahajan’s Strehl Ratio Formula

How does surface roughness in optical components influence wavefront aberration? Surface roughness is wont to result in speckle, and I shall now show analytically as well as by posting some numerical simulation results that a good rule of thumb answer to this question is the following:

After propagating a significant distance (i.e. so that Fraunhofer diffraction relates the field on the two planes), the field looks like the ideal field attenuated by the square root of the Strehl ratio (which takes a particularly simple form for this stochastic roughness problem) together with a widely spread, roughly uniform intensity speckle pattern (the numerical results will help make the meaning of this clearer).

This result intuitively makes sense: especially when the stochastic perturbation field has spatial frequency content that is uniformly distributed in a wide spatial frequency range, whereas the main field is a spherical or plane wave with low spatial frequency variation on it. When these two are Fourier transformed in the Fraunhofer diffraction integral (or to calculate the field at the focal plane when dealing with a focussing field), the stochastic part sprays in all directions: it comprises plane waves mostly skewed at high angles relative to the main axial direction. Therefore the perturbation field swiftly sunders itself from the main focussing or paraxial beam, leading to a roughly uniform speckle pattern over a region much bigger than the beam waist.

## Analytical Demonstration of the High Spatial Frequency Rule of Thumb

So now we compute the focal plane or Fraunhofer-diffraction amplitude with standard Fourier optics assumptions holding. The response at the lateral position $(x, y)$ on the focal plane is approximately given by (here $f$ is the focal length):

$$\tag{1}\Phi(x,y) = \int_{\mathbb{R}^2} \exp\left(i\,\frac{2\,\pi\,(x\,X+y\,Y)}{\lambda\,f} + i\,\tilde{\phi}(X,Y)\right)\Psi_0(X,Y) \,\mathrm{d}X\,\mathrm{d}Y$$

where $\Psi_0(X,Y)$ is the “ideal” field, i.e. unmarred by surface roughness. It could be an aberration free field, but it could also include for the kind of low spatial frequency aberration $\exp(i\,\phi_{LF}(X,\,Y))$ that can be accurately resolved into Zernike functions i.e.

$$\tag{2}\phi_{LF}\left(X,\,Y\right)=\sum\limits_{\nu=-\infty}^\infty\sum\limits_{\ell = 0}^{|\nu|} \mathcal{Z}_\ell^\nu\left(\frac{r}{R}\right) \,e^{i\,\nu\,\theta}$$

is accurate when truncated to only a few terms and where $(r,\,\theta)$ are the polar co-ordinates corresponding to the Cartesian $(X,\,Y)$, $R$ is the farfield beam radius and $\mathcal{Z}_\ell^\nu$ are the Zernike polynomials.

Now, for a given constant $(x, y)$ we make the following coordinate rotation in the $(X, Y)$ plane:

$$\tag{3}U = \frac{x\,X+y\,Y}{\sqrt{x^2+y^2}};\quad V= \frac{-y\,X+x\,Y}{\sqrt{x^2+y^2}}$$

so that equation (1) becomes:

$$\tag{4}\Phi(x,y) = \int_{\mathbb{R}^2} \exp\left(i\,U\,\frac{2\,\pi\,\sqrt{x^2+y^2}}{\lambda\,f} + i\,\tilde{\phi}^\prime(U,V)\right)\Psi_0^\prime (U,V) \,\mathrm{d}U\,\mathrm{d}V$$

where the primed functions show that their argument coordinates have been rotated. But note that, given reasonable assumptions about the statistics of the surface roughness, the moments and other important statistics of the rotated coordinate aberration function $\tilde{\phi}^\prime(U,V)$ do not change. Now we imagine integrating equation (4) first with respect to $V$ along a line of constant $U$. In this integration the only variation comes from the term $\exp(i\tilde{\phi}^\prime(U,V))\,\Psi_0^\prime (U,V)$. We assume that $\Psi_0^\prime (U,V)$ varies very slowly compared to the roughness, so that we can break each slice into the integral with respect to $V$ up into sub-intervals of the form $[V_j,\,V_{j+1}]$ where $V_{j+1} – V_j$ is both:

1. Small enough that $\Psi_0^\prime (U,V)$ is roughly constant over this interval; but
2. Wide enough that if we can assume that $\tilde{\phi}^\prime(U,V))$ is a true random variable and make suitable ergodic assumption about the random process it belongs to so, then we can interpret the integral:$$\tag{5}\int\limits_{V_{j+1}}^{V_j} \exp(i\tilde{\phi}^\prime(U,V))\, \mathrm{d}\, V$$
as approximately the length of the line between $V_{j+1}-V_j$ times the expected value $\int_{-\infty}^\infty \exp(i\,\phi)\,p(\phi)\,\mathrm{d}\,\phi$ of $\exp(i\tilde{\phi}^\prime(U,V))$ given whatever probability density function $p(\phi)$ defines the statistical distribution of $\phi$.

Given these conditions, the integral for constant $U$ is:

$$\tag{6}\int\limits_{-\infty}^{\infty} \exp(i\,\phi)\,p(\phi)\,\mathrm{d}\,\phi\,\times\,\int\limits_{-\infty}^{\infty}\,\Psi_0^\prime (U,V)\,\mathrm{d}\,V$$

Thus we see the original equation (1) can be written:

$$\tag{7}\Phi(x,y) \approx \int_{\mathbb{R}^2} \exp\left(i\,\frac{2\,\pi\,(x\,X+y\,Y)}{\lambda\,f} \right)\Psi_0(X,Y) \,\mathrm{d}X\,\mathrm{d}Y \times \int_{-\infty}^\infty \exp(i\,\phi)\,p(\phi)\,\mathrm{d}\,\phi$$

In other words, this is simply the focal plane field distribution that would arise from a lens with no surface roughness times a constant attenuating factor:

$$\tag{8}\Gamma=\int_{-\infty}^\infty \exp(i\,\phi)\,p(\phi)\,\mathrm{d}\,\phi$$

set by the surface roughness’s probability distribution (i.e. it is the distribution’s characteristic function evaluated at unity).

If the surface roughness is such that the optical path error imprinted on the field belongs to a Gaussian distribution with zero mean phase and with RMS phase error $\sigma_w$ wavelengths (i.e. in radians the RMS value is $\sigma = 2\,\pi\,\sigma_w$) then:

$$\tag{9}\Gamma=\frac{1}{\sqrt{2\,\pi}\,\sigma}\,\int_{-\infty}^\infty \exp\left(i\,\phi – \frac{\phi^2}{2\,\sigma^2}\right)\,\mathrm{d}\,\phi = \exp\left(-\frac{\sigma^2}{2}\right)$$

and so we recover Mahajan’s Formula $\Gamma^2 = \exp(-(2\,\pi\,\sigma_w)^2)$ for the Strehl ratio. The original Mahajan reference here is:

V. Mahajan, “Strehl ratio for primary aberrations in terms of their aberration variance”, J. Opt. Soc. Am. 73 (6): 860–861

but Mahajan’s formula was a semi-empirical fit found to apply for low spatial frequency aberrations represented by Zernike polynomials. Here we see it has a theoretically “exact” interpretation for stochastic surface roughness. In this new interpretation, the Mahajan formula would not seem to be greatly sensitive to the surface roughness’s probability density function. For example, consider uniformly distributed surface roughness; equation (9) then becomes:

$$\tag{10}\Gamma=\frac{1}{2\,\sqrt{3}\,\sigma}\,\int_{-\sqrt{3}\,\sigma}^{\sqrt{3}\,\sigma }\exp(i\,\phi)\,\mathrm{d}\,\phi = \operatorname{sinc}(\sqrt{3}\,\sigma)$$

which, by plotting these functions, can be shown to be very like the Gaussian case for $\sigma$ up to about 0.2 waves RMS.

Optical power is conserved (i.e. the power output from the lens is the same as that going through the focal plane) and the diffraction integral transform is unitary, so the attenuation at first seems weird. The question naturally arises as to where the rest of the power has gone to. However, recall that Fourier optics used above is grounded on the assumption that the focal plane field is tightly confined around the focus. Therefore, the analysis above has simply found the field that is near the focus. The rest of the power, a fraction $1-|\Gamma|^2$ of the incident power, must therefore go through the focal plane well away from the focus. That is, this power is sprayed randomly and widely over the focal plane as speckle.

## Numerical Study of the High Spatial Frequency Case

With the numerical procedure I’m about to show, I have found that the Mahajan formula works extremely accurately for Gaussian surface roughnesses begetting RMS aberrations up to 0.5 waves RMS extremely well. This corresponds to a Strehl ratio of about $5\times10^{-5}$, or a central peak loss of more than 40dB. This is far beyond the extent to which Mahajan experimentally found his rule of thumb to fit behaviours for Strehl ratios greater than about 0.1 and he claims that his formula is good for any greater than 0.1 Strehl ratio. My own experience (by doing full Maxwell equation simulations of lens systems and comparing the results to those from the Mahajan formula) is that the Mahajan formula is only reliable in the low spatial frequency case down to about 0.3 Strehl and is much less accurate than when it is applied to the present, high spatial frequency case in the way I do. My drawing below shows the numerically simulated system.

Figure 1: Schematic Simulated System

Here the full vector Maxwell’s equations are integrated numerically by split step Fourier-leapfrog algorithm of my own devising and whose accuracy is limited only by the transverse grid spacing and the axial step size. In both cases, the system input was a plane wave of 488nm freespace wavelength: first a linearly polarised wave, then a circular polarised one.

The surface roughness model was a Gaussian random process with exponential autocorrelation function:

$$\tag{11}\mathcal{R}(\mathbf{x}_1,\,\mathbf{x}_2) = \sigma^2 \exp\left(-\frac{|\mathbf{x}_1-\mathbf{x}_2|}{L}\right)$$

with a correlation length $L$ equal to $1.5\mu\mathrm{m}$. The simulation region’s sideways breadth was $100\mu\mathrm{m}\times100\mu\mathrm{m}$ and an example surface roughness topographic map with the above statistics is shown below where the RMS surface roughness is 150nm.

Figure 2: Sample Surface Roughness with $1.5{\rm \mu\,m}$ Correlation Length and $150{\rm nm}$ RMS roughness. Field of View is $100{\rm \mu\,m}\times 100{\rm \mu\,m}$.

Some results of the numerical simulations are shown below, and the central peak, with the same shape as the theoretical ideal (i.e. zero surface roughness) case, can clearly be seen, as can be the widely spread speckle power. In the linearly polarised case, electric field vector sets a preferred direction for the speckle, and so we see that it is not isotropic, whereas, for the circularly polarised case, there is no such preferred direction and the speckle is clearly isotropic.