1. Introduction

Hermite Gauss (HG) modes are well known solutions of the paraxial wave equation and are widely used in the laser community [1], from optical tweezers using gaussian beams (for which the 2018 Nobel prize was awarded), to interferometry and multi-mode fibres. HG modes are ideal for representing beams within spherical optical systems because their transformation through these systems as well as their propagation through free space can be described with simple ABCD matrices from ray optics. Even complex optical systems such as the LIGO and LISA interferometers can be approximately analyzed using HG modes and ABCD matrices [2]. Imperfections in optical systems are then often described by scatter matrices which transfers energy between modes. Once these scatter matrices are known, optical systems can be simulated very efficiently using standard matrix multiplication techniques [3]. Linearized versions of these scatter matrices are very well known for small imperfections such as small angular and transversal misalignments or offsets in waist locations and waist sizes (mode mismatch) [2,4].

A disadvantage of HG modes can be that numerical calculations of scatter matrices for general beam distortions are computationally expensive [5,6], unless analytical expressions, such as those we provide here, exist. While HG modes may not converge very well for light distributions containing high spatial frequencies, these higher spatial frequency components will rapidly diffract away, making this approach suitable for many interesting studies.

One case which will become more and more relevant as laser fields enter many different areas of science and technology are phase-distorted apertured HG modes. Indeed, many optical components such as mirrors, pupils, or even photodiodes themselves, can be described using this formalism—a transformation and circular clipping—if the propagated beam description is in terms of HG modes [7]. These distortions are well described by Zernike polynomials, an orthonormal set of eigenmodes defined on a circular aperture. These polynomials are widely used in astronomy to characterize wavefront errors in optical systems [8]. In this paper we provide analytical solutions of the scatter matrices for an arbitrary set of HG eigenmodes which pass through an aperture and experience wavefront distortions described by Zernike polynomials. An example of such a system is the LISA optical system in which a Gaussian mode passes through a beam expanding telescope which also apertures the outgoing beam to maximize the power at a distant receiving telescope of given diameter. The beam then is compressed by a reverse process to mix it with the local oscillator; another (nearly) Gaussian beam [9].

This paper is organized as follows: In the next chapter we define the HG modes and the Zernike polynomials as well as the parameters we use to describe them. In it we also outline the decomposition scheme used to model light, and the analytic expression for the projections between HG modes of different bases through an aperture. The third chapter applies this to the propagation of clipped Gaussian beams, comparing the results to those found in the literature obtained through classical scalar diffraction theory. The fourth chapter adds wavefront distortions in the form of Zernike polynomials to the apertured Gaussian beam, examining how these affect the transmitted field. Chapter five talks about a specific applications of the modal decomposition to the LISA mission, while the final chapter, six, summarizes our results and describes future directions for this research. We finish with the appendix in which we mathematically derive the analytic expression for the overlaps.

2. Hermite Gauss modes, Zernike polynomials, and diffraction at circular apertures

HG modes are solutions of the paraxial wave equation, so that their superposition can well represent light distributions anywhere the light is composed of plane wave vectors restricted to being nearly directed along a single optical axis. Generally, these are regions where all wave vectors make angles with the optical axis $\theta$ small enough so that $\sin \theta {\, \longrightarrow \,} \theta$ [1].

Once the optical axis is defined, the HG modes are uniquely described by the light’s wavelength and two parameters, the waist radius $w_0$ of the beam, and the location along the optical axis of said waist, $z_0$. We write out the modes in terms of physically meaningful quantities; the radius of the beam $w(z)$, the beam wavefront radius of curvature $R(z)$, and the Gouy Phase $\psi (z)$ which is an additional phase per mode that builds because of its deviation from a plane wave [1]:

These HG modes, for any choice of $w_0$ and $z_0$, represent a complete orthonormal set of basis functions for paraxial light. Propagation of a superposition of such modes along the optical axis is a simple matter of adjusting the width, curvature, and Gouy phase for all modes, changing the $z-z_0$ to $z+\Delta z - z_0$ with $\Delta z$ the propagation distance.

Wavefront distortions caused by an optical field are typically described by a phase or wavefront error (WFE) map $\Phi (x,y)$ [7]. In circular apertures these distortions can be expressed as a linear combination of the orthogonal set of Zernike polynomials:

2.1 Prior solutions to the propagation of truncated paraxial modes

The scattering of paraxial light modes through circular apertures has been extensively studied. In 1971 R.G. Schell and G. Tyras derived an analytic description of the irradiance pattern of the truncated Gaussian Beam, valid in regimes similar to the Fresnel region for the scattering [10]. A similar representation was later derived and examined by C. Campbell in 1987 [11]. Propagations of higher order Laguerre Gauss beams (another basis for paraxial light distinct from Hermite Gauss modes) through circular apertures in the Fresnel region was expressed again as a sum of Bessel functions by G. Lenz in 1996 [12]. However these expansions as infinite sums of Bessel functions are difficult to propagate through optical systems, and in general it’s complicated to examine what happens when there’s further diffraction of the light. For the cases of multiple apertures these diffraction methods can only be applied to diffraction at the initial aperture.

An alternative is to find the diffraction beam decomposition in terms of paraxial modes, since then the same technique can be used at each following aperture. This was one of the approaches taken in a 1972 paper, in which the authors looked at the overlap between outgoing Laguerre Gauss modes and a general truncated Laguerre Gauss mode that had been propagated to the Fresnel region [13]. By rewriting the propagated beam in the same form as the entering beam the same propagation technique could be repeated multiple times for multiple apertures, something more desirable for application.

Each of these methods begin with a diffraction integral having Kirchoff boundary conditions at the aperture, yielding Kirchoff’s diffraction formula:

While we utilize a method similar to some of the prior work, decomposing the propagated field into the same types of modes as the incoming field to allow for multiple iterations at multiple apertures, we avoid propagating the beam first. Instead we decompose the field in the plane of the aperture into the HG mode components, and then propagate the modes themselves. Besides being able to use this with multiple apertures, another advantage is that we no longer require different sets of solutions for different regions of space. The integral to be solved is also much simpler than the diffraction integral in Eq. (8). Finally, with this method we can also easily incorporate the effects of small wavefront distortions (Section 4).

2.2 HG-scatter matrix at the circular aperture

In order to decompose the light field at the aperture into HG modes, we assume at the very least we already had a description of the non-apertured beam in terms of these modes, potentially up to some additional wavefront error. For wavefront distortions described by some general phase map $\Phi (x,y)$ present in a given clipped light field $E(x,y)$, the amount of field exiting said aperture in HG mode $(m,\,n)$ is given by [15,16]

 

Fig. 1. Illustration of beam decomposition and reconstruction procedure. Graphic created by Simon Barke for this paper.

Download Full Size | PDF

2.3 Applications to propagation

As has been discussed in the prior section, once we have our initial field $E(x,y)$ at some plane decomposed into a sum of HG modes from some basis, propagating through space is as simple as adding the propagation distance to the $q$-parameter. The changes to width, curvature, and even Gouy phase can all be derived from this change. Calculating the complex amplitude coefficients for transmission through the next aperture or mirror is equivalent to solving expressions of the form:

3. Propagation of clipped Gaussian beams

We aimed to study the effectiveness of this method by using it to propagate a clipped Gaussian beam to compare with published solutions. We compare to numerical computations of the scattering integral of Eq. (8), so that results should coincide with each of the studies mentioned in section 2.1. Specifically, the larger number of clipping parameter examples provided for the clipped Gaussian beam in the 1971 Schell and Tyras study [10] make it ideal for comparison and testing of our method. In this paper Gaussian beams having the form

In Fig. 2 we compare our field distributions, truncated at varying orders in HG mode, with the field obtained via numerical approximations to Eq. (8). The comparison to the actual distribution in the clipping plane is shown at the top left of Fig. 2. The top right of Fig. 2 compares the on-axis intensity of the results from numerical integration (solid lines) of Eq. (8), utilizing similar approximations to those used in [10], to our expansions (dashed lines) in an optimal basis for representing the light through orders 0, 10, 20, and 30. There are several important features to note, the most apparent of which is that the more drastically clipped beams will require more modes for a better fit. This makes sense because the sharper the clipping the higher the spatial frequency components necessary to represent it. The oscillatory behavior of the more strongly clipped Gaussian beams along the optical axis is characteristic of beams passing through circular apertures.

 

