1.1.

Astrophysical Motivation

Integrated models of optical observatories are highly beneficial to their design and use.^1,2 Accurate observatory models permit powerful insight into predicting the as-built performance of a given instrument. However, the accuracy of these models is fundamentally limited by the assumptions made. To facilitate high-yield scientific observations, astronomical observatories are nominally designed to operate in the diffraction limit where wavefront aberrations are small. Diffraction-limited optical observatories are necessarily modeled with diffraction integrals derived from the Huygens–Fresnel principle to support the wave-like behavior of light. The paraxial and scalar assumptions that angles of incidence are small and that polarization is negligible³ are made to ease the computational burden on the model. The resultant Fresnel and Fraunhofer diffraction integrals are accurate, providing these conditions are met. If the performance of the observatory is limited by a factor outside the assumptions made, and then the model will be ignorant of it. An example of this is the linear and shift-invariant assumption imposed on diffraction models of astronomical observatories. Ray aberrations (e.g., coma and astigmatism) have a field dependence and consequently change across an observatory’s field of view. However, diffraction integrals assume shift invariance. This means that the aberrations do not change across the field of view and a separate ray trace model must be used to capture this effect.

Integrating optical models from different regimes in physics has become a popular method by which to overcome this limitation. Linking ray trace models to diffraction models in particular can overcome the paraxial and scalar assumption imposed by the Fresnel and Fraunhofer diffraction integrals. In the prior example, to capture the influence of optical aberrations, a new ray trace must be performed and the optical path difference of the rays must be translated to a diffraction model for each point of interest in the field of view. Similarly, diffraction integrals are incapable of determining the effects of optical polarization. For example, the Daniel K. Inoyue Solar Telescope supports a suite of polarimetric instrumentation that is sensitive to the influence of optical polarization. To support this regime of optical physics, the scalar assumption is not sufficient, so the polarization state is propagated along geometric ray paths using polarization ray tracing^4,5 to determine the influence of polarization aberrations on the optical beam.⁶ Modern space telescopes also require integrated models to accurately predict the instrument behavior. For example, the optical models of the James Webb Space Telescope incorporate the influence of dynamic thermal, structural, and optical effects simultaneously to produce an accurate library of the observatory’s jitter.⁷ High-contrast imaging instruments that use coronagraphs, designed to separate exoplanets from diffracted starlight, are in dire need of integrated physical optics modeling from their inception. These coronagraphs aim to discern targets that are orders of magnitude dimmer than their host star.^8–10 Understanding all sources of error is of paramount importance to the functionality of the instrument.

High-contrast imaging instruments have been successfully deployed on the ground (e.g., SCExAO,¹¹ MagAO-X,¹² NIRC2,¹³ GPI,¹⁴ and SPHERE¹⁵) and in space (e.g., NICMOS¹⁶ and NIRCam^17,18) to pursue the direct detection of extrasolar planets, debris, and protoplanetary disks. The Decadal Survey on Astronomy and Astrophysics 2020 (Astro2020) recommends pursuing these instruments for a future 6 meter diameter infrared/optical/visible (IROUV) flagship observatory for the progression of astrophysical sciences.¹⁹

Presently, the optical design of observatories is done in a ray-tracing engine^20,21 (e.g., CODE V and Zemax OpticStudio) because it is more suitable to optimizing the shapes of observatory mirrors. Upon reaching a diffraction-limited optical design, the system is then assumed to be well-represented by a paraxial diffraction model. The wavefront maps produced by the ray trace model of the observatory and the contributions from the imperfect polishing of the observatory mirrors are sent to a linearized physical optics propagator to examine the image plane electric field in the presence of diffraction from structure in the beam and phase errors on the optics.

Many tools have been developed to simulate the performance of high-contrast imaging instrumentation. Tiny Tim is one of the first of these widely used packages used to simulate the Hubble Space Telescope (HST) instrument point-spread functions (PSFs).²² The tool generates aberrated PSFs based on the instrument, observation, and dynamic aberrations for a given observing scenario, enabling highly accurate simulations of the observatory performance. However, it only considers aberrations that are conjugate to the exit pupil of the observatory. This limits the model’s ability to capture out-of-pupil effects, such as the Talbot effect and speckles from optical surfaces.³ To capture these effects, optical propagation packages integrate optical models of observatories by adding Fresnel diffraction to the PSF simulation, enabling the modeling of plane-to-plane diffraction effects. Open-source packages that currently support Fresnel diffraction include: PROPER,²³ Physical Optics Propagation in PYthon (POPPY),^24–26 High-Contrast Imaging in Python (HCIPy),²⁷ AOTools,²⁸ and prysm.^2,29 Using these tools, near-field diffraction that limits high-contrast imaging can be modeled, and focal-plane wavefront sensing methods can be tested.

These open-source physical optics propagation tools form the cornerstone of high-contrast imaging instrument modeling and design. The open-source framework means that the codes are accessible to anyone, so the physics are completely verifiable by the scientific community.³⁰ It is in the scientific community’s best interest to continue to develop open-source propagation physics modules to increase the scope of and further integrate our observatory models.