Fig. 2. Comparisons between our expansion of the light to various maximum mode order (defined as maximal $m+n$) and the distributions as obtained through the literature. Approximations to maximum mode order 0, 10, 20, and 30 are shown as dashed lines increasing in brightness as we increase the number of modes used to approximate the light.

Download Full Size | PDF

The bottom graphs of Fig. 2 shows comparisons between the same approximations within the planes $159a$ and $2000a$ away from the transmitting aperture. These show clearly that our approximations increasingly match a simple Gaussian approximation (using only the lowest order mode) the further we get from the aperture. Higher spatial frequency Fourier components of beams diverge/spread faster than the lower spatial frequency components, so near the optical axis we are left with the latter which are described well by lower order HG modes. At $2000a$ from the clipping aperture the lowest order mode is nearly sufficient to on its own represent the light. As a final note, as we move further from the optical axis, it appears that all approximations (and published results) approach the lowest order term, likely a vestige of the Gaussian envelope all the modes share.

We have included the coefficients for the optimal fit to represent the light up through order 30 (first 496 modes) in Table 2 of appendix C, and depicted its residual from the true distribution in the plane of the aperture (Fig. 3) as well as the predicted central 10 mrad of the far field (Fig. 4). Fitting with orders $\le 30$ was good enough to represent $99.73\%$ of the power in the $F=.5$ beam, $99.98\%$ for $F=1$, and $>99.99\%$ for $F=2$. Best fit basis parameters used for the representational modes were $w_0 = .251a$ for both $F=.5$ and $F=1$, while $w_0 = .25a$ for $F=2$, and in all three cases the wavefront was given no curvature at the aperture i.e. $z_0 = 0$. Note that as the bases are very similar, the fit residuals will be maximized at the same places, since this is representative of where the bases lacks the ability to well represent light. The most pronounced residual is found at the aperture edges, and with the stronger clipping we get the greatest residual. We see this quantitatively in the mode coefficients themselves. The coefficients drop off more sharply with order for $F=2$ which is more weakly clipped, reaching as low as $10^{-5}$ in magnitude for coefficients to modes with order as low as 6. Contrast this to the $F=.5$ case for which we must go to order 16 before the mode coefficients drop as significantly. Note that Table 2 coefficients will work for any clipped Gaussian problem with the same clipping number $F$, as long as we scale the best fit basis appropriately with aperture size.

 

Fig. 3. Log Absolute Residual of fit at the aperture. The min and max residuals all occur at the same place because the residuals are basis dependent and each of the representing basis were nearly identical.

Download Full Size | PDF

 

Fig. 4. Log of the Far Field Intensity Relative to the Far Field Intensity at the Optical Axis from both HG expansion (solid) and Fresnel Approximation (dashed) for the 3 clipped Gaussian beams. Both the paraxial expansion and Fresnel approximation break down the further we get from the optical axis, but we see these two match well for the first several diffraction rings.

Download Full Size | PDF

Looking at the central 10 mrad of the far field, we see the effects of the clippings in numerous rings of maxima and minima, something that would not be present in a purely Gaussian beam.

3.1 Strongly clipped beams

It is clear that the stronger the clipping, the worse our fit with the same number of modes. To probe the limits of HG modes’ representational capabilities we look at a more extreme case of clipping, comparable to the propagation scenarios examined by Huang and Ding [17]. They started with a diffraction integral similar to prior work, rewriting everything in terms of dimensionless variables, and reduced the problem as others had to a 1-D integral over a Bessel function times a Gaussian. They solved the integral by rewriting the Bessel function portion of the integrand as a sum of Gaussian functions, allowing for easier integration. They then examined an extreme case, where the parameter $B \equiv \left (\frac {a}{w_0}\right )^{2}$ was equal to 1. Contrasting this to the earlier work by Schell and Tyras, who went as low as $F = .5$, $B=1$ corresponds to $F \approx .32$. We test the validity of our method’s predicted field in the Fresnel region against numerical computation of the integral Huang and Ding had worked from, the same way they had tested their model, seen in Fig. 5.

 

Fig. 5. Comparison of Numerical Solutions to the Diffraction Integral obtained by Huang and Ding in [17] (black line) with our HG expansions to various mode order (dashed lines).

Download Full Size | PDF

We have compared over the same field regions as [17] and used the same convention that $z_0\equiv \frac {a^{2}}{\lambda }$ is now the Fresnel Distance of the aperture. We see a noticeable point-wise difference from our fit in the extreme near field close to the optical axis. However the overall distribution in the plane largely matches our approximation to mode order 30.

Due to extreme clipping, our expansion to order 30 only captures 99.4% of the light’s total power, compared to the best fit order 20 and 10 mode expansions, holding 99.14% and 98.48% respectively. This illustrates that while there can be point-wise differences, we have correctly accounted for the majority of the power; higher-than-30 order modes can only hold, and therefore contribute, an additional $.6\%$ of the total power. Appendix D. shows how the fractional RMS deviation from the true distribution, a more intuitive metric for the deviation in the distributions shape, can be solely determined from the captured power. For the order 30 fit the fractional RMS amplitude deviation is $.0775$. More noticeable is that we can capture the majority of the power in the lower order modes, having order $\le 10$, leaving out only the smaller high spatial frequency contribution from the clipping. As we move further from the source, even our lower order approximations accurately represent beam in the axial region at which a detector might be placed. In the language of Fourier optics, the higher spatial frequency terms that couple to higher order HG modes will have scattered towards infinity.

4. Adding in phase distortions

In this section we add wavefront distortion to the most strongly clipped beam. As mentioned in section 2, these wavefront error (WFE) maps are well represented by linear combinations of Zernike polynomials. Such an additional spatially varying phase distortion over the aperture (Fig. 6) has been added to the most strongly clipped Gauss beam, built from Zernike polynomials with coefficients chosen arbitrarily (Table 1).

Table 1. Zernike Coefficients of Distortion

View Table | View all tables in this article

This alters the complex coupling coefficients between modes as described before (See Eq. (12)); the coefficient of mode $(m,\,n)$ is now a sum, or using the same notation as introduced before

Figures 7–10 depict graphically the impact of the now complex HG-mode coefficients on the far field. Figures 7 and 8 show the intensity of the distribution in the far field over the central 2 mrad region as well as the residual of this from the distribution with no wavefront distortion present. Note that the main features in Fig. 7 are the diffraction rings shared by the clipped gaussian without the additional wavefront error. However, the higher residuals about these ring locations in Fig. 8 from the case with no additional wavefront shows that our additional phase has led to some slight displacement and azimuthal asymmetry in these rings.

 

Fig. 6. Additional Phase in Radians

Download Full Size | PDF

 

Fig. 7. Log far field Intensity in 2 mrad central cone.

Download Full Size | PDF

 

Fig. 8. Log Absolute Intensity Residual in central 2 mrad of far field due to phase map.

Download Full Size | PDF

Figures 9 and 10 show the phase and phase residuals instead. We’ve taken the logarithm of the residuals to make the differences more easily visible. While small differences are clearly visible, the residual is not nearly as pronounced as the phase distortion itself. Figure 9 shows that the phase remains largely constant between the diffraction rings, but changes as expected at each ring. However, the added wavefront distortion deforms the rings which leads to rather large differences between the original and the distorted field at those locations. This is shown in Fig. 10. For many applications including LISA, we are mostly interested in the phase changes and the resulting gradient in the central region of the far field which are significantly smaller.

 

Fig. 9. Phase in 2 mrad central cone of far field.

Download Full Size | PDF

 

Fig. 10. Log of Absolute Phase Difference due additional phase map in central 2 mrad cone of far field.

Download Full Size | PDF

These figures show the power of the method for a known but random wavefront error of a given telescope. Testing multiple random linear combinations of distortions can then be used to derive generic requirements for the wavefront quality of the optical system.

4.1 Connecting the far field to the WFE

A potentially better way is to calculate the contributions of individual Zernike distortions to for example the phase gradient in the far field to guide the polishing process. While a detailed study of this is well beyond the scope of this paper, we present two specific individual Zernike distortions, defocus and oblique astigmatism shown in Fig. 11, as initial examples.

 

Fig. 11. The Additional WFE maps being added.

Download Full Size | PDF

The defocus only leads to radial changes in the distribution. Indeed, looking at the effects of adding a more and more pronounced defocus, Fig. 12, to the transmitted beam smears the diffraction rings due to the clipping, eventually getting rid of the first ring altogether (along with the pronounced phase change that would have occurred at said ring).

What we expect from the astigmatism is for the development of azimuthal asymmetries, maintaining $180^{\circ }$ rotational symmetry in field and $90^{\circ }$ in intensity. Figure 13 shows the loss of the primary diffraction minima ring at points between a set $90^{\circ }$ interval, and in fact we see a smearing of the phase with features similar to the original oblique astigmatism itself.

 

Fig. 12. The changing [log of] intensity (top) and phase (bottom) of the distributions within a 2 mrad cone of the far field in response to an increasing defocus in the phase of the transmit beam (going from a maximal phase offset over the aperture of 1 mrad to 1 radian by factors of 10 as we move from left to right). Note the smoothing of the entire phase as we increase the aberration.

Download Full Size | PDF

 

Fig. 13. The changing log of intensity (top) and phase (bottom) of the distributions within a 2 mrad cone of the far field in response to an increasing oblique astigmatism in the phase of the transmitted beam (going from a maximal phase offset over the aperture of 1 mrad to 1 radian by factors of 10 as we move from left to right). Note the slower phase transition along a line at $45^{\circ }$ to the axis.

Download Full Size | PDF

5. Real world applications: LISA

Mentioned briefly before, the Laser Interferometer Space Antenna (LISA) is a future space based, laser-interferometric, gravitational-wave observatory aiming to measure much lower frequency gravitational waves than its LIGO counterpart. The observatory is composed of three spacecrafts (SC) in a triangular configuration, separated by $L = 2.5\,\rm {Gm}$ and exchanging lasers with each other. Thus, each SC pair acts as an arm of an extremely long interferometer.

Two $\sim 30\,$cm telescopes on each spacecraft are used to transmit $2\,\rm {W}$ of near infrared light ($\lambda = 1064\,$nm) towards and to receive around $700\,\rm {pW}$ from each of the other two spacecrafts. The main science signals will be beat signals between the received and a local laser beam. The phase evolutions of these beat signals are then used to extract the length changes due to gravitational waves using time delay interferometry. [9]

Apparent length changes can be caused by angular spacecraft motion $\delta \phi$ of the sending spacecraft which scans the large cone of the far field across the receiving telescope. A perfectly spherical beam in which the radius of curvature is identical to the distance would be immune to such spacecraft jitter. However, any phase front gradient $\nabla _t \Phi _0$ allows spacecraft jitter $\delta \phi$ to couple into the length measurement,

The spacecraft jitter in LISA is expected to be $\delta \phi < 10\,\rm {nrad}/\sqrt {\textrm {Hz}}$ which causes the beam at the far spacecraft ($L=2.5\,\rm {Gm}$) to laterally move randomly by $25\,\rm {m}/\sqrt {\textrm {Hz}}$. This can be used to derive the requirements for the phase gradient of the far field:

However, there will be additional phase distortions present in the transmit beam in addition to any caused by the clipping. Using the formalism presented here, we can calculate the phase gradient for any such WFE that can be described as a linear combination of Zernike polynomials. The top part of Fig. 14 shows a sample phase map constructed from such a linear combination with amplitudes falling off with the mode index and some random variations as shown in the lower part of Fig. 14. We then use our modal decomposition to derive the far field at the receiving spacecraft and calculate the phase gradient for a constant distance.

 

Fig. 14. Sample Phase map normalized to maximum amplitude of 1 with Zernike amplitudes declining as shown.

Download Full Size | PDF

In the case of LISA, the received field is expected to be actively centered within a few $10\,\rm {nrad}$ on the telescope [9] such that only the central part of the received wavefront is critical.

With this in mind, we can now take a sample WFE composed of Zernike polynomials, calculate the phase gradient in the critical central region in the far field and scale then the amplitudes until the transversal phase gradient is sufficiently small. Breaking down this process for the above phase map, we first calculate the mode coefficients for numerous scalings of the WFE. This gives us a description of the field in terms of the modes from which we can calculate the phase, and from there the phase gradient, of the far field. The far field phase gradients are depicted for five WFE scalings in Fig. 15 below.

 

Fig. 15. Logarithm of the Phase Gradients in the 200 nrad by 200 nrad nominal center of the far field. Note that the non-symmetric Zernike terms shift the beam’s center.

Download Full Size | PDF

Without distortion ($0\lambda$), the far field is nicely centered with virtually no gradient. For a telescope WFE of $\lambda /40$, the gradient in the 200 nrad by 200 nrad window mostly meets the LISA requirements. For larger WFE, the beams will have to be re-centered at the receiving spacecraft to stay within the blue regions of each figure. A second risk reduction strategy is to use the measured beam jitter, calibrate its jitter to length coupling and subtract it from the main data. However, a detailed discussion about mitigation techniques to minimize the jitter to length coupling is not the subject of this paper. But having a versatile and fast tool to calculate phase gradients in the far field as a function of telescope WFEs typically described by Zernike polynomials will help to study this problem in great detail in the future.

6. Conclusion

What has been shown so far is a small sample of the modeling power of analytic HG-mode expansion methods. In the future, we will study the allowed WFE for LISA in more detail, studying not only a more diverse set of WFEs but also adding for example the wavefront sensing and active alignment system to the analysis and verify the proposed jitter to length subtraction.

Appendix A. General expression

Utilizing the same notation as in section 2.3, each of $a,\,b,\,c,\,d,\,f,\,$ and $g$ are non-negative integers. As in that section, $a$ and $b$ label the outgoing, or representing HG mode. $f$ and $g$ represent the HG mode incident on or entering the aperture. $\pm c$ and $d$ represent the upper and lower indices respectively of the present Zernike polynomial, so as mentioned before should have $d-c$ even. Finally, $r$ is the radius of the aperture, $o$ is utilized to represent outgoing or representative beam parameters, and $i$ to represent the incoming parameters.

Although in section 2.3 the entire solution was presented in a single formula, for the sake of proving it we break into two cases. The first, with positive upper Zernike index, has that an overlap of zero when either $ {b+g \; {\textrm{or}} \; a+f+c \; {\textrm{are odd}}}$, and is otherwise:

(23)$$\begin{aligned}&{}_{w_o,z_o,\lambda_o} \langle {a,b}\vert {c,d} \vert {f,g}\rangle _{w_i,z_i,\lambda_i}(r,z) = \sigma \sum_{A=0}^{\left \lfloor\frac{a}{2}\right \rfloor}\sum_{B=0}^{\left \lfloor\frac{b}{2}\right \rfloor}\sum_{C=0}^{\left \lfloor\frac{d-c}{2}\right \rfloor}\sum_{D=0}^{\left \lfloor\frac{c}{2}\right \rfloor}\sum_{F=0}^{\left \lfloor\frac{f}{2}\right \rfloor}\sum_{G=0}^{\left \lfloor\frac{g}{2}\right \rfloor}\left[\frac{\left(\frac{a+b+d+f+g}{2}-(A+B+C+F+G)\right)!}{\left(\frac{a+b+c+f+g}{2}-(A+B+F+G)\right)!}\right. \\ &\times\frac{(-1)^{D}\left(-\gamma\right)^{A+B+C+F+G}r^{2C}W_o^{2(A+B)}W_i^{2(F+G)}(d-C)!\Gamma\left(D-(B+G)+\frac{b+g+1}{2}\right)\Gamma\left(\frac{a+c+f+1}{2}-(A+D+F)\right)}{8^{A+B+F+G}A!B!C!F!G!(a-2A)!(b-2B)!\left(\frac{d+c}{2}-C\right)!\left(\frac{d-c}{2}-C\right)!\left(f-2F\right)!\left(g-2G\right)!(2D)!\left(c-2D\right)!} \\ &\left.\times \left(1-{\textrm{e}}^{-\gamma r^{2}}\sum_{K=0}^{\frac{a+b+d+f+g}{2}-(A+B+C+F+G)}\frac{\left(\gamma r^{2}\right)^{K}}{K!}\right)\right] \end{aligned}$$