Commercial optical design codes offer the ability to make diffraction calculations based on ray data, but their physical optics simulation techniques are not as transparent or versatile as the open-source propagation codes that are used to design coronagraphs for astronomical observatories. The current open-source physical optics codes used for observatory modeling are also limited in their scope because of the Fresnel approximation, which is incapable of accurately modeling the field after fast-focusing and highly aspheric surfaces.^31,32 As observatories get larger, their optics may become faster (i.e., lower $F\#$) and more aspheric to fit within an available volume. Some coronagraph architectures capable of Earth-like exoplanet detection (e.g., phase-induced amplitude apodization or PIAA³³) employ mirrors that apodize the pupil with highly aspheric mirrors and require tailored propagators in order to be included in physical optics models.³¹ In the regime where the contribution of these surfaces is best represented by a ray trace, a diffraction calculation must be made to appropriately model the optical field at the image plane. To continue the development of integrated optical models, exploring the possibilities and limitations of new propagation techniques is desirable.

An example of the typical integrated modeling pipeline for astronomical observatories outfitted with coronagraphs is shown in Fig. 1. The observatory is typically designed and modeled in ray trace software (e.g., CODE V and Zemax OpticStudio) to accurately model the wavefront in the observatory’s exit pupil. Upon finalizing the design, the complex-valued exit pupil is decomposed into a functional representation (e.g., a set of polynomial coefficients, such as the Zernikes³⁵) and passed to the entrance pupil of a coronagraph model constructed in an open-source, Fourier-based physical optics propagator (e.g., POPPY, PROPER, and HCIPy). The front-end model computes the complex field distribution at the coronagraph mask and then propagates the field past the mask to the image plane. The field is taken by a model of the detector (e.g., EMCCD Detect³⁶ and Pyxel³⁷) to create a simulated raw science image that can be postprocessed (e.g., PyKLIP³⁸ and NMF imaging³⁹).

Fig. 1

Modeling flow inspired by the structural thermal optical performance (STOP) modeling process for the Roman Coronagraph.³⁴ This diagram illustrates the different modeling regimes required to create a simulated image of an observatory. We aim to further integrate this modeling pipeline by creating an open-source GBD platform to unify the ray trace model of the observatory with the diffraction model of the coronagraph.

To bridge the gap between commercial ray tracing engines and open-source physical optics propagation codes, we investigate the viability of a ray-based diffraction calculation called Gaussian beamlet decomposition (GBD) for modeling observatories with coronagraphs. Traditionally, GBD operates using the complex ray tracing algorithm described in the works by Greynolds⁴⁰ and Harvey et al.⁴¹ This technique has been previously implemented in FRED,^41,42 and possibly in CODE V,⁴³ but an exact method of its implementation in these software packages is not clearly available in the literature. An alternative approach called the transfer matrix algorithm was recently developed to improve GBD’s viability for precision diffraction simulation by Worku and Gross.^44–46 However, their implementation is not public and has not yet been formally evaluated as a tool to augment the modeling of astronomical observatories or high-contrast imaging instrumentation. To formally evaluate GBD as a modeling tool for astronomical instrumentation, a complete algorithm for its implementation is derived in this article.

1.2.

Gaussian Beamlet Decomposition

GBD is a method of physical optics propagation that approximates the propagated field as a finite sum of coherent Gaussian beams that each propagate along a ray path. This method has been implemented in optical design packages⁴⁷ to perform coherent calculations on non-paraxial systems. The operating principle of GBD is to decompose the field in the entrance pupil of an optical system [Fig. 2(a)] into a finite set of Gaussian beams [Fig. 2(b)]. Their coherent sum [Figs. 2(c) and 2(d)] approximates the initial field decomposition and can be propagated anywhere in the optical system along geometric ray paths.

Fig. 2

Illustration of the operating principle of GBD in one dimension. (a) The aperture function for traditional imaging systems as a top-hat function; (b) the decomposition of this function into a discrete set of Gaussian beams; (c) nine evenly spaced Gaussian profiles before their coherent summation. The coherent sum (black) of these nine beamlets shows that the beamlets are incapable of perfectly reconstructing the original aperture function (dashed red), specifically the sharp edge and uniform amplitude. (d) Using more beamlets results in a more accurate reconstruction of the aperture function. The amplitude ripple has virtually vanished, and the field near the aperture edges is almost entirely recovered.

Undersampling the field in the entrance pupil leads to artifacts in the decomposition. A characteristic amplitude ripple based on the period of the beamlet decomposition remains in the field. Due to the soft edges of the Gaussian beams, they cannot completely reconstruct the field of a sharp aperture edge. Using a larger number of smaller beamlets decreases the beamlet distribution period and increases the slope of the Gaussian beams, allowing for the mitigation of both of these effects [shown in Fig. 2(d)].

Upon decomposing the initial wavefront into a sufficient set of Gaussian beams, we can use their analytical linear propagation laws⁴⁸ (described in Sec. 2) to compute the coherent field at any arbitrary plane in the optical system. Fourier transform-based propagation methods derived from the Huygens–Fresnel principle typically assume that the field is scalar and the optical system is paraxial.³ This is usually appropriate for stellar coronagraphs, which operate on slowly focusing beams with diffraction-limited optics, but may not be for next-generation observatories whose large apertures may necessitate relatively fast telescope optics.

GBD computes the same complex optical field without making the paraxial assumption across the observatory. Rather, the coherent field of Gaussian beams is derived from the ray data directly. Doing so imposes the paraxial assumption about a single beamlet instead of the entire observatory, which is a much less stringent approximation. Gaussian beams are technically infinite in extent, but extremely localized around the beam waist. Consequently, the contribution of the field very far from the Gaussian is negligible. This locality enables the simulation of the optical system to generally be non-paraxial. By making the diffraction calculation directly from ray data, GBD circumvents the need for translating the wavefront to a physical optics propagator and imposing the paraxial assumption on the optical system. Instead, the ray trace model can be directly integrated into the diffraction model.

There are two main approaches that exist in the literature to implement GBD: the complex ray tracing method and the transfer matrix method. The complex ray tracing method was recently described by Harvey et al.⁴¹ in their seminal paper about implementing GBD in photon engineering’s non-sequential ray tracing software FRED. This method traces waist and divergence rays to compute the complex field at the plane of interest using Arnaud’s method of complex ray tracing.^40,49 Through FRED, the complex ray tracing method has seen widespread use for non-paraxial coherent beam analysis.

Another GBD approach, the transfer matrix method was developed by Worku and Gross^44–46 to mitigate GBD’s inability to simulate sharp-edge diffraction and add new utility to the technique. This formulation of GBD works by computing the differential ray transfer matrix for a given ray path and then using that data to solve the Gaussian beam solution to the general Collins integral.⁵⁰ Worku and Gross have leveraged the general Collins integral to provide alternative conditions to the Gaussian beam solution to modify the decomposition, such as truncated⁴⁶ and pulsed⁵¹ beamlet decomposition. The option of modifying the beamlets to overcome the limitations of GBD makes the transfer matrix method extremely attractive for use in high-contrast imaging where preservation of high-spatial frequency content is important. Of particular interest are mirror segment gaps and optomechanical structures that obscure the primary mirror. Therefore, we elect to investigate Worku and Gross’s transfer matrix method of GBD for the work presented in this article. Our goal is to publicize the transfer matrix method by developing an algorithm for its implementation and then use it to characterize GBD’s suitability for high-contrast imaging simulation. Our work can then be used as a platform with which to study the suitability of Worku and Gross’ modified GBD in future investigations.

1.3.

Hybrid Propagation Physics

The ray-based nature of GBD introduces problems in modeling the electric field when the rays are vignetted. Structure in the field where the initial decomposition occurs (typically, the entrance pupil) can be well-represented by Gaussian beams as long as the structure of interest is larger than the beamlets used in decomposition. Diffraction from structure in interemediate planes (between the pupil and focal planes) is challenging to represent if the beamlets diverge considerably. However, secondary beamlets can be traced from these intermediate structures to aid in the accuracy of the simulation.⁵² At the focal plane of a diffraction-limited system, the rays are highly concentrated while the diffracted field spreads out considerably (e.g., the Airy disk). Because of GBD’s reliance on ray tracing, it cannot represent diffraction from structure in the focal plane well without redecomposing the field.⁴² For the case of a Lyot-type coronagraph, all rays are vignetted at the focal plane mask (FPM) and the field decomposition is lost. To circumvent this, we compute the field before the FPM with GBD and propagate it through the remaining coronagraph with traditional diffraction integrals, where we expect the low-order aberrations to be small and the paraxial assumption to be valid. This hybrid method (shown in Fig. 3) enables the user to alter the propagation physics for the electric field based on where it is the most appropriate; GBD will simulate the fast beams in the fore-optics and paraxial diffraction will propagate the field through the coronagraph imaging optics. In practice, GBD would be used to propagate to the entrance pupil of the coronagraph for simulating systems with wavefront control, where it would then hand off the field to a paraxial diffraction model.

Fig. 3

Diagram demonstrating the hybrid propagation physics model. The observatory optics are best described by a ray-based propagation model, so GBD is used to compute the field at the image plane before the coronagraph FPM. The field array is then passed to a paraxial diffraction model, which propagates it past the FPM and through the remainder of the coronagraph. This propagation scheme permits the user to model the influence of the observatory optics non-paraxially, without losing accuracy after propagating past the FPM.

The end result of such a model allows for direct integration of the ray trace model with the physical optics model, without imposing the paraxial approximation on the observatory. Like Fresnel and Fraunhofer diffraction, GBD is an approximation to diffraction physics. The decomposition of the field into Gaussian beams does not have an analytical solution. Therefore, undesirable artifacts can be introduced into the field if the decomposition is not well-understood and the sampling is insufficient.⁵³ To better understand the impact of a GBD PSF on high-contrast imaging simulations, we develop a hybrid propagation model to compare GBD to an equivalent Fraunhofer diffraction model.

To our knowledge, the transfer matrix method of GBD has not seen widespread implementation. Given its obvious benefits, we believe that this is because the transfer matrix method is not well-understood by the scientific community. We aim to remedy this by presenting a vectorizeable algorithm for the transfer matrix method and open-sourcing our simulation platform for future investigators to use. This article is the first work, to our knowledge, to provide the explicit mathematics of GBD and provide our code as an object-oriented module for more widespread use.

Our GBD module was built in the Poke^55,56 Python package currently available on Github. Poke was originally developed to be a polarization ray tracing module to study the influence of polarization aberrations on astronomical coronagraphs.⁴ For this study, we expanded its capabilities to include GBD. Poke operates using ray tracer API’s to trace a raybundle through every surface in the optical system. The relevant ray data are stored in a Rayfront object and can be loaded into a Python environment and interacted with independent of the ray tracer that generated the data (Fig. 4). The Rayfronts can also be compiled into binary file types using the msgpack⁵⁷ package and distributed to any interested investigator, effectively open-sourcing the physical optics calculations done on ray data.

Fig. 4

Illustration of Poke’s use as an interface to open-source ray-based physical optics. Poke only requires an interface between a ray tracing engine (orange, left) to generate and save a Rayfront object, which includes all ray data necessary for the GBD calculation. Currently, Poke supports sequential systems in CODE V and Zemax, but we are working on adding support for open-source packages that have ray tracers, like ray-optics.⁵⁴ The GBD field that Poke computes can then be sent to any open-source diffraction modeling package (yellow, right) to complete the hybrid propagation model.

In this study, we conduct ray traces in Zemax, which are then saved as a Poke Rayfront object. Poke performs the GBD simulations using the saved ray data to generate the field at the focal plane of the telescope. These data are exported to a coronagraph model built using HCIPy.

In Sec. 2, we outline the mathematics of Gaussian beam propagation and differential ray tracing used to perform GBD simulations. In Sec. 3, we present the mathematical algorithm for the transfer-matrix method of GBD that we developed for this investigation. In Sec. 4, we compare the results of observatory PSFs produced by GBD with one produced using traditional diffraction methods and analyze the artifacts that remain in the field. In Sec. 5, we assess the suitability of GBD for high-contrast imaging models and establish a roadmap for our module’s continued development.

2.1.

Propagation of a Single Gaussian Beam

The equation for a single Gaussian beam is parameterized entirely by the complex beam parameter $Q(z)$:³

$Q{(z)}^{-1}$ is a convenient expression of the Gaussian beam because it fully encapsulates the information required to describe the transverse electric field of the beam as it propagates. The real part of $Q{(z)}^{-1}$ is related to the radius of curvature of the wavefront:

In the paraxial regime, $Q{(z)}^{-1}$ can be propagated using the ABCD ray transfer matrices of geometrical optics:⁴⁸

To account for system misalignments, $Q{(z)}^{-1}$ is a $2\times 2$ matrix $\mathbf{Q}{(z)}^{-1}$ that encodes the complex curvature in two orthogonal directions and how they couple into each other, allowing for the beamlet to be generally astigmatic:^53,58

This treatment allows for greater versatility in the beamlet propagation but requires that the elements of the ray transfer matrices are also $2\times 2$ matrices. The propagation formula shown in Eq. (5) has a similar matrix extension:

We can then express the propagated Gaussian beam with the following equation:

2.2.

Computing the Differential Ray Transfer Matrix

The ABCD ray transfer matrix is a useful and concise method for analyzing properties of ray paths along optical systems. In the regime of geometrical optics, a generally skew ray can be traced through a system using $4\times 4$ ABCD ray transfer matrices.⁵⁹ These matrices model simple optical elements (e.g., thin lenses) with ease by operating on an input column vector that represents a light ray. The simplest ray transfer matrix that describes a paraxial and orthogonal optical system is a $2\times 2$ operator that maps an input ($i$) spatial and angular coordinate to the appropriate output ($o$):

For simplicity, it is convenient to represent the radial position in the plane transverse to propagation $(x,y)$ and the corresponding direction in the dimension $(l,m)$ as a position and angle vector, respectively $(\mathbf{r},\mathit{\theta })$. The ABCD matrix can similarly be condensed into $2\times 2$ submatrices that operate on each spatial dimension, yielding a familiar notation:

This description is powerful because it communicates the elegance and simplicity of ray transfer matrices. All dimensions transverse to propagation are accounted for, but the calculus to propagate a ray is still the same. The ray transfer matrices for simple and paraxial optical elements (e.g., thin lens) are well known,⁵⁹ and were used in concert with GBD in a prior investigation with paraxial systems.⁵³ However, our aim is to use GBD to model non-paraxial optical system diffraction. Therefore, we need a method of computing the ray transfer matrix for an arbitrary skew ray path. A simple dimensional analysis of the matrix relation in Eq. (11) is a good place to start understanding how to construct an arbitrary ABCD matrix. Because the position element $\mathbf{r}$ of the ray vector must be in units of distance, and the angular element $\mathit{\theta }$ must be dimensionless, the units of the ABCD matrix elements are constrained. $\mathbf{A}$ and $\mathbf{D}$ must be dimensionless and transform the ray position and angle through the optical system, indicating that they represent magnification. $\mathbf{B}$ and $\mathbf{C}$ must have units of distance and inverse distance, respectively. $\mathbf{B}$ operates on an angle and is therefore a metric of propagation distance through an optical system given some ray angle. $\mathbf{C}$ operates on a position and is therefore an indicator of the amount of refraction a ray experiences given a position in the entrance pupil.

Stone and Forbes’⁶⁰ work in differential ray tracing for inhomogeneous media was instructive in terms of deriving a method to construct the ABCD matrix. They illustrate the construction of individual optical elements through ray derivatives in a generally $4\times 4$ matrix through derivatives of surface data. Using this method, the position and angular derivatives of an optical surface are taken and arranged in a matrix like in Eq. (11). The matrix product of the optical elements is then the final differential ray transfer (or ABCD) matrix. Their method is functional if the analytical expression of the optical elements are known but could be very computationally intensive if the optical system has many elements. Instead, we approximate the differential matrix by tracing additional rays and compute the finite difference of the ray coordinates and directions at the input and output of the optical system to approximate the derivative.

Our implementation of the differential ray tracing technique in Poke utilizes a user-specified ray tracer (e.g., Zemax OpticStudio) that propagates rays by computing Snell’s law at each surface. For GBD, the ray coordinates of interest are on the source plane (where the field decomposition is done, e.g., the entrance pupil) and the transversal plane. The transversal plane is the plane normal to the central ray of a Gaussian beam that intersects the point in which we wish to evaluate the field. Propagation to this plane is critical, because the solution to Gaussian beam propagation [Eq. (8)] is only valid between planes orthogonal to the propagation direction defined by the central ray. An element of the ray transfer matrix can be computed by determining the ratio of the differential ray data on the transversal plane to the differential ray data on the source plane. An example of computing the element $A_{yy}$ is given by the following equation:

Fig. 5

Diagram illustrating differential ray tracing in the 2D case for a simple lens system. The central ray (black) is propagated along with two parabasal rays with a differential addition in the $y$ direction (red) and another with a differential addition in the $l$ direction (blue). To determine the ABCD matrix, the ray data are computed on the transversal plane, which is normal to the central ray and intersects the point in which we wish to evaluate the field.

The ray transfer matrix for a non-orthogonal optical system has 16 unknowns [see ABCD matrix in Eq. (10)], and each ray yields 4 quantities. To solve for every element of the matrix, 4 linearly independent rays must be traced. The simplest ray set is geometrically orthogonal,⁴⁰ where copies of the central ray $(x,y,l,m)$ are modified by a differential quantity $(\delta )$ in each of the four ray coordinates, two in position ($\delta x$, $\delta y$) and two in slope ($\delta l$, $\delta m$). The differential ray set is given by the following equation:

With differential ray tracing, we are able to propagate a single Gaussian beam through an arbitrary optical system using the ABCD matrix. To perform GBD, we next need to understand how to decompose the field in the entrance pupil of the optical system.

2.3.

Entrance Pupil Spatial Decomposition

The final variable to constrain in GBD is how to appropriately decompose the field in the entrance pupil into a finite set of Gaussian beams. This problem is illustrated in 1D earlier in Fig. 2; however, for imaging systems, the decomposition is a 2D problem. The fundamental Gaussian mode does not represent a complete set,⁶¹ and therefore, the decomposition of the field is not unique. We must carefully consider how the beamlets are distributed in the entrance pupil for accurate diffraction calculations.

Various sampling schemes exist in the literature, with different strengths and weaknesses. The even Cartesian sampling scheme [shown in Fig. 6(a)] described by Harvey et al.⁴¹ is the most straightforward, where the beamlets lie evenly spaced along a Cartesian grid. The ray coordinates in the entrance pupil are then computed from an overlap factor (OF), which describes the overlap of the beamlets’ $1/e$ waist radii $w_o$:

Fig. 6

Illustration of the (a) even, (b) Fibonacci, and (c) polar sample schemes for decomposing the field at the entrance pupil of the optical system with approximately the same number of beamlets in each figure. Quantifying the ramifications of these sampling schemes is of paramount importance to accurate diffraction simulation.

The Fibonacci sampling scheme [shown in Fig. 6(b)] introduced by Worku and Gross places the beamlets along a Fibonacci spiral, which results in a more accurate decomposition for circular apertures.⁴⁵ The distribution of the beamlets is even along polar angles on the spiral. The polar distribution of the beamlets is given by a position $R$ and angle $\mathrm{\Theta }$:

3.1.

Propagating Rays to the Transversal Plane

A GBD simulation begins by running a ray trace though an optical system in the user’s preferred design code. GBD needs the ray data at the plane where the decomposition occurred (the source plane) and the plane where we choose to observe the field (the evaluation plane). We use the ray data in evaluation space to propagate the rays to the transversal plane, where Eq. (8) is valid.

We first need to derive the propagation distance ${\mathrm{\Delta }}_{\text{ray}}$ for the central and four differential rays. To do so, we find the intersection of the line defined by the ray we want to propagate ${\mathbf{k}}_{\text{ray},E}$ and the plane normal to the central ray of the Gaussian beam ${\mathbf{k}}_{\mathrm{cen},E}$, which intersects the evaluation point ${\mathbf{r}}_{\text{point},E}$. This plane is the transversal plane and is defined by the following equation:

The line along the ray is defined by the following equation:

The expression in Eq. (20) is written for the general case of propagating a ray from the evaluation plane to the transversal plane. It is applied to the central and differential rays by substituting ${\mathbf{r}}_{\text{ray},E}$ and ${\mathbf{k}}_{\text{ray},E}$ for the position and direction associated with the ray to propagate. We update the ray positons by invoking Eq. (19). The ray position on the transversal plane (${{\mathbf{r}}^{\prime }}_{\text{ray},E}$) is given in the following equation:

However, these coordinates are still expressed in the basis of the evaluation plane. We must next rotate the ray coordinates into the transversal plane coordinate system so that they are orthogonal to the propagation direction.

3.2.

Transformation to the Transversal Plane

In the evaluation plane coordinate system, we define the orthogonal basis vectors $\mathit{\alpha }$, $\mathit{\beta }$, and $\mathit{\gamma }$ (shown in black on Fig. 7 to the right of the evaluation plane), to be the directions along the $x$, $y$, and $z$ axes, respectively. Similarly, the transversal plane is defined by basis vectors ${\mathbf{l}}_E$, ${\mathbf{m}}_E$, ${\mathbf{k}}_{\mathrm{cen},E}$, (shown in purple on Fig. 7), which are the analogous $x$, $y$, and $z$ directions for the transversal plane. To compute these directions, we start by taking the cross product of the central ray ${\mathbf{k}}_{cen,E}$ with the surface normal of the evaluation plane ${\mathit{\eta }}_E$ (shown in black on the right of Fig. 7):

A feature of Eq. (22) is that the ${\mathbf{l}}_E$ is orthogonal to both the central ray and the surface normal. Similarly, to determine our final basis vector, which is mutually orthogonal to both the central ray and ${\mathbf{l}}_E$ basis vector, we compute their cross product:

These vectors form a complete basis that describe the coordinate system of the transversal plane but are expressed in the coordinate system of the evaluation plane. Therefore, we can construct an orthogonal transformation matrix that performs a rotation of basis from the evaluation plane to the transversal plane. This matrix is shown in the following equation:

The same must also be done for the field evaluation point of interest ${\mathbf{r}}_{\text{point},E}$ to evaluate the Gaussian field later, shown in the following equation:

Now that the ray data are expressed in the transversal plane coordinate system, we can use them to compute the differential ray transfer matrix used to propagate the Gaussian beamlet.

3.3.

Computing the Differential Ray Transfer Matrix

To compute the differential ray transfer matrix for a Gaussian beam, we require the five ray coordinates and directions derived in the previous section (${\mathbf{r}}_{\text{ray},T}$ and ${\mathbf{k}}_{\text{ray},T}$), as well as the data for the same rays on the source plane (${\mathbf{r}}_{\text{ray},S}$ and ${\mathbf{k}}_{\text{ray},S}$). The differential ray transfer matrix is calculated using the methods discussed in Sec. 2. Recall that an element of the ray transfer matrix is given by the difference of the central and differential ray coordinates at the source and evaluation plane [see Eq. (12)]. The resultant matrix is given by Eq. (28). The full matrix including the ray data used to compute each of the elements of Eq. (28) can be found in Appendix A:

The total ray transfer matrix can be organized into a tensor composed of $2\times 2$ submatrices shown in the left-hand side of Eq. (28). These are exactly the non-orthogonal representations of the ray transfer matrix discussed in Sec. 2 that we need to compute the propagated Gaussian field profile.

3.4.

Computing the Gaussian Field

The central ray serves as the coordinate origin for our Gaussian field evaluation on the transversal plane. With the position vectors effectively transformed in Eqs. (25) and (27), we can define our centered radial coordinate ${\mathbf{r}}_o$ as the simple difference of these two positions to center the coordinate system in the following equation:

Finally, we can call upon Eq. (8) to perform the Gaussian field evaluation, with some added phase factors to account for the free-space propagation of the Gaussian beamlet. The result is shown in the following equation:

The procedure outlined above is repeated for each beamlet and location on the evaluation plane. From a computational perspective, this procedure is somewhat daunting. Repeating this process means the computation complexity scales by the number of beamlets and number of points on the evaluation plane, which can quickly become prohibitive. Fortunately, the algorithm described is entirely vector and matrix operations, so it is simple to take advantage of broadcasted array operations in Python to vectorize the computation and parallel computing to accelerate it. In our implementation, the vector operations described above are broadcasted such that the field of a single beamlet at all points is computed simultaneously. This operation exists in a loop over the number of beamlets (for our preliminary efforts in computational acceleration, see Appendix B).

Note the subtlety that the propagation procedure does not depend on the elementary field of choice until the electric field is evaluated. The assumption of this propagation method is that the elementary field is paraxial about the area encompassed by the differential rays. In principle, as long as a field’s propagation formula is known analytically via the general Collins integral, it can be propagated with this method. Thus GBD is a special case of the method described in this section. Because of the history of modeling resonators with ray transfer matrices, we know the ABCD propagation laws for Laguerre–Gaussian⁶² and Hermite–Gaussian⁶³ beams, which form complete sets. This indicates that the beamlet decomposition method would be of higher accuracy when decomposing the field into modes of higher spatial order. The formula for flattened elliptical Gaussian beams is also known,⁶³ which are capable of mitigating the soft-edge effect imposed by the traditional beamlet decomposition.⁴⁷ This method also works with Worku and Gross’s half- and quarter-truncated Gaussian beamlets.⁴⁶ The algorithm described in this section is one of the key results of this study, because it can propagate the field of any known solution to the Collins integral to an arbitrary array of points. We build this algorithm into Poke to test its ability to mimic diffraction simulations with known results and perform a comparison against commercial software.

4.1.

Fiducial Coronagraph

Our goal in this study is to assess the feasibility of GBD to integrate ray models of observatories into physical optics models of coronagraphs accurately. To quantify this, we propagate the images produced by GBD and the analytical Airy function through a charge-2 VC.^64,65 The optical VC is the general case of the vector VC,⁶⁶ which has shown great promise for future missions to image Earthlike exoplanets.⁶⁷ These coronagraphs are excellent at rejecting low-order spatial modes while transmitting the remainder of the light. We expect that traditional GBD will have some difficulty in accurately modeling high-spatial frequency content, and that it will manifest in the focal plane of this fiducial coronagraph if the error is limiting. If not, then we can conclude that GBD is suitable for high-contrast imaging simulation.

The complex amplitude of the charge-$q$ VC FPM is given by

4.2.

Observatory PSF

First, we examine GBD’s ability to construct the Airy pattern. The Airy pattern represents the “ideal” diffraction-limited observatory image for a circular aperture. The analytical solution for this image is known and available in POPPY, so we have a point of comparison that is not limited by numerical simulation errors. Figures 9 and 10 show the PSF simulations generated by GBD using the even and Fibonacci sampling schemes. We also plot a comparison of the modulation transfer function (MTF) to better illustrate how the transfer of individual spatial frequencies is affected by GBD. An important parameter in these simulations was the degree in which the effective diameter of the entrance pupil was appropriately captured. In Figs. 2(c) and 2(d), we observe that GBD results in some energy spillover outside of the original aperture function, which would result in a PSF with an incorrect spatial extent. To mitigate this effect, we remove any beamlets within half a waist radius of the aperture boundary. This is a trade-off between too much energy outside of the aperture, which results in a smaller PSF and too little energy within the aperture, resulting in a larger PSF.

Fig. 9

Comparisons of the (a,b,e) PSF and (c,d,f) MTF with (a,c) even sampling and (b,d) the analytical Airy pattern. The PSFs are given in units of normalized intensity, such that the sum of the energy in the PSF is unity. (e,f) The fractional difference and (g) the azimuthally averaged radial profile. The MTF is plotted out to the cutoff frequency of the HST of around $25\mathrm{cy}/\mathrm{mm}$. The radial oscillation in the fractional difference PSF is indicative of the artifacts introduced by GBD. In the MTF, it is apparent that frequencies below $20\mathrm{cy}/\mathrm{mm}$ are well-maintained, but the higher spatial frequencies increase in fractional difference. The effect on the PSF is revealed in the radial profile, where the PSF appears to spread out at larger angular separations. The RMS difference of the PSF data is $2.3\times {10}^{-6}$.

Fig. 10

Comparisons of the (a,b,e) PSF and (c,d,f) MTF with (a,c) Fibonacci sampling for the same number of beamlets as Fig. 9(a) and (b,d) the analytical Airy pattern. The PSFs are given in units of normalized intensity, such that the sum of the energy in the PSF is unity. (e,f) The fractional difference and (g) the azimuthally averaged radial profile. The MTF is plotted out to the cutoff frequency of the HST of around $25\mathrm{cy}/\mathrm{mm}$. The presence of ripples in the fractional difference PSF is noticeably lower than Fig. 9, and the fractional difference MTF shows that frequencies out to $23\mathrm{cy}/\mathrm{mm}$ are well-maintained. The spreading present in Fig. 9 is less prevalent in the radial profile. The RMS difference of the PSF data is $1.6\times {10}^{-6}$, indicating that this distribution results in a more accurate simulation.

Figures 9 and 10 compare the same simulation of the Airy function where the only variable is how the entrance pupil was decomposed. In both the PSF and MTF dimensions, it is clear that the Fibonacci sampling is the superior decomposition method for this aperture. For more undersampled cases, there may be a tradeoff in accuracy. To understand this, in Fig. 11, we examine how well the analytical Airy function is reconstructed for the two sample schemes as a function of the number of beamlets.

Fig. 11

RMS difference between the sum-normalized GBD PSF and analytical Airy function for the Fibonacci (blue circles) and even (black squares) sample schemes. For the circular aperture, it is clear that the Fibonacci sample scheme is more accurate for every case given a circular aperture. However, the returns are diminishing as the number of beamlets increases.

The Fibonacci sampling clearly wins out in terms of performance given a fixed number of beamlets, which translates directly to computation time. This also means that using the Fibonacci sample scheme can yield a simulation of the same accuracy as the even sample scheme with fewer beamlets. By judicious choice in sample scheme, the computational complexity of GBD can be lessened or the simulation accuracy can be increased.

To further test GBD’s capability to model a more realistic observatory PSF, we add obscurations in the entrance pupil that corresponds to the secondary mirror and supporting spiders. We also tilt the secondary mirror by 0.05 deg to aberrate the beam. By doing so, we simultaneously test our algorithm’s ability to capture aberrations from optical misalignment and diffraction from structure in the aperture. The results of this simulation are shown in Fig. 12.

Fig. 12

Simulations of an aberrated PSF with secondary support structure in the entrance pupil of our HST model given in Table 1. (a) The GBD PSF and (b) Zemax Huygens PSF have very similar structure with (c) the dominant difference in structure being the features from the sharp-edge diffraction from the spiders. This simulation helps quantify GBD’s difficulty in modeling sharp-edge diffraction effects. The RMS difference of the PSF data is $7.5\times {10}^{-6}$.

The data in Fig. 12 were simulated in a comparison between our proposed algorithm and Zemax’s Huygens PSF analysis tool. The Huygens PSF is computed by propagating spherical waves along ray paths in Zemax to the image plane and coherently summing them. The Zemax documentation asserts that this method makes fewer assumptions than the FFT PSF, so it was chosen as the point of comparison for the aberrated PSF. In our GBD simulation, we vignette any Gaussian beam that has any of its differential rays vignetted. Our results show extremely similar structure in the PSF from both the aberration induced by the secondary mirror misalignment and the structure from the HST spiders. The fractional difference reveals that the first couple “rings” of the PSFs are almost identical, indicating that the low-order aberrations were sufficiently simulated by GBD. There is a larger fractional difference in the dimmer structures that result from the sharp edges of the secondary obscuration and spiders, which is a known challenge for GBD simulations. This result can be further improved by careful implementation of alternative beamlet profiles, such as Worku and Gross’s truncated beamlets,⁴⁶ or higher-order transverse electric field modes (e.g., Hermite and Laguerre–Gaussian). However, the result using fundamental Gaussian modes is encouraging since the RMS difference of the PSFs is within the same order of magnitude as the results shown in Figs. 9 and 10. Zemax is an industry standard optical modeling tool, but its propagation algorithms are closed-source. It is beyond the scope of this paper to validate the accuracy of Zemax’s PSF simulation tools, but it is encouraging that our results in Fig. 12 agree. Now that our algorithm’s ability to perform PSF simulations has been verified, we can analyze the degree to which it introduces artifacts in high-contrast imaging simulations.

4.3.

Coronagraph Response

In traditional VCs, the on-axis field from an unvignetted circular aperture should be entirely rejected. We expect from the observatory PSF simulations that GBD does not trace high-spatial frequency information well due to the soft edges and amplitude ripples introduced by the beamlet decomposition. Any meaningful errors from this step should pass through the coronagraph unperturbed. To formally assess GBD’s suitability for high-contrast imaging, we construct the PSF of the fiducial observatory with a circular aperture using GBD and compare it to a Fraunhofer model of the same system. Both PSFs are propagated with Fraunhofer diffraction through the VC and the coronagraph focal plane is compared to assess the presence, if any, of residual signal introduced by GBD. The parameters used in this simulation are given in Table 2, and the results shown in Fig. 13 are cropped to show the innermost $30\lambda /D$ to better display the structure near the inner working angle.

Table 2

Simulation parameters for the result shown in Fig. 13.

Fig. 13

Comparison of the coronagraph PSF from 3 to $30\lambda /D$ generated by a solely Fraunhofer diffraction model (a) and the proposed hybrid propagation scheme for a varying number of beamlets across the 2.4-m unobscured HST aperture with an OF of 2. (b) 200 beamlets, (c) 300 beamlets, (d) 400 beamlets, (e) 500 beamlets, and (f) 600 beamlets across the aperture. These data were generated by producing a PSF with the given propagation scheme and then propagating it through HCIPy’s VortexCoronagraph class with a topological charge of 2. To reflect the residual starlight’s astrophysical flux ratio “contrast,” the PSF data are normalized to the maximum of an off-axis PSF at $\approx 16\lambda /D$ that propagates through the VC. The mean contrast of the masked region ($C$) is shown on each image. The residuals from (b)–(f) the hybrid propagation are brighter than the equivalent Fraunhofer simulation and minimize at the 500 beamlet case. This is particularly apparent near the inner working angle, where we see a maximum signal of $\approx 1\times {10}^{-8}$ for the 200 beamlet case decrease to $\approx 5\times {10}^{-9}$ in the 500 beamlet case.

Figure 13 shows the residuals that propagate through to the final image plane, which arise from inaccuracies in the propagation model. The hybrid propagation model generally minimizes in average contrast until the 600 beamlet cases [Fig. 13(f)], where there is a slight increase. The residuals from traditional Fraunhofer diffraction leave some low-energy features near the core of the PSF, which can be improved with greater sampling but the rest of the field is largely below ${10}^{-11}$ contrast. The decreasing residual energy with the number of beamlets used in the hybrid propagation scheme is encouraging, but we appear to have found the point in which increasing the number of Gaussian beams no longer increases simulation accuracy. The algorithm outlined in Sec. 3 is carried out as-written with some computational acceleration done by taking advantage of Python’s ability to vectorize matrix operations and other Python packages to accelerate the exponential calculation, discussed in detail in Appendix B. We have not formally explored parallel processing packages on central processing units (CPUs) or graphical processing units (GPUs) to accelerate this computation, but expect that these could make higher-sampled simulations more feasible to minimize the artifacts that remain in GBD PSFs for high-contrast imaging simulations. Now that we have established an open-source platform for GBD, Worku’s modified GBD⁴⁶ could be developed to minimize the number of beamlets required to minimize the artifacts in the coronagraphic focal plane. Exploring higher order spatial modes and the astigmatic fundamental mode of Gaussian beams could also yield greater accuracy for the same or less computational complexity. The result in Fig. 13 using traditional GBD shows that we can reduce the residuals to below ${10}^{-9}$ contrast near the inner working angle. It is also worth noting that this result suggests that with sufficient sampling, GBD can presently be used to simulate the PSFs for systems with less stringent contrast floors. To minimize the residuals from the propagation technique, other decomposition methods must be explored. We have created a platform for the algorithm’s development in an open-source environment and will outline a road map for the most pressing optimizations that can be conducted by future investigations to improve GBD’s accuracy.

5.1.

Future Work

The simulations presented in this article, while accurate, necessitated highly sampled simulations that were computationally intensive. The longest simulation used as many as 360,000 Gaussian beamlets across 2,560,000 pixels, which creates a phase array that is $\sim 100$ TB in size. This simulation was done on an AMD Ryzen9 3950X processor and took $\approx 24\mathrm{h}$ to complete. For GBD to be a practical diffraction technique for high-contrast instrument sensitivity analysis, it would be ideal to generate a field of the same accuracy in a shorter amount of time. Preliminary efforts in accelerated computing are discussed in Appendix B, where we were able to reduce simulation runtime from taking $\approx 5\mathrm{h}$ to complete to $\approx 10\min$. We intend to submit a follow-up study that outlines a more memory-efficient algorithm that we can use to more rapidly explore the degrees of freedom available in GBD (e.g., number of beamlets, OF, and etc.) for optimal diffraction simulation.

Modified GBD is already known to increase simulation accuracy for fewer beamlets⁴⁶ and as such is the next natural step in the development of our GBD module. Worku and Gross’s work on truncated Gaussian beams showed high accuracy in reconstructing the field after 2D polygonal apertures and “squeezed” the half-truncated beamlets in the azimuthal direction to increase the accuracy of a field after a circular aperture. Hexagonal apertures are of particular interest to astronomy because of their use in telescopes, such as the W.M. Keck Observatory and James Webb Space Telescope. Understanding how the truncated beamlets are able to reconstruct the field after segmented apertures is another important step in developing GBD for high-contrast imaging simulations. It is also worth repeating the notion in Sec. 3 that nothing about the proposed algorithm requires the beamlets to be strictly Gaussian. Therefore, considering alternative beamlets to use in the field decomposition may increase the accuracy beyond what traditional GBD and modified GBD are capable of. GBD has previously demonstrated the ability to model plane-to-plane diffraction effects (e.g., the spot of Arago),⁴¹ and its ability to model diffraction from surface polishing errors should be investigated for a more comprehensive modeling pipeline.

The beamlet decomposition algorithm would also benefit from more consideration into leveraging how parallelizeable it is. The contribution from each beamlet at each pixel is computed independently. With thoughtful consideration to the structure of our code and wide library of parallel processing packages available to us in Python on CPU’s and GPUs,^68,69 the beamlet computation could be rapidly accelerated. Preliminary experiments on GPUs have shown a decrease in runtime by a factor of $\approx 16\times$. Several investigators have been exploring the Jax Python package for its ability to perform automatic differentiation and parallel computing in support of physical optics simulations.^70–72 Automatic differentiation could be extremely useful for future beamlet decomposition algorithms to improve the accuracy of the ABCD matrix computation. Parallelization and vectorization is the most natural path forward for our beamlet decomposition algorithm because of the independence of the beamlet operations. Jax makes this process simple with the “pmap” and “vmap” functions while also allowing for just-in-time compilation for accelerated computing.

The use of a beamlet decomposition algorithm completely integrates a diffraction model with a ray model of an optical system, resulting in a more physically complete modeling pipeline. Consequently, other ray-based analyses can be integrated directly into the diffraction model. Polarization ray tracing⁵ is a natural extension to beamlet decomposition simulations because it can trace the complex amplitude of individual beamlets for generally vectorial field propagation through optical systems. This capability was demonstrated by Worku and Gross⁴⁴ in the context of high-numerical aperture microscope objectives but has the potential to be a powerful simulation tool for astronomical telescopes, including the next generation giant segmented mirror telescopes (Thirty Meter Telescope, Extremely Large Telescope, and Giant Magellan Telescope)⁴ and the Astro2020-recommended IROUV space observatory, which may be sensitive to polarization aberrations.

5.2.

Open-Source Science and Engineering

This research was inspired by POPPY, a Python physical optics module originally developed to simulate the James Webb Space Telescope. Our goal is to expand the capabilities of POPPY by investigating new propagation physics. We are developing the GBD module in a self-contained package with interfaces to other popular diffraction codes (POPPY and HCIPy) to provide ray information to diffraction models. The code used to prototype the GBD method can be found in the Poke repository on GitHub.⁵⁵ This repository is intended to be purely experimental as we develop the propagation physics for other high-contrast imaging packages that have substantive support.