In the second case, for Zernike terms with negative superscripts, we have that the overlap is zero when b + g or a + f + c are even, otherwise:

(24)$$\begin{aligned}&{}_{w_o,z_o,\lambda_o}\langle {a,b}\vert {-c,d} \vert {f,g}\rangle _{w_i,z_i,\lambda_i}(r,z) = \sigma \sum_{A=0}^{\left \lfloor\frac{a}{2}\right \rfloor}\sum_{B=0}^{\left \lfloor\frac{b}{2}\right \rfloor}\sum_{C=0}^{\left \lfloor\frac{d-c}{2}\right \rfloor}\sum_{D=0}^{\left \lfloor\frac{c-1}{2}\right \rfloor}\sum_{F=0}^{\left \lfloor\frac{f}{2}\right \rfloor}\sum_{G=0}^{\left \lfloor\frac{g}{2}\right \rfloor}\left[\frac{\left(\frac{a+b+d+f+g}{2}-(A+B+C+F+G)\right)!}{\left(\frac{a+b+c+f+g}{2}-(A+B+F+G)\right)!}\right. \\ &\times\frac{(-1)^{D}\left(-\gamma\right)^{A+B+C+F+G}r^{2C}W_o^{2(A+B)}W_i^{2(F+G)}(d-C)!\Gamma\left(D-(B+G)+\frac{b+g+2}{2}\right)\Gamma\left(\frac{a+c+f}{2}-(A+D+F)\right)}{8^{A+B+F+G}A!B!C!F!G!(a-2A)!(b-2B)!\left(\frac{d+c}{2}-C\right)!\left(\frac{d-c}{2}-C\right)!\left(f-2F\right)!\left(g-2G\right)!(2D+1)!\left(c-2D-1\right)!} \\ &\left.\times\left(1-{\textrm{e}}^{-\gamma r^{2}}\sum_{K=0}^{\frac{a+b+d+f+g}{2}-(A+B+C+F+G)}\frac{\left(\gamma r^{2}\right)^{K}}{K!}\right)\right] \end{aligned}$$

In all of these we have defined the following variables:

(25)$$\begin{aligned}\delta z_{i/o} \equiv z-z_{i/o}\qquad z_{ri/ro} \equiv \frac{\pi w_{i/o}^{2}}{\lambda_{i/o}}\qquad W_{i/o} \equiv w_{i/o}\sqrt{1+\left(\frac{\delta z_{i/o}}{z_{ri/ro}}\right)^{2}}\qquad \psi_{i/o} \equiv \arctan\left(\frac{\delta z_{i/o}}{z_{ri/ro}}\right) \\ \gamma = \frac{1}{w_o^{2}\left(1+i\frac{\delta z_o}{z_{ro}}\right)}+ \frac{1}{w_i^{2}\left(1-i\frac{\delta z_i}{z_{ri}}\right)}\qquad \sigma = \frac{c!\sqrt{a!b!f!g!}2^{a+b+f+g+1}{\textrm{e}}^{i\left(k_o\delta z_o-k_i\delta z_i+\left(f+g+1\right)\psi_i-\left(a+b+1\right)\psi_o\right)}}{\pi r^{d}W_o^{a+b+1}W_i^{f+g+1}\sqrt{\gamma^{a+b+d+f+g+2}}} \end{aligned}$$

where $z$ is the aperture position on the optical axis.

We mention that for cases where the radius of the aperture is far smaller than the width of both incoming and outgoing basis at the aperture, we have $\left \lvert {\gamma r^{2}} \right \rvert \ll 1$ so that it may be preferable to simplify:

(26)$$\begin{aligned}1-{\textrm{e}}^{-\gamma r^{2}}\sum_{K=0}^{M}\frac{\left(\gamma r^{2}\right)^{K}}{K!} &= 1-{\textrm{e}}^{-\gamma r^{2}}\left({\textrm{e}}^{\gamma r^{2}}-\sum_{K=M+1}^{\infty}\frac{\left(\gamma r^{2}\right)^{K}}{K!} \right) \\ &= {\textrm{e}}^{-\gamma r^{2}}\sum_{K=M+1}^{\infty}\frac{\left(\gamma r^{2}\right)^{K}}{K!}{\, \longrightarrow \,} {\textrm{e}}^{-\gamma r^{2}}\frac{\left(\gamma r^{2}\right)^{M+1}}{(M+1)!} \end{aligned}$$

While this may not appear simple or elegant, it is exact and can be implemented into programs providing the overlap between incoming and outgoing modes through an aperture, without the need for numerical integration. There is reason to believe that this method could therefore give a significant computational speed advantage over other beam propagation methods, but a direct comparison is outside the scope of this work.

Appendix B. Derivation

In the following we use the same shorthand or encapsulated variables as from Appendix A. We first evaluate the positive Zernike term case. We start by expanding out each of the Hermite and Zernike polynomials

(27)$$\begin{aligned}{}_{o}\langle {a,b}\vert {c,d} \vert {f,g}\rangle _{i} &= \int_0^{r} \rho \textrm{d} \rho\int_0^{2\pi} \textrm{d} \phi \tilde{U}_{a,b}^{*}\left(\rho,\phi,z;w_o,z_o,\lambda_o\right)Z^{c}_d\left(\frac{\rho}{r},\phi\right)\tilde{U}_{f,g}\left(\rho,\phi,z;w_i,z_i,\lambda_i\right) \\ &=\frac{\sqrt{a!b!f!g!}{\textrm{e}}^{i\left((f+g+1)\psi_i+k_o\delta z_o-(a+b+1)\psi_o-k_i\delta z_i\right)}}{\pi W_o^{a+b+1} W_i^{f+g+1}r^{d}\sqrt{2^{a+b+f+g-2}}}\sum_{A=0}^{\left\lfloor\frac{a}{2}\right\rfloor}\sum_{B=0}^{\left\lfloor\frac{b}{2}\right\rfloor}\sum_{C=0}^{\left\lfloor\frac{d-c}{2}\right\rfloor}\sum_{F=0}^{\left\lfloor\frac{f}{2}\right\rfloor}\sum_{G=0}^{\left\lfloor\frac{g}{2}\right\rfloor} \\ & \left[\frac{(-1)^{A+B+C+F+G}\left(d-C\right)!2^{\frac{3}{2}\left(a+b+f+g-2(A+B+F+G)\right)}W_i^{2(F+G)}W_o^{2(A+B)}r^{2C}}{A!B!C!F!G!(a-2A)!(b-2B)!\left(\frac{d+c}{2}-C\right)!\left(\frac{d-c}{2}-C\right)!(f-2F)!(g-2G)!}\right. \\ & \times \int_0^{r} \textrm{d} \rho \rho^{a+b+d+f+g-2(A+B+C+F+G)+1}{\textrm{e}}^{-\gamma\rho^{2}} \\ & \left.\times\int_0^{2\pi} \textrm{d} \phi \cos(c\phi)\left(\sin\phi\right)^{b+g-2(B+G)}\left(\cos\phi\right)^{a+f-2(A+F)}\right] \end{aligned}$$

where we replace the $\cos (c\phi )$ with $\sin (c\phi )$ for negative upper Zernike index. We have two separate types of integrals, but we begin with integrals of the form

(28)$$G(a,b) \equiv \int_0^{2\pi} \textrm{d} \phi \sin^{a}\phi \cos^{b}\phi$$

because both $\cos (c\phi )$ and $\sin (c\phi )$ can be written as a sum of powers of $\sin$ and $\cos$.

B.1. Integral 1 solution

We skip how we got to the solutions form and simply prove it to be correct, we find with the same nomenclature

(29)$$\begin{aligned} G(a,b) & = 2\frac{\Gamma\left(\frac{a+1}{2}\right)\Gamma\left(\frac{b+1}{2}\right)}{\Gamma\left(\frac{a+b}{2}+1\right)}\qquad & \textrm{when}\; a\; \textrm{and} \; b \; \textrm{are even} \\ & = 0\qquad & \textrm{when either}\; a \;\textrm{or}\; \textrm{b}\; is\; odd \end{aligned}$$

To prove these cases we start off by showing $G(a,b) = G(b,a)$,

(30)$$G(a,b) = \int_0^{2\pi} \textrm{d} \theta \sin^{a}\theta \cos^{b}\theta$$

making the variable change $\phi = \frac {\pi }{2}-\theta$ gives us

(31)$$\begin{aligned} =-\int_{\frac{\pi}{2}}^{-\frac{3\pi}{2}} \textrm{d} \phi\sin^{a}\left(\frac{\pi}{2}-\phi\right)\cos^{b}\left(\frac{\pi}{2}-\phi\right)=\int_{-\frac{3\pi}{2}}^{\frac{\pi}{2}} \textrm{d} \phi\cos^{a}\phi\sin^{b}\phi \\ = \int_{-\frac{3\pi}{2}}^{0} \textrm{d} \phi\cos^{a}\phi\sin^{b}\phi + \int_0^{\frac{\pi}{2}} \textrm{d} \phi\cos^{a}\phi\sin^{b}\phi \end{aligned}$$

changing variables in the first integral to $\gamma = \phi + 2\pi$ and in the second integral $\gamma = \phi$ gives us

(32)$$\begin{aligned} =\int_{\frac{\pi}{2}}^{2\pi} \textrm{d} \gamma \cos^{a} (\gamma-2\pi) \sin^{b} (\gamma-2\pi) + \int_0^{\frac{\pi}{2}} \textrm{d} \gamma \cos^{a}\gamma \sin^{b}\gamma \\ =\int_0^{2\pi} \textrm{d} \gamma \sin^{b}\gamma \cos^{a} \gamma = G(b,a) \end{aligned}$$

concluding the first part of our proof, $G(a,b) = G(b,a)$. The next part is to show it must be zero for odd values of $a$ or $b$. We start off again

(33)$$G(a,b) = \int_0^{2\pi} \textrm{d} \theta \sin^{a} \theta \cos^{b}\theta$$

and make the change $\phi = -\theta$ which gives

(34)$$=-\int_0^{-2\pi} \textrm{d} \phi \sin^{a}(-\phi)\cos^{b}(-\phi) = \int_{2\pi}^{0} \textrm{d} \phi (-1)^{a}\sin^{a}\phi \cos^{b}\phi\,.$$

By making the final change to $\gamma = \phi +2\pi$ gives

(35)$$\begin{aligned}=(-1)^{a}\int_0^{2\pi} \textrm{d} \gamma \sin^{a}(\gamma-2\pi)\cos^{b}(\gamma-2\pi) \\= (-1)^{a}\int_0^{2\pi} \textrm{d} \gamma \sin^{a}\gamma\cos^{b}\gamma = (-1)^{a} G(a,b) \\ {\, \Longrightarrow \,} G(a,b) = (-1)^{a}G(a,b) \end{aligned}$$

telling us that we must have $G(a,b) = 0$ when $a$ is odd, and by the last step, because $G(a,b) = G(b,a)$, this tells us $G(a,b) = 0$ when $b$ is odd as well. The final step requires us solve when both are even. We begin by assuming that the prior formula,

(36)$$G(2k,2j) = 2\frac{\Gamma\left(k+\frac{1}{2}\right)\Gamma\left(j+\frac{1}{2}\right)}{\Gamma\left(k+j+1\right)}$$

is true for all values from 0 to $k$ in both indices with $j\le k$. We know it’s true for $j=k=0$, which gives $2\pi$ for both formulas. We next proceed using our inductive assumption:

(37)$$\begin{aligned} G(2(k+1),2j) = \int_0^{2\pi} \textrm{d} \theta \sin^{2(k+1)}\theta \cos^{2j}\theta \\= \left.\frac{-1}{2j+1}\sin^{2k+1}\theta \cos^{2j+1}\theta\right|_0^{2\pi}+\frac{2k+1}{2j+1}\int_0^{2\pi} \textrm{d} \theta \sin^{2k}\theta \cos^{2(j+1)}\theta \\ =\frac{2k+1}{2j+1}\int_0^{2\pi} \textrm{d} \theta \sin^{2k}\theta\cos^{2j}\theta\left(1-\sin^{2}\theta\right)\\= \frac{2k+1}{2j+1}\left(G(2k,2j)-G(2(k+1),2j)\right) \end{aligned}$$

So that rearranging gives us

(38)$$\begin{aligned} G(2(k+1),2j)\left(1+\frac{2k+1}{2j+1}\right) = \frac{2k+1}{2j+1}G(2k,2j) \\ {\, \Longrightarrow \,} G(2(k+1),2j) = \frac{2k+1}{2(k+j+1)}G(2k,2j) \end{aligned}$$

which by the inductive assumption

(39)$$\begin{aligned} =2\frac{2k+1}{2(k+j+1)}\frac{\Gamma\left(\frac{2k+1}{2}\right)\Gamma\left(\frac{2j+1}{2}\right)}{\Gamma\left(k+j+1\right)} \\ =2\frac{\Gamma\left(\frac{2(k+1)+1}{2}\right)\Gamma\left(\frac{2j+1}{2}\right)}{\Gamma\left((k+1)+j+1\right)} \end{aligned}$$

proving the induction for $k+1$. Because $G(2k,2j) = G(2j,2k)$ this implies it for $j+1$ as well completing the induction. So that going back to our initial problem, we had

(40)$$\int_0^{2\pi} \textrm{d} \phi \cos(c\phi)\left(\sin\phi\right)^{b+g-2(B+G)}\left(\cos\phi\right)^{a+f-2(A+F)}$$

for ${}_{w_0,z_0,\lambda _0}\langle {a,b}\vert {c,d} \vert {f,g}\rangle _{w_i,z_i,\lambda _i}$ and

(41)$$\int_0^{2\pi} \textrm{d} \phi \sin(c\phi)\left(\sin\phi\right)^{b+g-2(B+G)}\left(\cos\phi\right)^{a+f-2(A+F)}$$

for ${}_{w_0,z_0,\lambda _0}\langle {a,b} \vert {-c,d} \vert {f,g}\rangle _{w_i,z_i,\lambda _i}$. Expanding out

(42)$$\begin{aligned} \cos(c \theta) = \frac{{\textrm{e}}^{ic\theta}+{\textrm{e}}^{-ic\theta}}{2} = \frac{\left(\cos\theta+i\sin\theta\right)^{c} +\left(\cos\theta -i\sin\theta\right)^{c}}{2} \\ =\sum_{j=0}^{\left\lfloor\frac{c}{2}\right\rfloor}\frac{(-1)^{j} c!}{(2j)!(c-2j)!}\cos^{c-2j}\theta\sin^{2j}\theta \end{aligned}$$

and through similar steps we can show

(43)$$\sin (c\theta ) = \sum_{j=0}^{\left\lfloor\frac{c-1}{2}\right\rfloor}\frac{(-1)^{j} c!}{(2j+1)!(c-2j-1)!}\cos^{c-2j-1}\theta \sin^{2j+1}\theta$$

yielding

(44)$$\begin{aligned} \int_0^{2\pi} \textrm{d} \phi \cos(c\phi)\left(\sin\phi\right)^{b+g-2(B+G)}\left(\cos\phi\right)^{a+f-2(A+F)} \\ = \sum_{D=0}^{\left\lfloor\frac{c}{2}\right\rfloor}\frac{(-1)^{D} c!}{(2D)!(c-2D)!}\int_0^{2\pi} \textrm{d} \phi\left(\sin\phi\right)^{b+g-2\left(B+G-D\right)}\left(\cos\phi\right)^{a+c+f-2\left(A+D+F\right)} \\ = \sum_{D=0}^{\left\lfloor\frac{c}{2}\right\rfloor}\frac{(-1)^{D} c!}{(2D)!(c-2D)!}\frac{\Gamma\left(\frac{b+g+1}{2}-(B+G-D)\right)\Gamma\left(\frac{a+c+f+1}{2}-(A+D+F)\right)}{\Gamma\left(\frac{a+b+c+f+g}{2}+1-(A+B+F+G)\right)} \end{aligned}$$

if both $b+g$ and $a+c+f$ are even, and 0 otherwise. For the other’s form

(45)$$\begin{aligned} \int_0^{2\pi} \textrm{d} \phi \sin(c\phi)\left(\sin\phi\right)^{b+g-2(B+G)}\left(\cos\phi\right)^{a+f-2(A+F)} \\ =\sum_{D=0}^{\left\lfloor\frac{c-1}{2}\right\rfloor} \frac{(-1)^{D} c!}{(2D+1)!(c-2D-1)!}\int_0^{2\pi} \textrm{d} \phi \left(\sin\phi\right)^{b+g+1-2(B+G-D)}\left(\cos\phi\right)^{a+c+f-1-2(A+D+F)} \\ =\sum_{D=0}^{\left\lfloor\frac{c-1}{2}\right\rfloor}\frac{(-1)^{D}c!}{(2D+1)!(c-2D-1)!}\frac{\Gamma\left(\frac{b+g}{2}+1-(B+G-D)\right)\Gamma\left(\frac{a+c+f}{2}-(A+D+F)\right)}{\Gamma\left(\frac{a+b+c+f+g}{2}+1-(A+B+F+G)\right)} \end{aligned}$$

when both $b+g$ and $a+c+f$ are odd but 0 otherwise.

B.2. Back to our problem

Using this we have

(46)$$\begin{aligned} {}_{w_0,z_0,\lambda_0}\langle {a,b}\vert {c,d} \vert {f,g}\rangle _{w_i,z_i,\lambda_i}=\chi\left. \sum_{A=0}^{\left\lfloor\frac{a}{2}\right\rfloor}\sum_{B=0}^{\left\lfloor\frac{b}{2}\right\rfloor}\sum_{C=0}^{\left\lfloor\frac{d-c}{2}\right\rfloor}\sum_{D=0}^{\left\lfloor\frac{c}{2}\right\rfloor}\sum_{F=0}^{\left\lfloor\frac{f}{2}\right\rfloor}\sum_{G=0}^{\left\lfloor\frac{g}{2}\right\rfloor}\right[\int_0^{r} \textrm{d} \rho \rho^{a+b+d+f+g-2(A+B+C+F+G)+1}{\textrm{e}}^{-\gamma \rho^{2}} \\ \left. \times\frac{(-1)^{A+B+C+D+F+G}(d-C)!W_i^{2(F+G)}W_o^{2(A+B)}r^{2C}}{8^{A+B+F+G}A!B!C!F!G!(a-2A)!(b-2B)!\left(\frac{d+c}{2}-C\right)!\left(\frac{d-c}{2}-C\right)!(f-2F)!(g-2G)!(2D)!(c-2D)!}\right. \\ \times\left.\frac{\Gamma\left(\frac{a+c+f+1}{2}-(A+D+F)\right)\Gamma\left(\frac{b+g+1}{2}-(B-D+G)\right)}{\Gamma\left(\frac{a+b+c+f+g}{2}+1-(A+B+F+G)\right)}\right] \end{aligned}$$

where we have denoted

(47)$$\chi = \frac{c!\sqrt{a!b!f!g!}2^{a+b+f+g+1}{\textrm{e}}^{i\left((f+g+1)\psi_i +k_o\delta z_o-(a+b+1)\psi_o-k_i\delta z_i\right)}}{\pi W_o^{a+b+1}W_i^{f+g+1}r^{d}}$$

Note that for this case, as well as the case with negative upper Zernike index, we have that the outer integral was 0 unless $a+b+c+f+g$ was even, and because $d-c$ is even this also means that this is only nonzero when $a+b+d+f+g$ is even. This tells us the total exponent of $\rho$ can be written in the form $2m+1$ for some $m$, yielding the next integral to solve as

(48)$$\int_0^{r} \textrm{d} \rho \rho^{2m+1}{\textrm{e}}^{-\gamma \rho^{2}}$$

B.3. Integral 2 solution

Defining

(49)$$H(m) \equiv \int_0^{r} \textrm{d} \rho \rho^{2m+1}{\textrm{e}}^{-\gamma \rho^{2}}$$

we have

(50)$$=\left. -\frac{\rho^{2m}{\textrm{e}}^{-\gamma \rho^{2}}}{2\gamma}\right|_{\rho=0}^{r}+\frac{m}{\gamma}\int_0^{r} \textrm{d} \rho \rho^{2(m-1)+1}{\textrm{e}}^{-\gamma\rho^{2}} = \frac{-r^{2m}{\textrm{e}}^{-\gamma r^{2}}}{2\gamma}+\frac{m}{\gamma}H(m-1)$$

This form, valid as long as $m\ge 1$ hints

(51)$$H(m) = P(k)\equiv\frac{m!}{\gamma^{k}(m-k)!}H(m-k)-\frac{m!r^{2(m+1)}{\textrm{e}}^{-\gamma r^{2}}}{2}\sum_{T=1}^{k}\frac{1}{(m+1-T)!(\gamma r^{2})^{T}}$$

for any $k<m$This is clearly true for $k=0$ and for $k=1$ we get back the relation we just saw. We assume $H(m) = P(k)$ is true for some $k+1\le m$ which gives, using that $m-k\ge 1$ so that our initial condition is still valid,

(52)$$\begin{aligned} H(m) &= P(k) = \frac{m!}{\gamma^{k}(m-k)!}H(m-k)-\frac{m!r^{2(m+1)}{\textrm{e}}^{-\gamma r^{2}}}{2}\sum_{T=1}^{k}\frac{1}{(m+1-T)!(\gamma r^{2})^{T}} \\ &=\frac{m!}{\gamma^{k}(m-k)!}\left(-\frac{r^{2(m-k)}{\textrm{e}}^{-\gamma r^{2}}}{2\gamma}+\frac{(m-k)}{\gamma}H(m-(k+1))\right)-\frac{m! r^{2(m+1)}{\textrm{e}}^{-\gamma r^{2}}}{2}\sum_{T=1}^{k}\frac{1}{(m+1-T)!(\gamma r^{2})^{T}} \\ &=\frac{m!H(m-(k+1))}{\gamma^{k+1}(m-(k+1))!}-\frac{m!r^{2(m+1)}{\textrm{e}}^{-\gamma r^{2}}}{2(m+1-(k+1))!\left(\gamma r^{2}\right)^{k+1}}-\frac{m! r^{2(m+1)}{\textrm{e}}^{-\gamma r^{2}}}{2}\sum_{T=1}^{k}\frac{1}{(m+1-T)!(\gamma r^{2})^{T}} \\ &=\frac{m!}{\gamma^{k+1}(m-(k+1))!}H(m-(k+1))-\frac{m! r^{2(m+1)}{\textrm{e}}^{-\gamma r^{2}}}{2}\sum_{T=1}^{k+1}\frac{1}{(m+1-T)!(\gamma r^{2})^{T}} = P(k+1) \end{aligned}$$

finishing the proof that $H(m) = P(k)$ for all $k\le m$. Using this and setting $k=m$ we have

(53)$$H(m) = \frac{m!}{\gamma^{m}}H(0)-\frac{m!r^{2(m+1)}{\textrm{e}}^{-\gamma r^{2}}}{2}\sum_{T=1}^{m} \frac{1}{(m+1-T)!(\gamma r^{2})^{T}}$$

As for the final evaluation,

(54)$$H(0) = \int_0^{r} \textrm{d} \rho \rho {\textrm{e}}^{-\gamma \rho^{2}} = \left.-\frac{{\textrm{e}}^{-\gamma \rho^{2}}}{2\gamma}\right|_{\rho=0}^{r} = \frac{1}{2\gamma}-\frac{{\textrm{e}}^{-\gamma r^{2}}}{2\gamma}$$

so that

(55)$$H(m) = \frac{m!}{2\gamma^{m+1}}-\frac{m! r^{2(m+1)}{\textrm{e}}^{-\gamma r^{2}}}{2}\sum_{T=1}^{m+1}\frac{1}{(m+1-T)!(\gamma r^{2})^{T}}$$

By making the change of variables in the summation to $K = (m+1)-T$ we rewrite this as

(56)$$\begin{aligned} H(m) = \frac{m!}{2\gamma^{m+1}}-\frac{m! r^{2(m+1)} {\textrm{e}}^{-\gamma r^{2}}}{2}\sum_{K=0}^{m} \frac{1}{K!\left(\gamma r^{2}\right)^{m+1-K}} \\ =\frac{m!}{2\gamma^{m+1}}\left(1-{\textrm{e}}^{-\gamma r^{2}}\sum_{K=0}^{m}\frac{\left(\gamma r^{2}\right)^{K}}{K!}\right) \end{aligned}$$

B.4. Incorporating the second integral into the final expression

Using this with our last expression we have the final solution

(57)$$\begin{aligned} {}_{w_o,z_o,\lambda_o}\langle {a,b}\vert {c,d} \vert{f,g}\rangle _{w_i,z_i,\lambda_i} =\chi\left. \sum_{A=0}^{\left\lfloor\frac{a}{2}\right\rfloor}\sum_{B=0}^{\left\lfloor\frac{b}{2}\right\rfloor}\sum_{C=0}^{\left\lfloor\frac{d-c}{2}\right\rfloor}\sum_{D=0}^{\left\lfloor\frac{c}{2}\right\rfloor}\sum_{F=0}^{\left\lfloor\frac{f}{2}\right\rfloor}\sum_{G=0}^{\left\lfloor\frac{g}{2}\right\rfloor}\right[\frac{\int_0^{r} \textrm{d} \rho \rho^{a+b+d+f+g-2(A+B+C+F+G)+1}{\textrm{e}}^{-\gamma \rho^{2}}}{A!B!C!F!G!(a-2A)!(b-2B)!} \\ \left. \frac{(-1)^{A+B+C+D+F+G}(d-C)!W_i^{2(F+G)}W_o^{2(A+B)}r^{2C}\Gamma\left(\frac{a+c+f+1}{2}-(A+D+F)\right)\Gamma\left(\frac{b+g+1}{2}-(B-D+G)\right)}{8^{A+B+F+G}\left(\frac{d+c}{2}-C\right)!\left(\frac{d-c}{2}-C\right)!(f-2F)!(g-2G)!(2D)!(c-2D)!\Gamma\left(\frac{a+b+c+f+g}{2}+1-(A+B+F+G)\right)}\right] \\ =\frac{\chi}{2\gamma^{\frac{a+b+d+f+g}{2}+1}}\sum_{A=0}^{\left\lfloor\frac{a}{2}\right\rfloor}\sum_{B=0}^{\left\lfloor\frac{b}{2}\right\rfloor}\sum_{C=0}^{\left\lfloor\frac{d-c}{2}\right\rfloor}\sum_{D=0}^{\left\lfloor\frac{c}{2}\right\rfloor}\sum_{F=0}^{\left\lfloor\frac{f}{2}\right\rfloor}\sum_{G=0}^{\left\lfloor\frac{g}{2}\right\rfloor} \\ \left[\frac{\Gamma\left(\frac{a+b+d+f+g}{2}+1-(A+B+C+F+G)\right)}{\Gamma\left(\frac{a+b+c+f+g}{2}+1-(A+B+F+G)\right)}\left(1-{\textrm{e}}^{-\gamma r^{2}}\sum_{K=0}^{\frac{a+b+d+f+g}{2}-(A+B+C+F+G)}\frac{(\gamma r^{2})^{K}}{K!}\right)\right. \\ \left. \frac{(-1)^{D}(-\gamma)^{A+B+C+F+G}(d-C)!W_i^{2(F+G)}W_o^{2(A+B)}r^{2C}\Gamma\left(\frac{a+c+f+1}{2}-(A+D+F)\right)\Gamma\left(\frac{b+g+1}{2}-(B-D+G)\right)}{8^{A+B+F+G}A!B!C!F!G!(a-2A)!(b-2B)!\left(\frac{d+c}{2}-C\right)!\left(\frac{d-c}{2}-C\right)!(f-2F)!(g-2G)!(2D)!(c-2D)!}\right] \end{aligned}$$

or redefining

(58)$$\sigma \equiv \frac{\chi}{2\gamma^{\frac{a+b+d+f+g}{2}+1}} = \frac{c!\sqrt{a!b!f!g!}2^{a+b+f+g}{\textrm{e}}^{i\left((f+g+1)\psi_i+k_o\delta z_o-(a+b+1)\psi_o-k_i\delta z_i\right)}}{\pi W_o^{a+b+1}W_i^{f+g+1}r^{d}\sqrt{\gamma^{a+b+d+f+g+2}}}$$

this gives back the general expression from the first section of the appendix,

(59)$$\begin{aligned}{}_{w_0,z_0,\lambda_0}\langle {a,b}\vert {c,d} \vert{f,g} \rangle _{w_i,z_i,\lambda_i}= \left.\sigma\sum_{A=0}^{\left\lfloor\frac{a}{2}\right\rfloor}\sum_{B=0}^{\left\lfloor\frac{b}{2}\right\rfloor}\sum_{C=0}^{\left\lfloor\frac{d-c}{2}\right\rfloor}\sum_{D=0}^{\left\lfloor\frac{c}{2}\right\rfloor}\sum_{F=0}^{\left\lfloor\frac{f}{2}\right\rfloor}\sum_{G=0}^{\left\lfloor\frac{g}{2}\right\rfloor}\right. \\ \left[\frac{\Gamma\left(\frac{a+b+d+f+g}{2}+1-(A+B+C+F+G)\right)}{\Gamma\left(\frac{a+b+c+f+g}{2}+1-(A+B+F+G)\right)}\right. \\ \times\left(1-{\textrm{e}}^{-\gamma r^{2}}\sum_{K=0}^{\frac{a+b+d+f+g}{2}-(A+B+C+F+G)}\frac{(\gamma r^{2})^{K}}{K!}\right) \\ \times\frac{(-1)^{D}(-\gamma)^{A+B+C+F+G}W_i^{2(F+G)}W_o^{2(A+B)}r^{2C}}{8^{A+B+F+G}A!B!C!F!G!(a-2A)!(b-2B)!} \\ \left.\times\frac{(d-C)!\Gamma\left(\frac{a+c+f+1}{2}-(A+D+F)\right)\Gamma\left(\frac{b+g+1}{2}-(B-D+G)\right)}{\left(\frac{d+c}{2}-C\right)!\left(\frac{d-c}{2}-C\right)!(f-2F)!(g-2G)!(2D)!(c-2D)!}\right] \end{aligned}$$

The same work, but using the expressions for the negative upper Zernike index from the Integral 1 Proof subsection within this section of the appendix (which leads to changing the fourth summation upper limit and a few of the factorial expressions), leads to the second general expression from the first portion of this appendix.

C. Appendix: Coefficient table

Table of coefficients to build clipped Gaussian beams of section 3.

Table 2. Table of the Non-Zero Mode Coefficients. These happen to be real because our initial distributions were real and these were fit to modes at their waist. These coefficients do not sum to unity because they represent a fit to a unit power beam that has been clipped.

View Table | View all tables in this article

D. Appendix: Fractional RMS deviation of an HG-fit

We define the Fractional RMS of a Beam $U$ from it’s fit $U_{fit}$ as

(60)$$FRMS = \sqrt{\frac{\iint \left \lvert {U-U_{fit}} \right \rvert^{2}}{\iint \left \lvert {U} \right \rvert^{2}}}$$

Because the HG modes are complete over the 2-D plane we can then write out

(61)$$U = \sum_{(m,n)\in A\cup B}a_{m,n}\textrm{HG}_{m,n}$$

where $A$ is defined as the set of mode indices that we have been able to include in our fit, i.e. $A\subset \mathbb {N}_0^{2}$ (with $\mathbb {N}_0 = \mathbb {Z}^{+}\cup \{0\}$) is determined by

(62)$$U_{fit} = \sum_{(m,n)\in A}a_{m,n}\textrm{HG}_{m,n}\qquad B = \mathbb{N}_0^{2}\setminus A$$

Now due to orthogonality properties of HG modes we have the power in any beam

$R = \sum _{m,n} c_{m,n}\textrm {HG}_{m,n}$ is given by

(63)$$\begin{aligned}P_R &= \iint R R^{*} = \iint \left(\sum_{m,n,j,k}c_{m,n}c_{j,k}^{*}\textrm{HG}_{m,n}\textrm{HG}_{j,k}^{*}\right) \end{aligned}$$

(64)$$\begin{aligned}&=\sum_{m,n,j,k}c_{m,n}c_{j,k}^{*}\iint \textrm{HG}_{m,n}\textrm{HG}_{j,k}^{*} = \sum_{m,n,j,k}c_{m,n}c_{j,k}^{*}\delta_{m,j}\delta_{n,k} \end{aligned}$$

(65)$$\begin{aligned}&= \sum_{m,n}\left \lvert {c_{m,n}} \right \rvert^{2} \end{aligned}$$

or is just the squared absolute sum of the beam coefficients. So suppose the fraction of the total power we’ve managed to represent in our beam is given by $F$, i.e.

(66)$$F = \frac{\sum_{(m,n)\in A} \left \lvert {a_{m,n}} \right \rvert^{2}}{\sum_{(m,n)\in A\cup B}\left \lvert {a_{m,n}} \right \rvert^{2}}=\frac{\sum_{(m,n)\in A} \left \lvert {a_{m,n}} \right \rvert^{2}}{\sum_{(m,n)\in A}\left \lvert {a_{m,n}} \right \rvert^{2}+\sum_{(m,n)\in B}\left \lvert {a_{m,n}} \right \rvert^{2}}$$

Then we get for our fractional RMS

(67)$$\begin{aligned}FRMS &= \sqrt{\frac{\iint \left \lvert {\sum_{(m,n)\in A\cup B}a_{m,n}\textrm{HG}_{m,n}-\sum_{(m,n)\in A}a_{m,n}\textrm{HG}_{m,n}} \right \rvert^{2}}{\iint\left \lvert {\sum_{(m,n)\in A\cup B}a_{m,n}\textrm{HG}_{m,n}} \right \rvert^{2}}} \end{aligned}$$

(68)$$\begin{aligned}&=\sqrt{\frac{\iint\left \lvert {\sum_{(m,n)\in B}a_{m,n}\textrm{HG}_{m,n}} \right \rvert^{2}}{\sum_{(m,n)\in A\cup B}\left \lvert {a_{m,n}} \right \rvert^{2}}} = \sqrt{\frac{\sum_{(m,n)\in B}\left \lvert {a_{m,n}} \right \rvert^{2}}{\sum_{(m,n)\in A\cup B}\left \lvert {a_{m,n}} \right \rvert^{2}}} \end{aligned}$$

(69)$$\begin{aligned}&= \sqrt{\frac{\sum_{(m,n)\in A\cup B}\left \lvert {a_{m,n}} \right \rvert^{2}-\sum_{(m,n)\in A}\left \lvert {a_{m,n}} \right \rvert^{2}}{\sum_{(m,n)\in A\cup B}\left \lvert {a_{m,n}} \right \rvert^{2}}} = \sqrt{1-\frac{\sum_{(m,n)\in A}\left \lvert {a_{m,n}} \right \rvert^{2}}{\sum_{(m,n)\in A\cup B}\left \lvert {a_{m,n}} \right \rvert^{2}}} \end{aligned}$$

(70)$$\begin{aligned}&=\sqrt{1-F} \end{aligned}$$

Telling us that if the fraction of the total power our fit holds is $F$ then the fractional RMS is equal to $\sqrt {1-F}$, giving us an idea of how much our fit can deviate from the real distribution.

Funding

National Aeronautics and Space Administration (80NSSC18K0552).

Acknowledgments

We’d like to acknowledge Ada Uminska for her part in helping review this paper.

Disclosures

The authors declare no conflicts of interest.

References

1. A. E. Siegman, “Lasers,” (University Science Books, 1986), Ch. 16, pp. 626–648

2. C. Bond, D. Brown, A. Freis, and K. Strain, “Interferometer Techniques for Gravitational-Wave Detection,” arXiv:0909.3661 (2015), Ch. 9, pp. 114–130.

3. C. Bond, D. Brown, A. Freis, and K. Strain, “Interferometer Techniques for Gravitational-Wave Detection,” arXiv:0909.3661, (2015), Ch. 11.9, pp. 182–184.

4. The Virgo Collaboration, “The Virgo Physics Book, Vol. II: Optics and Related Topics”, (2005), Ch. 3.3, pp. 202–211.

5. S. Schuster, “Tilt-to-Length coupling and diffraction aspects in satellite interferometry,” (2016) Ch. 7.3.3, pp. 162.

6. C. Bond, D. Brown, A. Freis, and K. Strain, “Interferometer Techniques for Gravitational-Wave Detection,” arXiv:0909.3661, (2015), Ch. 9.16, pp. 134.

7. C. Bond, D. Brown, A. Freis, and K. Strain, “Interferometer Techniques for Gravitational-Wave Detection,”, arXiv:0909.3661 (2015), Ch. 11, pp. 166–186.

8. C. Bond, D. Brown, A. Freis, and K. Strain, “Interferometer Techniques for Gravitational-Wave Detection”, arXiv:0909.3661, (2015), Ch. 11.7, pp. 172–174.

9. P. Amaro-Seoane, H. Audley, S. Babak, J. Baker, E. Barausse, P. Bender, E. Berti, P. Binetruy, M. Born, D. Bortoluzzi, J. Camp, C. Caprini, V. Cardoso, M. Colpi, J. Conklin, N. Cornish, C. Cutler, K. Danzmann, R. Dolesi, L. Ferraioli, V. Ferroni, E. Fitzsimons, J. Gair, L. G. Bote, D. Giardini, F. Gibert, C. Grimani, H. Halloin, G. Heinzel, T. Hertog, M. Hewitson, K. Holley-Bockelmann, D. Hollington, M. Hueller, H. Inchauspe, P. Jetzer, N. Karnesis, C. Killow, A. Klein, W. Klipstein, N. Korsakova, S. L. Larson, J. Livas, I. Lloro, N. Man, D. Mance, J. Martino, K. McKenzie, S. T. McWilliams, C. Miller, G. Mueller, G. Nardini, G. Nelemans, M. Nofrarias, A. Petiteau, P. Pivato, E. Plagnol, E. Porter, J. Reiche, D. Robertson, N. Robertson, E. Rossi, G. Russano, B. Schutz, A. Sesana, D. Shoemaker, J. Slutsky, C. Sopuerta, T. Sumner, N. Tamanini, J. I. Thorpe, M. Troebs, M. Vallisneri, A. Vecchio, D. Vetrugno, S. Vitale, G. Wanner, H. Ward, P. Wass, W. Weber, J. Ziemer, and P. Zweifel, “LISA Laser Interferometer Space Antenna, A proposal in response to the ESA call for L3 missions,” (2017).

10. R. G. Schell and G. Tyras, “Irradiance from an Aperture with a Truncated-Gaussian Field Distribution,” J. Opt. Soc. Am. 61(1), 31 (1971). [CrossRef]  

11. C. Campbell, “Fresnel Diffraction of Gaussian Laser Beams by Circular Aperture,” Opt. Eng. 26(3), 263270 (1987). [CrossRef]  

12. G. Lenz, “Far-Field Diffraction of Truncated Higher-Order Laguerre-Gaussian Beams,” Opt. Commun. 123(4-6), 423–429 (1996). [CrossRef]  

13. K. Tanaka, M. Shibukawa, and O. Fukumitsu, “Diffraction of a Wave Beam by an Aperture,” IEEE Trans. Microwave Theory Tech. 20(11), 749–755 (1972). [CrossRef]  

14. J. D. Jackson, Classical Electrodynamics: Third Edition, (John Wiley & Sons, 1999), Ch 10.5, pp. 478–491.

15. F. Bayer-Helms, “Coupling coefficients of an incident wave and the modes of spherical optical resonator in the case of mismatching and misalignment,” Appl. Opt. 23(9), 1369–1380 (1984). [CrossRef]  

16. C. Bond, D. Brown, A. Freis, and K. Strain, “Interferometer Techniques for Gravitational-Wave Detection,” arXiv:0909.3661, (2015), Ch. 11.3, pp. 168.

17. J. Huang and D. Ding, “An Alternative Analytical Description of Truncated Gaussian Beams,” J. Opt. A: Pure Appl. Opt. 10(1), 015003 (2008). [CrossRef]