Method Comparison for Simulating Non-Gaussian Beams and Diffraction for Precision Interferometry

In the context of simulating precision laser interferometers, we use several examples to compare two wavefront decomposition methods—the Mode Expansion Method (MEM) and the Gaussian Beam Decomposition (GBD) method—for their precision and applicability. To assess the performance of these methods, we define different types of errors and study their properties. We specify how the two methods can be fairly compared and based on that, compare the quality of the MEM and GBD through several examples. Here, we test cases for which analytic results are available, i.e., non-clipped circular and general astigmatic Gaussian beams, as well as clipped circular Gaussian beams, in the near, far, and extremely far fields of millions of kilometers occurring in space-gravitational wave detectors. Additionally, we compare the methods for aberrated wavefronts and their interaction with optical components by testing reflections from differently curved mirrors. We find that both methods can generally be used for decomposing non-Gaussian beams. However, which method is more accurate depends on the optical system and simulation settings. In the given examples, the MEM more accurately describes non-clipped Gaussian beams, whereas for clipped Gaussian beams and the interaction with surfaces, the GBD is more precise.


Introduction
In classic optics textbooks, e.g.[1,2,3,4], diffraction is defined as the phenomenon occurring when a wave is obstructed while propagating.This is, for example, the case when a Gaussian beam is clipped by an aperture.The phenomenon has been known for centuries, and there are various methods for the computation of diffracted wavefronts and their propagation.The most classic approach is the evaluation of diffraction integrals such as the Fresnel-Kirchhoff diffraction formula or the Fraunhofer diffraction equation, for instance, described in [5].
The propagation of diffracted light closely relates to the propagation of arbitrary wavefronts for which there exists no analytic propagation equation.Such arbitrary wavefronts include all clipped and diffracting beams, but likewise, half or quarter Gaussian beams were suggested to be applied in the GBD to optimize the simulation of sharp beam edges after passing through a hard aperture with an arbitrary shape [35].
Despite the number of publications using these methods, the number of publications describing and comparing the MEM and GBD in detail is currently low.Optimal settings for the methods are, therefore, often unknown, and the limitations of these methods are unclear.
For the MEM, Borghi et al. [16] derived the optimal decomposition parameters using Laguerre-Gaussian (LG) modes for circular symmetric fields of 1 mm radius, and particularly for top-hat beams.Yan Rong et al. [36] extended Borghi's optimal rule to an arbitrary radius of aperture.Liu et al. [19] presented the optimal decomposed beam waist for plane waves clipped by an arbitrary radius of aperture with Hermite-Gaussian (HG) modes.
Regarding the GBD, the publication status seems fairly sparse.About 30 years after the proposal of the basic concept of decomposing arbitrary electric fields into fundamental Gaussian beams, in Graynolds' overview article [21], he provided a revisit on the GBD, described the method's history and development over time and gave detailed implementation steps for the field decomposition, tracing, and computation of the resulting field in the target plane.In [37], various examples of modeling complex optical phenomena with the GBD were demonstrated, including interference and diffraction.However, these papers do not include a discussion of parameter settings, i.e. what waist size or waist location should be chosen for the grid beams, what overlap the grid beams should have or what type of grid would be ideal.
As previously mentioned, recent publications have further advanced the development of the GBD method.However, they again do neither address the question of parameter settings, nor do they compare the performance of the GBD with the MEM.Within this paper, we specify our experience values for parameter settings in the GBD, when it is used for simulating simple cases such as non-clipped and clipped Gaussian beams, for which we have analytical results available to compare with.We then directly compare the performance of the MEM and the original GBD method as introduced by Greynolds.This comparison is performed for strongly varying propagation distances, ranging from the common case of a few millimeters in the very near field, up to millions of kilometers.With this exceptionally large propagation distance, we test the applicability of the methods for the context of space interferometry and particularly space gravitational wave detectors like LISA and Taiji, for which the properties of the electric field needs to be characterized for distances up to 3 million kilometers.Additionally, we qualitatively test the MEM and GBD in decomposing and propagating aberrated wavefronts, for which we do not have an analytic result to compare with.Finally, we test the propagation through an optical setup by reflecting the decomposed fields from curved mirrors with different curvatures.
It should be noted that the MEM and GBD are based on the paraxial approximation and are only suitable for specific scenarios.This includes cases where the Gaussian beam waist is significantly larger than the wavelength.The MEM and GBD underly the paraxial approximation, because higher order or fundamental Gaussian beams are derived based on the paraxial approximation, this extends to their propagation using the ray transfer matrix method as well.
As the examples used in this context meet the assumptions of the paraxial approximation, the analytical methods are considered reliable and used as references.
After this introduction, we start this paper in Section 2 with an introduction of the MEM and GBD, specify the free parameters in both decomposition methods and define different types of errors to evaluate the precision of electric field estimates obtained by each method.We discovered that while the MEM performs well in accurately resolving far field wavefronts even with small mode orders, it is not able to ideally resolve high-frequency spatial oscillations in the near-field.The MEM's accuracy naturally improves with higher mode orders, it also improves with the propagation distance, which is typically within the ranges of interest.Similarly, we observed that the GBD's typical settings are inadequate for accurately resolving high-frequency spatial oscillations in the near field, while comparably smoother far fields are resolved with higher accuracy.The GBD's accuracy improves with increasing grid sizes.
We then discuss in Section 3 what mode order and grid size we chose to achieve a fair comparison of the methods.Afterwards we directly compare the performance of the MEM and the original GBD method as introduced by Greynolds in Section 4. This comparison is performed firstly in free space in Section 4.1, for strongly varying propagation distances, ranging from the common case of a few millimeters in the very near field, up to millions of kilometers.With this exceptionally large propagation distance, we test the applicability of the methods for the context of space interferometry and particularly space gravitational wave detectors like LISA and Taiji, for which the properties of the electric field needs to be characterized for distances up to 3 million kilometers.Our demonstration indicated that both techniques can accurately resolve the field at extremely far distances without requiring re-decomposition at 3 an intermediate point.A comparison between the direct methods revealed that the MEM outperforms the GBD when decomposing and propagating non-clipped circular and general astigmatic Gaussian beams.On the other hand, for clipped circular Gaussian beams, the GBD exhibited greater precision.However, these results may be dependent on various factors, such as the software and its implementation of both methods, the computer used, the operating system, and compilers.While in space gravitational wave detectors, the wavefronts will be aberrated, for instance, due to the telescope, and it also propagates through complex optical setups, therefore, we qualitatively test the MEM and GBD in decomposing and propagating aberrated wavefronts in Section 4.2, for which we do not have an analytic result to compare with.Based on our results, it can be inferred that both methods are generally suitable for decomposing and propagating aberrated wavefronts.Afterwards, we test the propagation through an optical setup by reflecting the decomposed fields from curved mirrors with different curvatures in Section 4.3.Our findings verified that when propagating through an optical setup that involves interactions with surfaces, the GBD method outperforms the MEM.Finally, we summarize the work and findings throughout this paper in Section 5.
All the simulations described throughout this paper have been performed using the software library IfoCAD [38], which has both an implementation of the MEM and GBD.For all simulations performed for this paper, circular symmetric Gaussian beams were used for the decompositions.However, the IfoCAD algorithm allows these beams to become simple or general astigmatic during propagation through the local setup.For this, the methods described in [39] are being used.The propagation of all modes and Gaussian beams follows the methods discussed in [9].

Wavefront Decomposition Methods
Within this section, we introduce the MEM and GBD methods in detail and individually test their performance using a simple exemplary case.We start with the MEM in Section 2.1 and continue with the GBD in Section 2.2.Both subsections are structured the same way: In Section 2.1.1 and Section 2.2.1 we specify the mathematical properties and implementation of the methods.In Section 2.1.2and Section 2.2.2, we define three kinds of different error definitions for estimating the decomposition precision.In Section 2.1.3and Section 2.2.3, we discuss how the MEM and GBD parameters should be chosen for minimal error.This is known only for the MEM for typical setups, while we can state only our experience values for the GBD.Finally, we individually test in Section 2.1.4and Section 2.2.4 the performance of each method on examples.

MEM: method description
The MEM is a well-known method defined, for instance, in [4], which describes the decomposition of an arbitrary wavefront into higher-order Laguerre-Gaussian (LG) modes or Hermite-Gaussian (HG) modes.LG modes are radially symmetric and therefore defined in cylindrical coordinates, while HG modes are defined in rectangular coordinates due to their axial symmetry.A conversion between both types of modes is known and described in detail, for instance, in [40,41].Throughout this paper, we therefore focus fully on decompositions using HG modes which are defined as HG mn (x, y, z; w 0d ) = c mn w(z) where w 0d is the waist of the fundamental mode HG 00 which is used as a parameter for all HG modes, and the beam radius w(z), the radius of curvature R(z) and the Gouy phase ζ(z) have the same definitions as for fundamental Gaussian beams.The coefficients c mn are normalization constants: The function H m (•) is the m th Hermite polynomial given by [4] H m (x) = (−1) m e x 2 d m dx m e −x 2 . ( An important property of these HG modes is that they are orthonormal and complete and thereby form a basis [4]: HG * mn (x, y, z; w 0d )HG kl (x, y, z; w 0d ) dxdy = δ mk δ nl , where HG * mn (x, y; w 0d ) is the complex conjugated Hermite Gaussian mode, and w 0d is the waist of the fundamental mode HG 00 which is used as parameter for all HG modes.The Kronecker delta function δ mk equals 1 for m = k and equals 0 if m k.This implies that any wavefront E(x, y) can be decomposed into a superposition E ∞ (x, y) of these modes [4]: where E ∞ (x, y) is a mathematically exact representation of E(x, y): where ≡ indicates the functions are equivalent for every point (x, y).The complex coefficients a mn with |a mn | 2 = P mn are usually referred to as mode overlap and describe how much beam power P mn is stored in each mode.They can be calculated by the inner product [4] a mn = HG * mn (x, y; w 0d )E(x, y)dxdy .
In real computations, it is not possible to use either an infinite number of modes in the decomposition (Eq.( 5)), or an infinite overlap integral (Eq.( 7)) for the determination of the mode overlap.Replacing the infinite surface integral in Eq. ( 7) with a finite one is uncritical provided the surface is chosen to be sufficiently large because electric fields of interest are usually fading out towards higher radial distances.Therefore, by choosing appropriately large integration boundaries, the introduced error becomes negligible.However, the error made by working with a finite mode order N is often non-negligible.Consequently, the decomposed field is no longer an exact representation of the input field E(x, y), but only an approximation: Here, we refer to N as the maximum mode order of the MEM.In Eq. ( 8), we use a triangular summation of the modes by summing n only up to N − m rather than N [41].This ensures, that within any decomposition, all polynomials up to the given order N are considered, and no polynomial orders larger than N are included.Consequently, in any MEM with mode order N, there are (N + 1)(N + 2)/2 HG modes superimposed.For radially symmetric fields E(x, y), the overlap a mn is set to zero if either the index m or n is odd.This means that for any mode order N the number ν of modes used in the MEM is given by

Error definitions for the MEM
The finite mode order decomposition given in Eq. ( 8) is not exact and will, therefore, have an error.It can be learned that for any given mode order, the decomposition error of the MEM depends on the mode order and the waist size w 0d chosen for the modes in the decomposition.Using the norm the normalized mean squared error (NMSE) is given by HG 2 mn (x, y; w 0d )dxdy using ∥E(x, y)∥ 2 = P, where P is the power of the initial beam.For any input field E(x, y), this normalized mean squared error ε NMSE (N, w 0d ) depends solely on the mode order N and waist size w 0d chosen during the decomposition, and it has the property of being propagation distance independent (cf.Eq. ( 7), and [16]).
The NMSE is defined via infinite surface integrals, which are replaced by numerical integrals over finite surfaces in optical simulations.This means, that in simulations a discretized NMSE (DNMSE) For radial surfaces and assuming radially symmetric beams, this is given by and for non-radially symmetric beams on rectangular surfaces by Here, r = x 2 + y 2 denotes the radial distance, N R , N X , N Y are the numbers of sampling points, and ∆r, ∆x and ∆y are the step sizes in the different directions, such that the maximal distances are Due to the assumed radial symmetric beams, we substituted R 2 dxdy by 2πr∆r for Eq.(13).Only the numerator is discretized in the DNMSE because the NMSE is normalized by the initial beam's total power P, which is usually known.
The discretized NMSE is a numerical representation of the NMSE, and thereby propagation distance independent provided that enough sampling points are chosen and the radial distance R is sufficiently large.However, this implies that with non-ideal settings, such as too few sampling points or a too small radial range, the error is indeed propagation distance dependent, which we highlight in Eq. ( 13) by the explicitly stated z-dependency.
One disadvantage of all shown errors is that they give only an integrated information and no distribution over a plane.We, therefore, define a sampling point dependent error ε rel (x i , y i , z) which we name relative error: For radially symmetric beams, we sample along the x-axis by setting y = 0. Therefore, the reduced 1D version of Eq. ( 15) can be written as: Throughout this paper, we use these relative errors to visualize the performance of the MEM, as well as for a qualitative comparison of the MEM with the GBD.To quantify the resulting information and judge the total error in the finite surface of interest, we define the summed relative error, for rectangular target surfaces and no assumed symmetry, or circular symmetric beams on a circular surface This summed error definition is, for typical settings in optical simulations, fairly independent of the chosen number of sampling points since Eq. ( 17) and Eq. ( 18) are representations of a discretized integral and therefore represent the surface under by the given function.However, the summed relative error is, unfortunately, not normalized.We use all introduced error types throughout this paper to study the performance of the MEM and to compare it with the GBD.

MEM settings
When a wavefront is decomposed using the MEM, there are only two parameters that need to be chosen: the waist w 0d of the modes used in the decomposition, and the maximum mode order N.For any maximum mode order, the choice of w 0d directly affects the magnitude of the mode overlap a mn and hence the resulting error ε NMSE (N, w 0d ).One can then for instance choose the decomposition waist w 0d such that the mode overlap of a specific mode is maximal (e.g.used in [13]) or, alternatively, such that the error made in the decomposition is minimal.Throughout this paper, we use the latter criterion, minimizing the error ε NMSE (N, w 0d ), following the examples of [16,19,36].
Which waist is optimal for the decomposition depends on the properties of the initial wavefront E(x, y, z).For an arbitrary wavefront, the optimal decomposition waist is therefore unknown.However, for the common special case of circular symmetric wavefronts originating from clipping at an aperture of radius R a the optimal waist was found to be [16,36,19,42]: This choice results in the minimum NMSE ε NMSE (N, w 0d ) in Eq. ( 12) for a given mode order N [16,36,19,42].

Example: MEM performance for a clipped Gaussian beam
We will now demonstrate the performance of the MEM in an example, for which the electric field is analytically known.In this example we investigate a clipped Gaussian beam.We assume that a Gaussian beam impinges orthogonally and perfectly aligned to a circular aperture with radius R a = 0.5 mm.The waist of the incident Gaussian beam has a radius of w 0 = 2 mm which is located in the aperture plane.The resulting circular symmetric clipped Gaussian beam is decomposed by an MEM with varying mode orders: N = 10, 20, ...50.For every mode order, the waist size w 0d of the modes is calculated using Eq.(19).We compute the electric fields amplitude and phase in various propagation distances and compare the results with the numerical evaluation of an analytic formula developed by Campbell in [43] for clipped Gaussian beams in the Fresnel region.Analogously, we use the analytical method of Tanaka et al. cf.[44,Eq.(1)∼ Eq.( 6)] for the Fraunhofer region.The distinction between the Fresnel and Fraunhofer region, is done using the Fresnel number F given by where λ denotes the wavelength of the beam, d the propagation distance after the clipping aperture, and R a the radius of the aperture.The near field refers to propagation distances, which make F larger than 1.If the Fresnel number is smaller than 1, the beam has propagated to the far field.In this example, we use the propagation distances d of 5 mm, 20 mm, 100 mm, i.e.F is 46.9925, 11.7481, 2.3496 in the near field, and d = 1000 mm with F = 0.2350 in the far-field.The number of sampling points for requesting the complex electric field is set to 3001.For convenience, all parameters of this example are listed in Table 1.The amplitude, phase, and relative error distribution calculated via Eq.( 16) are shown in Fig. 1 for different propagation distances after the clipping aperture with different mode orders.The lateral ranges chosen for this figure are unusually large.The corresponding results for a smaller lateral range are, therefore, shown in Fig. 2. All introduced error types have been calculated for both choices of lateral ranges and listed in Table 2.
The large lateral ranges used for Fig. 1 were chosen such that they are large enough to make the DNMSE propagation distance independent.In this case, the deviation between the input beam power and the MEM beam power was less than 2%.The incident beam power P is calculated simply from the Gaussian beam power passing through the aperture radius R a : with R a = 0.5 mm, w 0 = 2 mm and P 0 being the full power of the Gaussian beam prior to clipping.The resulting normalized power (P/P 0 ) of the clipped Gaussian beam is 0.1175.The power of the MEM beams were computed by the numerical sum R i=0 2π |E(r i )| 2 r i ∆r, which is a numerical representation of the denominator of Eq. ( 13).This procedure resulted in a slight variation of the MEM beam power in the different propagation distances.The deviations between the input beam power and the MEM beam power are 1.15879%, 1.15949%, 1.17214% and 1.65811% for propagation distances 5 mm, 20 mm, 100 mm, 1000 mm, respectively.Ideally, the lateral range would be chosen from the spot size of the clipped beam in the various propagation distances.Yet, particularly for clipped and diffracted beams, there is not one uniquely defined spot size, but rather a number of different concurring options, which are often not analytically known.Although a detailed discussion on spot sizes of clipped beams is beyond the scope of this paper, we want to compare the beams spot size with the chosen lateral ranges.In the near field behind the aperture, i.e. with Fresnel number F ≫ 1, the spot size of the clipped beam is still roughly equal the aperture radius and so in our example (F = 46.9925) the spot size is approximately 0.5 mm.Therefore, in row one of Fig. 1, the lateral range we are showing is 1.5 mm, which is approximately 3 times the spot size.For the propagation distance 1000 mm, which is in the far-field, we can use [45,Eq.(8)] to estimate that the spot size of the clipped beam to be 0.8751 mm.This is consistent with the spot size of a Gaussian beam with 0.5 mm waist at propagation distance 1000 mm, which The analytical methods for the near and far field are Campbell [43] and Tanaka et al. [44], respectively.Lateral distances for each propagation distance are chosen large enough to cover all the power.For these large lateral ranges, the MEM is effectively failing, generating zero amplitude and phases from lateral ranges that are about 3 times the spot size of the highest mode in the decomposition.Campbell [43] and Tanaka et al. [44], respectively.One can see that the further the beam propagates, or the higher the mode order is, the better is the performance of the MEM.
Table 2.The MEM errors, including the NMSE, discretized NMSE and the summed relative error, defined in Eq. ( 12), Eq. ( 13) and Eq. ( 18) respectively, are calculated for increasing mode orders at different propagation distances.The NMSE and discretized NMSE, which are propagation distance independent, when the lateral ranges R are large enough, are numerically equivalent.For smaller lateral ranges R, the discretized NMSE is propagation distance dependent.The summed relative error decreases with increasing mode orders for any propagation distance, and for a given mode order it increases (but not consistently) with increasing propagation distance.The number of sampling points X is 3001. is 0.8419 mm, and thereby slightly smaller than the clipped beam, as expected.Yet, in our computation, 3 times 0.8751 mm did by far not fulfill the propagation distance independent DNMSE, such that we had to extend the lateral range to 180 mm instead.These large lateral distances with a constant MEM beam power result in the expected propagation distance independency of the discretized normalized mean square error ε DNMSE , as shown in the fourth column of Table 2, except of some minor variations.We can see from Fig. 1 that for propagation distances of 20 mm and larger, there exists a maximum lateral range, for which the phase is correctly approximated by the MEM (see the zero lines for larger lateral distance).This is due to the finite size of the modes used in the decomposition.For a propagation distance z, the spot size of the higher order modes along x and y is [46] with w(z) being the spot size of the fundamental mode HG 00 used in the decomposition, and the Rayleigh z r = πw 2 0d /λ.Using these equations, we can estimate the spot size of the highest mode used in the MEM.The MEM can then only resolve fields in the range of maximally three times this spot size of the highest used mode.We can show this on the example of a propagation distance of 1000 mm (lowest row in Fig. 1) and a mode order of 50.For this, we find w x,50 0 = w y,0 50 = 32 mm, resulting in a maximal resolvable range of approximately 96 mm, which fits precisely the observation in the phase graph.
It might be expected that a higher mode order automatically implies that a larger lateral range can be resolved.However, that is not necessarily the case, as can be seen in Fig. 1 for a propagation distance of 20 mm.Here, this inverts: the higher the mode order, the smaller the resolvable lateral range.This is a consequence of using Eq. ( 19) to compute the optimal waist size, which decreases with increasing mode order.
Outside the maximal resolvable lateral range, the MEM is failing and generates zero amplitudes and phases.Consequently, the relative error is approximately 1 outside the maximal resolvable lateral range.This is clearly visible in Fig. 1.However, the large lateral range is not a choice usually taken in simulations since the spot properties are barely visible in these lateral ranges.Instead, simulations are usually performed with smaller lateral ranges in the target plane, such as shown in Fig. 2. Here, the lateral ranges cover only 0.4, 0.2, 0.06 and 0.022 times the ranges shown in Fig. 1, respectively.For instance, for the propagation distance 5 mm, the lateral distance shown in Fig. 2 is 0.6 mm, compared to the calculated spot size of 0.5 mm.For the propagation distance of 1000 mm, the computed spot size is 0.8751 mm, in comparison of the 4 mm lateral distance shown in Fig. 2.These lateral changes were chosen simply for good visibility of the amplitude and phase profiles without any hard criterion.
From the first row of Fig. 2, one sees that the MEM with the given settings and mode orders up to N = 50 insufficiently resolves the high-frequency spatial oscillation in the very near field behind the aperture.However, the further the beam propagates, the better is the performance of the MEM, such that after 1000 mm (i.e. at F = 0.235, the wavefront is well represented even with a mode order of 10.So while the NMSE is propagation distance independent and therefore constant for any choice of N we see from the left and center columns of Fig. 2 how the precision of the MEM increases with increasing propagation distance.This means that the error radially transmits outwards and might therefore be in a radial distance of no interest in the application.This also shows that it is not always necessary to choose high mode orders, particularly for far field simulations.Instead, the mode order should be chosen as a compromise between different criteria.The primary criterion is the increasing computational effort with increasing mode order.A second criterion is that the evaluation of the sum of Hermite-Gaussian modes with high polynomial orders is a typical mathematical challenge, resulting in numerical errors for high mode orders.Finally, the optimal decomposed beam waist calculated according to [16] and Eq. ( 19) decreases with increasing mode order and needs to be sufficiently large to not violate the paraxial approximation.Consequently, the mode order should be chosen carefully under consideration of the intended precision and the costs and risks if the mode order is chosen too high.
We can now compare the different errors for the case of the large lateral ranges.In that case, the DNMSE (fourth column of Table 2) is propagation distance independent and deviates from the analytically computed NMSE (third column) only slightly, with a maximum deviation of 9.32% (at 1000 mm with mode order 50).It can be seen that both the NMSE and DNMSE decrease with increasing mode orders for all propagation distances, as expected.In contrast to the DNMSE, the summed relative error ε rel shown in column 6 is not propagation distance independent, but increases for any mode order with the propagation distance.However, for any propagation distance, the summed relative error decreases with the increasing mode order.For the smaller lateral ranges, the DNMSE (fifth column of Table 2) is propagation distance dependent, as indicated before, but decreases with increasing mode orders for any given propagation distance, just like the NMSE.It generally decreases also with increasing propagation distances for a given mode order, but not strictly monotonously.Similarly, the summed relative error shown in column 7 is also propagation distance dependent and decreases with increasing mode orders for any propagation distance.For any given mode order, it increases as the beam propagates, due to the increasing step sizes ∆r.
Concerning the various error definitions, we find that the DNMSE is not propagation distance independent in typical simulation scenarios because the lateral range is then chosen too small.The relative error we have introduced here is a useful quantity that allows qualitatively judging the performance of the MEM directly from a graph.It allows, for instance, to directly see in Fig. 2), that the accuracy of the MEM increases with increasing mode order.The same finding is also found by the DNMSE or NMSE, but only in numbers that cannot be visualized comparably.In cases, where the relative errors cannot be clearly distinguished from the graph, like e.g. in the first row of Fig. 2, the summed relative error can help quantify physical dependencies (like the performance change with mode order or propagation distance).
In conclusion, we find that all defined types of errors have their individual strengths and weaknesses, such that a comparison of the performance of the MEM with the different error types can be helpful.Concerning the MEM itself, we find that it is not ideally resolving the high spatial oscillations of the diffracted beam in the near field but describes the beam accurately in the far field even if only low mode orders are used.

GBD: method description
Similar to the MEM, the GBD is a wavefront decomposition method.However, it decomposes any wavefront into fundamental Gaussian beams on a grid, as illustrated in Fig. 3.There are two supported shapes for the decomposition grid at the moment: square or hexagonal.Both grid shapes are depicted in Fig. 3.The quantities defining the grid are the edge length L, called window size, and the number of fundamental Gaussian beams along each dimension g.The lattice constant d g , called grid distance, is defined as d g = L g .The images show the grid of fundamental Gaussian Illustration of the GBD method.Shown on the left is the square grid structure, shown on the right is the hexagonal grid structure, both for even and odd g, respectively.One can see the (virtual) grid in blue, with the center marked in red.The origin points of the individual grid beams are marked in cyan, their waists with radius w 0g are shown as the dashed circles.The waist scaling factor f ws was set to 1.2 in the diagrams.
beams, depicted here by dashed circles that denote their waist radius w 0g .The waist radius is defined as where f ws is the so called waist scaling factor.For f ws = 1, the waists exacly touch each other.For larger f ws , the overlap of the fundamental Gaussian beams increases, for smaller f ws , it decreases, as shown in grid beams along each dimension is denoted by g with g × g = G being the total number of grid beams placed within the window in Fig. 3.We refer to this total number of grid beams as grid size.The definition of the window size and waist scaling factor shown here is the same as those in the IfoCAD, however, other software could have different definitions.
The hexagonal grid can be directly constructed from the square grid without changing the number of points or underlying math.Therefore, the algorithm to compute the GBD can be left unchanged when switching between grid geometries.To construct the hexagonal grid, columns with an even index are shifted up by d g 4 relative to the square grid position, columns with an odd index are shifted down by the same amount to create a hexagonal point structure.A rescaling by the factor of √ 3 2 along the horizontal direction is required to create an equidistant separation between the nearest neighbours of points, thus forming equilateral triangles.The mapping can be described by the function where x ′ i j and y ′ i j are the coordinates of the i j-th grid point calculated from the coordinates of the square grid point x i j and y i j and x 0 is the x-coordinate of the grid center.However, this transformation shrinks the window-size in x-direction due to the rescaling, resulting in a new window-size of √ 3 Section 2.2.4.Therefore, the mathematical description below will focus on the square grid, because the basic theory remains the same for both grids.The goal of the GBD is, to represent the wavefront as a superposition of the fundamental Gaussian grid beams weighted by coefficients: where E is the continuous wavefront to be decomposed, E i j are the electric fields of the grid beams with unity intensity, (x i j,0 , y i j,0 ) are their origin points and b i j are the complex weighing coefficients.To determine the coefficients, the above equation is evaluated at a discrete set of sampling points (x k , y l ), where k and l describe the location in the grid: one is the column, the other is the row.The resulting linear equation system is then solved for b i j .There must be at least as many sampling points as there are coefficients, which is G, one sampling point per grid beam.For better precision, you can also choose more sampling points.But because only the minimal number of required sampling points was used in the paper's simulations, we'll focus our explanations on this case.Both pairs of indices are compressed into a single sequential index to be able to write the linear equation system in matrix form.The E(x k , y l ) and the b i j can be written as column vectors, ⃗ W s and ⃗ b respectively.The E i j (x k − x i j,0 , y l − y i j,0 ) := m i jkl can be interpreted as a matrix M with each row corresponding to a sampling point containing the electric fields of each grid beam at this sampling point.These matrix entries describe, how strongly each grid beam influences the value of the superimposed electric field in the sampling point.Therefore, the GBD can be expressed as which can be solved by well understood methods such as the QR decomposition, which is adopted in IfoCAD.The following equations show in detail, how ⃗ W s , ⃗ b and M are composed for an equal number of grid beams and sampling points of G.
The GBD can be computationally expensive, if large grid sizes are chosen.For example, if G = 1000 × 1000, 10 6 beams will be superimposed, so both ⃗ W s and ⃗ b have 10 6 entries, which makes M of size 10 6 × 10 6 .Another approximation is employed to reduce the complexity of the problem further.Gaussian beam intensities drop off rapidly with increasing distances from the center.At points a few waist sizes apart, their contribution is near zero.Therefore, if the distance between the sampling point and the grid beam origin is larger then 3w 0g , the beam's contribution to the electric field at the sampling point is neglegible.The corresponding element in M can be set to zero.Consequently, M becomes a sparse matrix, and a software implementation making use of this can reduce both memory consumption and computational effort.Nonetheless, the high dimensionality of the equations should be kept in mind, when choosing the grid size of a GBD.The mathematical form of the sparsification can be expressed by where d i jkl is the distance between grid beam origin (x i j,0 , y i j,0 ) and sampling point (x k , y l ).

Finite GBDs and their error
To judge the quality and performance of the GBD to define a comparable set of errors like for the MEM.A NMSE for the GBD equivalent to Eq. ( 11) could be defined but not analytically evaluated like in the MEM (Eq.( 12)).Whether this error would be propagation distance independent is, therefore, not clear.However, we can still define the discretized NMSE comparable with the MEM error defined in Eq. ( 13) and Eq. ( 14) as follows Similarly, we define the 2D and 1D version of the relative error, and the summed relative error for the GBD according to Eqs. ( 15) to ( 18): We have now defined the same type of errors for the GBD as for the MEM and can use these for comparison.However, we do not know any major characteristics of the given errors when applied to a GBD.We, therefore, study and discuss their characteristics on the examples given throughout this paper.

GBD settings
For the MEM, it is known that for the stated set of applications, the relation defined in Eq. ( 19) can be used to achieve minimal error in the decomposition.For the GBD, we could not find any comparable information.It is, therefore, not clear how the parameters of the GBD should be chosen for minimal error.We can therefore only state fairly general information and the typical settings we choose in our simulations.
For the current implementation of the GBD in IfoCAD, the parameters which can be chosen explicitly, are the waist scaling factor f ws , the number g of grid beams along each primary axis of the square grid, and the window size L. The waist radius w 0g of the grid beams and grid distance d g are then determined using Eq.(24).Therefore, there are three parameters that influence the precision of the decomposition, from which the one dimensional number g of grid beams roughly compares with the mode order N of the MEM.One intuitively expects for a fixed windows size: the larger the number of grid beams g, the higher the precision, though this is valid only within a certain range.We show this property in Section 2.2.4 below.Furthermore, we investigate in Section 3 how to choose the grid length g and mode order N if both methods are being directly compared.We, therefore, understand the number g of grid beams as the primary handle for the precision of the GBD.The remaining two parameters (the waist scaling factor f ws and the window size L) are secondary handles, and we discuss their settings below.
Within this paper, we have two different types of examples: either a non-clipped Gaussian beam is being decomposed, or a wavefront that is clipped by an aperture.In the first example, the window size needs to be at least three times larger than the waist size, or else the beam would be clipped by the window during the decomposition, resulting in unintended and unphysical diffraction.On the other hand, the window should not be chosen too large to avoid an unnecessarily high number of grid beams with zero amplitudes and no influence on the final result.Comparable arguments hold for the second type of examples.Here, the window size needs to be sufficiently larger than the aperture.If the window was chosen to be smaller than the aperture size, the beam would obviously be clipped by the window, not the aperture.If the window size was chosen to equal the aperture size or be only slightly larger than it, the GBD would not be able to resolve the step function in the electric field, resulting in a GBD-beam with a considerable residual electric field amplitude outside the window.Only if the window is sufficiently oversized (compared to the aperture) the grid beams can resolve the step function in the electric field amplitude that originates from the clipping aperture and thereby accurately decompose the entire wavefront of interest.Within the examples of this paper, the window size is chosen between 1 to 1.5 times the diameter of the aperture.The window size equal to the diameter of the aperture is only used for the examples of non-clipped Gaussian beams (see Section 4.1), as in this case, the diameter of the aperture is already sufficiently large and does not clip the beam.Please note, the window size is defined as a full width, comparing rather to the diameter of the aperture.
The waist scaling factor should be chosen such that the grid beams have a non-negligible overlap.If the waist scaling factor was chosen too small, the GBD could accurately resolve the incident electric field in the grid points, but due to the lacking overlap of the grid beams, the GBD-beam would effectively have 'holes' between the grid points.On the other hand, if the waist scaling factor is chosen very large, a high number of grid beams contribute to the electric field in every sampling point, thereby significantly increasing the computational effort.Within this paper, the range of f ws chosen in all examples is between 3/2 and 10/3.
Unfortunately, we do not know an analytic relation between grid size, window size and waist scaling factor that forms an ideal choice for typical decompositions.The waist scaling factor and window sizes chosen within this paper are also not strictly optimized for the given examples, but simply follow the given logic.

Example: GBD performance for a clipped Gaussian beam
In this subsection, we show two examples on the performance of the GBD.In both cases, we decompose a Gaussian, which is clipped by a circular aperture.We assume normal incidence and the Gaussian beam to be optimally centred on the aperture.In the first example, we investigate the performance of the GBD with increasing grid size and compare and test the different error definitions.The second example illustrates the behavior of square and hexagonal grid shapes for the same grid size.
Example 1: comparing different grid sizes For the MEM, it is known from analytic equations that the precision of the decomposition increases monotonously with increasing mode order N.This is then also observed in simulations, as long as numerical errors are sufficiently small.For the GBD, one might likewise want to assume that for a fixed window size, the precision of the GBD increases with an increasing number of grid beams.However, we showed in Section 2.2.1 and Eq. ( 24) that the waist w 0g of the grid beams scales inversely with the number g of grid beams.This means the more grid beams are used, the smaller the grid beams' waist will get, provided the waist scaling factor is not adapted.Therefore, an increasing number of grid beams can quickly result in a violation of the paraxial approximation.It can, therefore, not be generally expected that an increasing number of grid beams is expected to increase the precision of the decomposition.In this example, we test the hypothesis of an increase in precision with an increasing number of grid beams.Intentionally, we work with a fixed window size L and a fixed waist scaling factor f ws and increase the grid size G up to values that cause the waist sizes to be in the order of the wavelength, thereby violating the paraxial approximation assumption.With this, we test in one simple example how the precision changes with the grid size, and we test slightly beyond settings that would normally be chosen.
The parameter settings of this example are listed in Table 3.In this example, the beam parameter, aperture size, shape and alignment, the propagation distances and sampling points are all chosen to be the same as in the MEM example in Section 2.1.4.The grid sizes are 100 × 100, 200 × 200, 500 × 500 and 1000 × 1000 respectively, using a square grid with a window size of 1.5 mm and a waist scaling factor f ws = 1.5 and the grid beam waist w 0g are calculated by Eq. (24).As shown in Table 3, the resulting waist sizes are critically small, up to a clear violation of the paraxial approximation in the case of 1000 × 1000 grid beams.
The amplitude, phase, and relative error are plotted for a large and small lateral range in Fig. 5 and Fig. 6, respectively.The corresponding errors are summarized in Table 4. Since the electric field of interest is the very same as in the MEM example, we use the very same lateral ranges as in Section 2.1.4.Table 4.The GBD errors, including the discretized NMSE and the summed relative error, defined in Eq. (31) and Eq. ( 36) respectively, are calculated for increasing grid sizes at different propagation distances.The discretized NMSE for both lateral ranges are propagation distance dependent.The summed relative errors for smaller lateral ranges decrease with increasing grid size at any propagation distance.The number of sampling points X is 3001.
propagation distance (mm) grid size DNMSE (for larger ranges) the summed relative error (for larger ranges) the summed relative error (for smaller ranges)    [43] and Tanaka et al. [44], respectively.Lateral distances for each propagation distance are chosen large enough to cover all the power.For any propagation distance, the performance of the GBD becomes better with the increasing grid size.5, but for smaller lateral ranges.Shown is the performance of a GBD with different grid sizes for an incoming circular Gaussian beam with a 2 mm waist being clipped by a 0.5 mm radius circular aperture.Shown are the amplitude (absolute value), phase and relative error ε rel (N R , R, z) at different propagation distances z = 5 mm, 20 mm, 100 mm (near field) and 1000 mm (far field) after the clipping aperture.The analytical methods for the near and far field are Campbell [43] and Tanaka et al. [44], respectively.One sees that the further the beam propagates, or the higher the grid size is, the better is the performance of the GBD.very descriptive because the lateral ranges are simply too large to judge the quality of the decomposition.However, Fig. 6 shows that the GBD describes the beam well, particularly in the far field, even if only 100 × 100 grid beams are used.Additionally, it can be seen that the further the beam propagates, or the higher the grid size, the better the performance of the GBD.This can also be seen from both the DNMSE (in column 4) and the summed relative error (in column 6), as they decrease with increasing grid sizes, although it is not strictly monotonous for increasing propagation distances with a given grid size in column 6.This is particularly interesting here, given that the high grid sizes imply unadvisably small waist sizes were used.The observed increasing precision of the GBD with these high grid sizes was, therefore, not naturally given.Additionally, the summed relative error in column 6 also changes nonmonotonically as the diffracted beam propagates.
From column 3 in Table 4, it can be seen that the DNMSE of GBD is propagation distance dependent not only for small lateral ranges but also in the case of large lateral ranges.This is unlike in the case of the MEM.Additionally, it can be seen that the DNMSE decrease with increasing grid sizes for any given propagation distance.This holds again for both choices of lateral ranges.However, it is not consistently given that for any choice of grid size, there is a strictly monotonous decrease of the DNMSE with increasing propagation distance.This is again different from the behavior of the MEM.However, it is currently unclear whether this originates from the method itself or its implementation.Likewise, also the summed relative error ε rel does not show a strictly monotonous decrease with increasing grid size.For example, the summed relative error in column 5 at a propagation distance of 5 mm increases in the step from a grid size of 100×100 to 200×200, and decreases for any further increase of the grid size.Finally, the summed relative error in column 5 is increasing with the propagation distance, unlike the DNMSE, which was mostly decreasing with the propagation distance.The reason for these observations could again originate from the method itself or be numerical precision or implementation problems, or the nearly zero-valued denominator in the error computation.Despite the non-monotonous behavior of the DNMSE and the summed relative error, we still use both for the total performance evaluation of the GBD in order to allow a direct comparison with the MEM.
Example 2: comparing grid shapes The shape of the grid can affect the accuracy of the GBD.We, therefore, repeat the previous example with the very same settings and a grid size of 500 × 500 beams, but compare this time the performance of the GBD with a square and a hexagonal grid.As introduced in Section 2.2.1, the window size in the horizontal direction is rescaled by a factor of √ 3/2 for the hexagonal grid, which is 1.5 mm× √ 3/2 in this example.We use the same waist scaling factor f ws = 1.5 and the waist radius w 0g for both square and hexagonal grid in IfoCAD.The resulting amplitude, phase and relative error of the square and hexagonal are plotted in Fig. 7.The corresponding DNMSE and the summed relative error are listed in Table 5.
Table 5.The DNMSE and the summed relative error for different grid shapes at different propagation distances given grid size 500 × 500.It can be seen that the hexagonal grid caused smaller errors for short propagation distances of 5 mm and 20 mm, while the square grid generated more precise results at higher propagation distances.The number of sampling points X is 3001.Both Fig. 7 and Table 5 show that at propagation distances i.e. 5 mm and 20 mm, the simulations performed with a hexagonal grid gradually show slightly better results than the simulations using a square grid.However, at larger distances of 100 mm or 1000 mm, the square grid resulted in higher accuracy.However, this is only one example, and we cannot draw generalized conclusion from it.  .Performace of a GBD with hexagonal and square grid shapes.Shown are the amplitude (absolute value), phase, and relative error distributions at different propagation distances from a circular aperture with 0.5 mm radius.The incoming circular-symmetric Gaussian beam was centered to the aperture and had its 2 mm waist located in the aperture plane.The analytical methods for the near and far field are Campbell [43] and Tanaka et al. [44], respectively.The number of sampling points X is 3001.

Fair Comparison
So far, the MEM and GBD have been introduced, their settings discussed, and their performance was individually tested on an example.In the next step, we want to directly compare the two methods, for which we need to define criteria to evaluate which method performed better or whether they performed equally well.Particularly, we need to define which mode order of the MEM should be compared with what grid size of the GBD and why this is chosen.This is discussed in the following.

Criteria for a fair comparison
There are two aspects which should be considered in the comparison of the MEM and GBD: accuracy and computational efforts.We judge the accuracy by the introduced errors, and the computational effort by the run time of the simulation.We consider a fair comparison a case, where either both methods achieve the same accuracy, and then the run time is used to judge the performance, or, both methods are set to have approximately the same run time, and the performance is judged by the achieved accuracy.Within this study, we use the later criterion and choose as run time the elapsed real time (not the CPU time).
From the basic principles of the two methods,the computational effort comprises three parts: the decomposition time, propagation time, and the superposition time.For the MEM, the decomposition time is the time spent on the integration to calculate the coefficients of the higher-order modes.For the GBD, it is the time span needed for the QR decomposition solving Eq. ( 27).Therefore, the decomposition time depends on the mode order and grid size, respectively, but also on the properties of the input field.The propagation time is very short in comparison because ray tracing methods, including the propagation of the Gaussian beam parameters with the ABCD matrix formalism, are highly efficient and computationally low demanding.For the MEM, all modes even share the same axis, such that only one ray needs to be traced, which then represents the beam axis of all modes.Finally, the superposition time depends on the number of sampling points in the target plane, as well as the number of modes or grid beams, respectively.The computational effort is therefore dominated by the decomposition time and the superposition time, and is naturally affected by other criteria such as the efficiency of the original implementation of the methods into the used software tool, the computational power of the used computer, and possibly even the used operation system.
In this study, we compare the performance of the MEM with GBD using IfoCAD (Version of 2022/10, git commit adf19a5b) to find mode orders and grid sizes which result in similar computational effort.All simulations shown here for testing the computational effort were performed on a MacBook Pro 2020 with 8 GB RAM and a 2.3 GHz processor with 8 cores.
Computational effort of the MEM and GBD In order to find settings that result in comparable computational effort for both methods, we performed a dedicated simulation where we varied the mode order and grid size for a fixed number of sampling points in the target plane.For this simulation, we chose again the case of a clipped Gaussian beam, with the same settings as in the previous examples: a circular Gaussian beam with waist 2 mm located in the aperture center incident onto this circular aperture with radius of 0.5 mm.We set N ranging from 10 to 100 with a step of 10 for the MEM and set g ranging from 100 to 1000 with a step of 100.The target plane was 5 mm away from the aperture and we used 3001 sampling points to compute the electric field for x ∈ −3, 3, y = 0. Here, the MEM and GBD were run with parallelization.We, therefore, distinguish two different times: the elapsed real time, which is the time the user needs to wait for a result, and the CPU time, which is the actual computation time and larger than the elapsed real time due to the used parallelization.We focus here on the elapsed real time and use this as the primary criterion.The computational effort of the MEM and GBD are summarized in Table 6 and Table 7 contain therefore information on the elapsed time during the decomposition, propagation, superposition.The total time is given both as elapsed real-time (column 5) and CPU time (column 6).
For the MEM we see in Table 6, that the decomposition time is by far the dominant time consumer in the given example, making out more than 95% of the total elapsed real-time.For the GBD, this is different: Table 7 shows that the superposition time is dominant over the decomposition time by a factor of more than 3.The propagation time is, as expected, insignificant.For the GBD with high grid sizes, some computational effort accumulates due to the large number of grid beams that need to be traced.For instance, for a grid size of 1000 × 1000, 1 million grid beams need to be propagated.Finally, the different superposition times are noteworthy.We understand this as a consequence of the different numbers of beams that need to be computed and evaluated at every target grid point.For instance, a mode order of 100 implies that 51 × 52/2 = 1326 modes are used in the MEM (cf.Eq. ( 9)).In the case of the GBD, a grid size of 500 × 500 grid beams implies that in every sampling point of the target plane, 250,000 beams are being tested, whether they contribute to the electric field.Even though only a few electric fields are indeed being superimposed in the end, the test itself for the high number of grid beams costs considerable time in the current implementation.
If we compare now the total elapsed real time for different settings of the MEM and GBD, we see several pairings of grid size G and mode order N, that can be used to achieve comparable computational effort.For instance, {N = 10, G = 100 × 100}, {N = 20, G = 200 × 200}, {N = 50, G = 400 × 400}, and we chose {N = 50, G = 400 × 400} for all comparisons using 3001 sampling points within this paper.
However, this choice depends on the number of sampling points used in the target plane since the computational effort of the GBD is dominated by the superposition time (which significantly depends on the number of sampling points in this plane), while the MEM is not.Therefore, the shown comparison should be repeated if a different number of sampling points is used.Within this paper, we use always the shown 3001 sampling points for all 2-dimensional cross sections of the electric field.However, we also show figures for the full cross-section (x and y for a fixed propagation distance z) and use 101 × 101 sampling points in these cases.Therefore, we repeat the above simulation for x ∈ −3, 3, y =∈ −3, 3 with 101 points each axis, i.e. 10201 sampling points total.The elapsed real time and the CPU time, are summarized in Table 8 and Table 9. We, therefore, choose {N = 50, G = 300 × 300} for simulations with 10201 sampling points, i.e. in cases where y is not set to 0 within this paper.Indeed, these parameters result only in roughly comparable computational effort, and better matching could be found if intermediate values were used.However, this is not necessary and is not the aim here, since the computational effort depends on additional simulation parameters and will change for other setups, particularly the properties of the wavefront that is to be decomposed, as well as the sampling grid in the target plane.However, our experience showed, that the given choice was consistently resulting in a comparable computational effort for all simulations compared within this paper.
Finally, we can compare the elapsed real time and the CPU time, and thereby the parallelization of both methods.For comparable parameters, i.e., {N = 50, G = 400 × 400} (Table 6, Table 7) and {N = 50, G = 300 × 300} (Table 8 and Table 9), the CPU time of GBD is longer.This shows that the GBD is more strongly parallelized than the MEM in the used IfoCAD version.

Method Comparison
In this section, we compare the performance of the MEM and GBD directly for various scenarios, which include non-clipped and clipped Gaussian beams in free space (Section 4.1), aberrated wavefronts (Section 4.2), and reflection from optical components (Section 4.3).

Non-clipped Gaussian beams and clipped Gaussian beams in free space
In every comparison within this subsection, we compute the introduced errors to evaluate the quality of each method.Unfortunately, this requires the electric field E to be known in every target plane, which strongly restricts the number of possible test cases.We, therefore, test in Section 4.1.1 the performance of the MEM and GBD for non-clipped circular and general astigmatic Gaussian beams for which the analytic representation of the electric field is widely known.In Section 4.1.2,we further investigate the case of circular symmetric clipped Gaussian Beams in the near-, far-, and extreme far-field, for which we can use again the analytic representations provided by Campbell [43] and Tanaka et al. [44] in the Fresnel and the Fraunhofer region.With extreme far-field, we refer to propagation distances of a few million kilometers, which occur in space gravitational wave detectors such as LISA [6] and Taiji [7].

Non-clipped Gaussian beams
Circular Gaussian beam The simplest case to compare the MEM and GBD is using non-clipped Gaussian beams, for which the electric field is analytically known in any propagation distance.We, therefore perform a first comparison of the MEM and GBD on the example of a non-clipped circular-symmetric Gaussian beam with the parameters listed in Table 10.It is noteworthy that we still define an aperture here.This is due to the IfoCAD version used here, which requires an aperture to be defined for both methods.However, with a radius being 4 times larger than the Gaussian waist radius, it is effectively not clipping the beam, given that the clipped power is approximately 1.3 × 10 −12 % of the full beam power.Instead, the aperture radius effectively defines the lateral range used in the numerical evaluation of the integral in Eq. ( 7).In the GBD, we are then using the very same aperture radius to decompose the exact same input field as in the MEM case.We chose not to overscale the window further and set the window size (which is a full width) to equal the aperture diameter, to not place unnecessarily many grid beams in regions without field amplitude.
Concerning the optical setup, we define a circular Gaussian beam with a waist radius of 1 mm which is centered in this aperture.The mode order, grid size, and sampling points are chosen according to Section 3. The waist w 0d of modes used in the MEM is calculated from Eq. ( 19); the GBD grid beam waist w 0g is calculated correspondingly by Eq. ( 24).After the decomposition, the MEM beam and GBD beam propagate for z r /1000, z r , 1000 z r and 3 million kilometers (3 Gm), with z r = 2.9526 m, where z r refers to the Rayleigh range of the incident Gaussian beam.We then speak of the far-field, if z ≫ z r [4].The resulting amplitude, phase, and relative error are shown in Fig. 8 10.The analytical results refer to the complex electric field of the non-clipped Gaussian beam in this case.It can be seen that MEM behaves better than GBD.
the corresponding discretized NMSE, and the summed relative errors are summarized in Table 11.The chosen lateral ranges of 3 times the local spot sizes cause all amplitudes shapes to appear identically.This also holds for the shape of the phase profile.However, this is not immediately visible due to phase wrapping particularly in the far-and extreme far-field, which can be resolved by using a phase-tracking algorithm.Figure 9 shows the unwrapped phase for the propagation distance of about 3 km, i.e. 1000 z r .From Fig. 8, it can be seen clearly that both MEM and GBD can accurately represent the circular Gaussian beam.However, in all shown propagation distances, the MEM is more accurate than the GBD (right hand side column of  graphs), which can also be seen from the discretized NMSE ε DNMSE (column 3) and summed relative error ε rel (column 4) in Table 11.The DNMSE error of the MEM is notably small, considering that the waist of the modes used in the MEM determined through Eq. ( 19) is 0.8 mm, which closely aligns with the waist of the non-clipped Gaussian beam.While it is possible to obtain a 0 error by choosing the waist of the modes equal to the non-clipped beam waist of 1 mm, such a comparison would be meaningless in this particular scenario.However, this agrees with the analytically calculated NMSE (Eq.( 12)) which is 1.3989×10 −14 .The discretized NMSE of the MEM result is found to be even slightly smaller than this, with a residual propagation distance dependency.Both of these properties originate from the finite radial range used in the decomposition and error computation.
In conclusion, we find that in this example both methods accurately resolve the incident wavefront, and that the MEM has exceptional precision and is, therefore, more accurate than the GBD.
A particular challenge in this simulation was the ultra-large propagation distance of 3 million kilometers.This gives natural rise to numerical precision problems and was a major concern and one of the initial reasons for why we tested this extreme far field, which is relevant for space-based gravitational wave detectors.The given example shows no apparent signs of numerical limitation.This was achieved by a separation of the optical pathlength (i.e. the ikz-term in the Gaussian beam), from the residual phase contributions [47].
Additionally, it was expected that the GBD could not propagate the beam into this extreme far-field without a re-decomposition in an intermediate plane.This expectation originated from the chosen small grid in the original decomposition plane, compared to the very large spot size in the target plane.After all, in the given example, the spot size of the total beam is 1 mm in the decomposition plane, for which we chose a window size of 8 mm.Since the beam is not re-decomposed, this is also being used in the target plane, where the Gaussian beam radius has increased to 1.0160×10 3 km (cf Fig. 8).Yet, despite that the grid is unfit for the target plane, it is clearly visible in Fig. 8 that the beam is well represented.A re-decomposition was not needed in the given example.
General astigmatic Gaussian beam After testing the performance in the case of a non-clipped circular symmetric Gaussian beam, we now test the performance of these methods on a general astigmatic Gaussian beam.In that case, two beam waists and a complex angular orientation θ need to be defined [48].The parameters used for the performance comparison are listed in Table 12.In this case, the lateral parameter y is not set to 0. As shown in Section 3, we therefore use a mode order of N = 50 and compare this with a GBD grid size of 300×300.The resulting electric field profiles of the MEM are plotted in Fig. 10 and for the GBD in Fig. 11 for three different propagation distances: z r1 /100, z r1 , 100z r1 , with z r1 = 2.9526 m being the XZ-plane Rayleigh range.We do not show the results of 3 Gm in this case because the sampling points of 101 × 101 is too low for such an extreme-far distance and we have already shown that both MEM and GBD are accurate at 3 Gm for the non-clipped circular Gaussian beam.We have computed the electric field in lateral distances (x, y) out to 2 times the spot size.This choice results in good visibility of the amplitude profiles, however, it slightly masks the magnitude of ellipticity in the resulting images.
In the first column of both figures, one sees the known characteristics of a general astigmatic Gaussian beam: its elliptical amplitude pattern which has a rotating orientation during propagation.Both methods equally well resolve this characteristic.The second column, showing the phase of the general astigmatic beam, is mostly smooth at propagation distances less than 3 m (top and central row), except of two lines of phase jumps.In the far field (lowest row), there is a pattern of circular shapes.This pattern is an aliasing effect originating from the high wavefront curvature, the resulting high number of phase jumps in combination with the low sampling rate.We demonstrate this effect and how it can be resolved for a cross-section of the phase profile in Fig. 12.Finally, the third column in Fig. 10 and Fig. 11 shows the accuracy of each method.The DNMSE and summed relative error of the MEM and GBD at different propagation distances are shown in Table 13.It can be seen from this table, as well as from the third column of both Fig. 10 and Fig. 11 that the MEM is again consistently more accurate than the GBD.Additionally, one can see that both the DNMSE and the summed relative error of the MEM increase with the propagation distances, while for the GBD, the DNMSE nonmonotonically 29  All the simulation parameters are listed in Table 12.It can be seen from the third column of this figure that MEM shows a good agreement with the initial general astigmatic Gaussian beam.12.It can be seen from the third column of this figure that GBD shows a good agreement with the initial general astigmatic Gaussian beam.The left-hand side image shows the phase using 101 sampling points, the centered image for 3001 sampling points, and the right-hand side image again for 3001 sampling points after a phase-tracking algorithm unwrapped the data.One sees that the small side-maxima, which are seen as circles in Fig. 10 and Fig. 11, disappear when a higher sampling rate is used, and indeed a smooth but strongly curved phase profile is restored in the far field.
decrease with the propagation distances, and the summed relative error again changes inconsistently.

Clipped Gaussian beam
In this subsection, we directly compare the performance of the MEM and GBD for a circular Gaussian beam clipped by a circular aperture.The parameter settings are summarized in Table 14, for convenience.In this case, the diffracted beam propagates 5 mm, and 100 mm in the near field, 1000 mm in the far field, and 3 Gm in the extreme far field, with Fresnel numbers of 46.9925, 2.3496, 0.2350 and 7.83208 × 10 −11 , respectively.The resulting amplitude, phase and relative error profiles are depicted in Fig. 13, and the corresponding errors are summarized in Table 15.
We have chosen again lateral ranges that allow good visibility of the amplitude profile.For 3 Gm, we provide two lateral ranges: approximately 2 times the spot size and a smaller lateral range of 400 m.In fact, we did not attempt to compute the exact spot size of the clipped Gaussian beam, but simply estimate it from two boundary cases: the spot size of a top hat and Gaussian beam.The spot size of a top hat beam with 0.5 mm radius can be estimated using [45, Eq.( 8)] resulting in 2625.3 km after 3 Gm.For a Gaussian beam with 0.5 mm waist, the spot size at a propagation distance of 3 Gm would be 2032.1 km, and the spot size of the clipped Gaussian beam is expected to be between these two values.We use here the spot size of the top-hat for the lateral ranges for simplicity.In Fig. 13, it can be seen that both methods describe the clipped Gaussian beam well, however, in all shown propagation distances, the GBD behaves better than the MEM.Quantitatively, we can also reach the same conclusion from the discretized NMSE and the summed relative errors listed in Table 15.Especially in the near field with propagation distances 5 mm and 100 mm, the differences in the errors of the MEM and GBD is significant.This is the same finding as in Section 2.1.4,that the MEM is insufficiently resolving the high-frequency spatial oscillation in the very near field.With the increase of the propagation distance, at 1000 mm (third row), the differences between both methods become smaller, and in the extreme far-field 3 Gm (fourth and lowest row), they narrow further.We can also 32  14.The analytical methods for the near and far field are Campbell [43] and Tanaka et al. [44], respectively.The graphs in the third column show the relative error, which indicates that the GBD performs better than the MEM.
see that both methods become more accurate with increasing propagation distances.Like before, the extreme far-field electric field was computed with the GBD in one step and did not require a re-decomposition in an intermediate plane.
The 0.5 mm aperture radius is a typical value in laboratory experiments and, therefore, fits well with the shown near-field propagation distances.However, it is not a realistic value for the extreme far-field simulation case originating from space gravitational-wave detectors.The aperture diameter in LISA-like missions is usually between 20 cm and 40 cm.Therefore, here we show another example, where both the beam waist diameter and aperture diameter are 30 cm, for propagation distance 3 Gm.Other parameters are chosen from Table 14, except the window size, which is 35 cm in this case, and the lateral range is chosen 3 times the spot size estimated by [45,Eq.(8)].The resulting amplitude, phase, and relative error are shown in Fig. 14  In this figure, it can be seen that the order of magnitude of the amplitude is 10 −8 (left-hand side image), and the relative error shown in the right-hand side image indicates that the GBD is again more accurate than the MEM.

Aberrated wavefronts
Aberration is a typical phenomenon in optics and a case in which the beam cannot be propagated analytically through the setup.In LISA-like missions, aberration occurs for instance when the beam propagates through the telescope.The beam that is launched towards the remote spacecraft is, therefore, not a perfect clipped Gaussian beam but has additional wavefront distortions.Likewise, the received beam is not a perfect top hat beam.These aberrations affect the readout noise, and therefore need to be studied.
We study here how the MEM and GBD decompose and propagate aberrated wavefronts.Unfortunately, we do not have an analytic solution to compare the results with.We can, therefore, only qualitatively compare the results of both methods without being directly able to judge which one is more accurate.Instead, we test whether the methods generate qualitatively agreeing results.
Mathematically, the wavefront aberration can be described by adding an additional phase term Ω a in the complex electric field of the original beam, i.e.: E a (r; 0) = E(r; 0) exp (ikΩ a ) , where E a is the complex electric field of the beam with aberration, E is the complex electric field of the original beam, and Ω a is the additional phase distribution caused by aberration.Such a phase term is usually described by Zernike polynomials [49]: where c m n represents coefficients and Z m n are Zernike Polynomals, where n − m ≥ 0 is an even number.To demonstrate that both the MEM and GBD can describe wavefront aberrations, we show their results for the well-known effects that are caused by individual Zernike polynomials up to fourth order (cf.[50,Fig. 4], [51, Fig. 3], and [52,Fig. 9.8]).The effect of optical aberration is often represented by calculating the point spread function (PSF) [50], i.e. the intensity profile at a distance z which is usually computed by Fourier transformation: which is also known as the response of the pupil function after Fourier transform (FT) under a certain distance z [53].
In the far field, the beam source can usually be regarded as a point source compared to the propagation distance, and the PSF equals then the intensity profile computed by Fraunhofer diffraction.In this section, we use the MEM and GBD to represent the effect of optical aberration, by calculating the amplitude profile of a diffracted Gaussian beam with aberration at the Fraunhofer region, instead of directly calculating the PSF.
In this example, we calculate amplitude profiles of the first 4 orders of Zernike polynomials by the MEM and GBD for an aperture with 1 mm radius at λ = 1064 nm, the coefficient of Zernike polynomials c m n for all m and n is set to 10. Parameter settings of the example are listed in Table 16.The propagation distance is 5 km (F = 1.8797 × 10 −4 ), therefore, it is in the Fraunhofer region.The MEM results are nearly identical and are not shown here to avoid unnecessary duplications.Instead, we show the difference between the MEM and GBD results in Fig. 16.It can be seen from this figure that the difference between MEM and the GBD is in the order of 10 −7 to 10 −6 .If this is related to the individual amplitudes of approximately 10 −5 to 10 −4 , we would speak of relative deviations in the order of a few percent.Even though we cannot judge which of the two methods is more accurate, we can conclude that either method can generally be used for decomposing and propagating aberrated wavefronts.Please note, we intentionally plot here only the difference between the results, and perform only a qualitative relative deviation, because the computation of a relative difference like |E GBD − E MEM |/|E MEM | causes divisions by zero with the shown lateral ranges.

Reflection from optical components
In the previous subsections, the behavior of the MEM and GBD in free space propagation have been compared.However, a major difference between both methods arises when the decomposed beams are propagated through an optical setup.This is an important test case because of the differences in the decomposition methods: while all MEM modes share the very same axis, the GBD grid beams are distributed on a grid.This means when an MEM beam interacts with a surface (i.e.reflects or refracts), it effectively interacts only with the intersection point and its local curvature, while the GBD beam probes the surface in multiple intersection points.We illustrate this difference for a simple test case: the reflection of a Gaussian beam from a spherically curved mirror with varying curvature.We then expect the GBD to show increasing levels of spherical aberration with increasing mirror curvature, while the MEM is expected not to resolve the occurring spherical aberration.This means we decompose a Gaussian beam with an MEM  and a GBD and reflect the original Gaussian as well as the MEM and GBD representations of this Gaussian from a spherically curved mirror.We then expect for increasing mirror curvature an increasing level of deviation between the GBD beam and the reference Gaussian, while the MEM is not expected to show this behavior.The results of this simple test are illustrated in Fig. 17 and for a stronger curvature case in Fig. 18, which confirm the expected behaviour.
For this simulation, we assumed the Gaussian beam with a 1 mm waist radius located in its origin which propagated by z = 10 mm before it impinged orthogonally and centered onto the spherically curved mirror which had a diameter of x mm.The mirror was, therefore, sufficiently oversized to reflect the full beam.After reflection, the complex electric field was calculated at an observation plane, which was 5 mm away from the mirror and orthogonal to the beam.All simulation parameters are summarized in Table 17.
Please note, we use the wording 'relative error' in Fig. 17 for consistency with the previous sections.However, this should be rather seen as a relative deviation, given that the Gaussian beam, which serves as a reference, cannot be trusted to be physically correct in this example because it is itself insensitive to spherical aberration.Table 17.Parameters list for reflection from optical components.Like in Section 4.1.1,we define here an aperture only for implementation reasons.However, the beam is effectively not being clipped.Therefore, the complex electric fields of the MEM and GBD in the observation plane can be directly compared with the results for a fundamental Gaussian beam.waist of the grid beam used in the GBD 0.0333 mm grid shape grid shape of the GBD square C curvatures of the mirror 0 mm −1 ,−0.002 mm −1 , −0.02 mm −1 , −0.1 mm −1 size the diameter of the mirror 1 cm×1 cm X number of sampling points 3001

Summary and Conclusion
In this paper, we have compared the two wavefront decomposition methods MEM and GBD for different test cases.To judge the performance of either of the methods and to allow a direct comparison of both, several different types of error estimates have been introduced: the normalized mean square error NMSE, its discrete analog DNMSE, a relative error, and its sum.The properties of all these errors were discussed and compared.
We found for the MEM that even though the well-known NMSE is propagation distance independent, the DNMSE is usually not so because the lateral ranges chosen in the target plane are too small.To achieve the propagation independency of the DNMSE, the lateral ranges in the given example were so large that the MEM with mode orders up to 50 was not able to resolve the necessary lateral range due to the finite spot sizes of the involved modes.
While the NMSE and its discretized analog are commonly known and used errors, they do not allow visualizing the error distribution over the cross-section of the field of interest.For this, the relative error can be used.Finally, the summed relative error is a useful addition to the relative error, quantifying the graphical findings.
To allow a direct comparison, we have tested the simulation runtime for different settings in a typical test case.We showed that the GBD runtime significantly depends on the number of sampling points in the target plane, unlike the MEM.A direct comparison of the precision of these methods for comparable runtime, therefore, depends significantly on the number of sampling points chosen in the target plane.Naturally, particularly this finding depends on the chosen software tool and the implementation of the methods.All simulations performed throughout this paper were performed using the software library IfoCAD (runtime tests with Version of 2022/10, git commit adf19a5b).However, the introduced method of finding settings that allow a fair comparison is independent of computer systems, software  17.The relative deviation of the GBD from the Gaussian beam and MEM increases when the curvature of the mirror increases.choices, and implementation details.Additionally, the findings underline that the GBD method in any software implementation should not only be optimized in the decomposition but also for the superposition in the target plane.
For the individual performance of the MEM, we found that it is not ideally resolving the high-frequency spatial oscillations in the near-field, but resolves the far field wavefronts accurately even with small mode orders.The accuracy of the MEM naturally improves with higher mode orders, but in the usual ranges of interest also with the propagation distance.
For the individual performance of the GBD, we likewise found that the high-frequency spatial oscillations in the near field are insufficiently resolved with typical settings.The comparably smooth far fields, in comparison, are resolved with higher accuracy.Naturally, the accuracy of the GBD becomes better with increasing grid sizes.However, this can quickly result in grid beam waists that are so small that they violate the paraxial approximation.We intentionally tested the performance of the GBD in such an imperfect case and found that the precision was not impaired by the non-ideal, extremely small waist sizes.That means the relative error and its sum decreased for GBDs with increasing grid sizes, despite the use of smaller and smaller waist sizes that violated the paraxial approximation.
We have directly compared the MEM and GBD for cases where the electric field amplitude and phase are analytically known in different propagation distances.We have performed this test for typical propagation distances in the near and far-field and additionally in the extreme far-field of millions of kilometers, as needed for space gravitational wave detectors.We showed that both methods can resolve the field at this extreme far distance without the need for a re-decomposition at an intermediate distance.The direct method comparison showed a better performance of the MEM for the decomposition and free-beam propagation of non-clipped circular and general astigmatic Gaussian beams.In the cases of clipped circular Gaussian beams, the GBD showed higher accuracy.However, these findings might well depend on the settings, the used software with its implementation of both methods, and the computer, operation system, and compilers.
Additionally, we compared the MEM and GBD representation of aberration and showed a qualitative agreement between the results.Finally, we showed in one example that the GBD is a superior method for propagation through an optical setup, where interactions with surfaces occur.While the MEM decomposes the initial field with modes that all share the very same beam axis, the GBD decomposes into fundamental Gaussian beams on a grid.Consequently, the MEM beam probes the curvature of a surface only at one intersection point between its beam axis and the surface, while the grid beams of a GBD beam probe the surface curvature in a grid of intersection points.We have shown this difference by a qualitative comparison of the electric fields of a Gaussian beam and its MEM and GBD representations after reflection from a curved mirror.We have shown that the GBD beam showed the expected spherical aberration unlike the Gaussian or MEM beam.
We can generally conclude that both methods are useful for decomposing and propagating non-Gaussian beams.Once the fields are decomposed, the propagation in free space or through an optical setup is computationally trivial.Which method is more accurate depends on the test case and simulation settings.However, for the propagation through optical layouts, where the beam interacts with surfaces, and particularly if non-spherical surfaces exist in the setup, the GBD with its grid of beams has a clear advantage over the MEM.

Figure 1 .
Figure1.Performance of a MEM with different maximum mode orders N = 10, 20, 30, 40 and 50, for an incoming circular Gaussian beam with a 2 mm waist being clipped by a 0.5 mm radius circular aperture.Shown are the amplitude (absolute value), phase and relative error ε rel (N R , R, z) at different propagation distances z = 5 mm, 20 mm, 100 mm (near field) and 1000 mm (far field) after the clipping aperture.The analytical methods for the near and far field are Campbell[43] and Tanaka et al.[44], respectively.Lateral distances for each propagation distance are chosen large enough to cover all the power.For these large lateral ranges, the MEM is effectively failing, generating zero amplitude and phases from lateral ranges that are about 3 times the spot size of the highest mode in the decomposition.

Figure 2 .
Figure2.MEM performance for the very same test case shown in Fig.1but for smaller lateral ranges.Shown is the performance of a MEM with different maximum mode orders N = 10, 20, 30, 40 and 50, for an incoming circular Gaussian beam with a 2 mm waist being clipped by a 0.5 mm radius circular aperture.Shown are the amplitude (absolute value), phase and relative error ε rel (N R , R, z) at different propagation distances z = 5 mm, 20 mm, 100 mm (near field) and 1000 mm (far field) after the clipping aperture.The analytical methods for the near and far field are Campbell[43] and Tanaka et al.[44], respectively.One can see that the further the beam propagates, or the higher the mode order is, the better is the performance of the MEM.

Fig. 4 .Figure 4 .
Figure 4. Illustration of the waist scaling factor f ws .If f ws = 0.5, the grid beams do not touch; If f ws = 1, the grid beams precisely touch; if f ws = 4 3 , the grid beams intersect; if f ws = 2, the grid beam waist radius and grid distance are equal.
in the MEM example shown in Fig.1, the GBD results with the large lateral ranges shown in Fig.5are not

Figure 5 .
Figure 5. Performance of a GBD for large lateral ranges with different grid sizes G = 100 × 100, 200 × 200, 500 × 500, and 1000 × 1000 for an incoming circular Gaussian beam with a 2 mm waist being clipped by a 0.5 mm radius circular aperture.Shown are the amplitude (absolute value), phase and relative error ε rel (N R , R, z) at different propagation distances z = 5 mm, 20 mm, 100 mm (near field) and 1000 mm (far field) after the clipping aperture.The analytical methods for the near and far field are Campbell[43] and Tanaka et al.[44], respectively.Lateral distances for each propagation distance are chosen large enough to cover all the power.For any propagation distance, the performance of the GBD becomes better with the increasing grid size.

Figure 6 .
Figure 6.Same as Fig.5, but for smaller lateral ranges.Shown is the performance of a GBD with different grid sizes for an incoming circular Gaussian beam with a 2 mm waist being clipped by a 0.5 mm radius circular aperture.Shown are the amplitude (absolute value), phase and relative error ε rel (N R , R, z) at different propagation distances z = 5 mm, 20 mm, 100 mm (near field) and 1000 mm (far field) after the clipping aperture.The analytical methods for the near and far field are Campbell[43] and Tanaka et al.[44], respectively.One sees that the further the beam propagates, or the higher the grid size is, the better is the performance of the GBD.

Figure 7
Figure 7. Performace of a GBD with hexagonal and square grid shapes.Shown are the amplitude (absolute value), phase, and relative error distributions at different propagation distances from a circular aperture with 0.5 mm radius.The incoming circular-symmetric Gaussian beam was centered to the aperture and had its 2 mm waist located in the aperture plane.The analytical methods for the near and far field are Campbell[43] and Tanaka et al.[44], respectively.The number of sampling points X is 3001.

Figure 8 .
Figure 8.The amplitude (absolute value), phase, and relative error distributions at propagation distances, z r /1000 , z r (near-field), 1000 z r (far-field) and 3 million kilometers (extreme far-field).The simulation parameters are listed in Table10.The analytical results refer to the complex electric field of the non-clipped Gaussian beam in this case.It can be seen that MEM behaves better than GBD.

Figure 9 .
Figure 9.The wrapped and unwrapped phase of propagation distance 2.9526 km for the non-clipped circular Gaussian beam.

Figure 10 .
Figure 10.MEM representation of a general astigmatic Gaussian beam at different propagation distances.From left to right: amplitude (absolute value), phase and relative error distribution of the MEM.From top to bottom: propagation distance 29.526 mm, 2.9526 m and 295.26 m respectively.All the simulation parameters are listed in Table12.It can be seen from the third column of this figure that MEM shows a good agreement with the initial general astigmatic Gaussian beam.

Figure 11 .
Figure 11.GBD representation of a general astigmatic Gaussian beam at different distances.From left to right: amplitude (absolute value), phase, and relative error distribution of the GBD.From top to bottom: propagation distance 29.526 mm, 2.9526 m and 295.26 m respectively.All the simulation parameters are listed in Table12.It can be seen from the third column of this figure that GBD shows a good agreement with the initial general astigmatic Gaussian beam.

Figure 12 .
Figure 12.Cross-section of the phase obtained with the MEM along y = 0 and z = 295.26m.The left-hand side image shows the phase using 101 sampling points, the centered image for 3001 sampling points, and the right-hand side image again for 3001 sampling points after a phase-tracking algorithm unwrapped the data.One sees that the small side-maxima, which are seen as circles in Fig.10and Fig.11, disappear when a higher sampling rate is used, and indeed a smooth but strongly curved phase profile is restored in the far field.

Figure 13 .
Figure 13.Amplitude (absolute value), phase and relative error distribution of a clipped circular Gaussian beam with propagation distances 5 mm, and 100 mm in the near field, 1000 mm in the far field, and 3 Gm in the extreme far field (from top to bottom).The lateral range in the fourth row is 2 times the local spot size.The simulation parameters are listed in Table14.The analytical methods for the near and far field are Campbell[43] and Tanaka et al.[44], respectively.The graphs in the third column show the relative error, which indicates that the GBD performs better than the MEM. .

Figure 14 .
Figure 14.Amplitude (absolute value), phase and relative error distribution of a clipped circular Gaussian beam at propagation distances 3 Gm.Both the beam waist diameter and aperture diameter are 30 cm.The right-hand side image shows the relative error, which indicates that the GBD performs better than the MEM.

For a good resolution
of our results, we used 201 × 201 = 40401 sampling points and chose for the MEM a mode order of N = 50 and for the GBD a grid size G = 150 × 150.The results computed by the GBD are shown in Fig. 15, and show the expected amplitude profiles (compare e.g.(cf.[50,Fig. 4], [51, Fig. 3], and [52, Fig. 9.8]).

Figure 15 .
Figure 15.Amplitude profile of the first 4 orders of wavefront aberrations described by Zernike polynomials, which is generated by the GBD.Aberrations of order n=1 are called Tilt, n=2 are Astigmatism and Defocus, n=3 are Coma and Trefoil, and n=4 are Tetrafoil, 2nd Astigmatism, and Spherical Aberration.The shown images cover 200 m × 200 m each.

Figure 16 .
Figure 16.Amplitude difference profile between the MEM and GBD of the first 4 orders of wavefront aberrations described by Zernike polynomials.The shown images cover 200 m × 200 m each.

Figure 17 .
Figure 17.Amplitude, phase distribution and relative error after the beam reflected from mirrors with different curvatures.The graphs in each row represent the amplitude, phase and relative error respectively from left to right, and each row represents the different curvature of mirrors 0 mm −1 ,−0.002 mm −1 , and −0.02 mm −1 from top to bottom.The simulation parameters are listed in Table17.The relative deviation of the GBD from the Gaussian beam and MEM increases when the curvature of the mirror increases.

Figure 18 .
Figure 18.Amplitude profile of the initial Gaussian Beam, MEM and GBD through longitudinal sections near the focal point after reflected from a concaved spherical mirror with the curvature of −0.1 mm −1 .The Amplitude is scaled by log(Amplitude+1).Only the GBD shows signs of spherical aberration.

Table 1 .
Parameters list of the MEM example.

Table 3 .
Parameters list of the GBD example.

Table 6 .
Computational effort of the MEM with increasing mode order.Shown are the decomposition, propagation, superposition, and the resulting total time.The number of sampling points in the target plane is 3001.All times are elapsed real-time, except for the last column, which shows total CPU time.

Table 7 .
Computational effort of the GBD with increasing grid size.Shown are the decomposition, propagation, superposition, and the resulting total time.The number of sampling points in the target plane is 3001.All times are elapsed real-time, except for the last column, which shows total CPU time.

Table 8 .
Computational effort of the MEM with increasing mode order.Shown are the decomposition, propagation, superposition, and the resulting total time.The number of sampling points in the target plane is 101 × 101.All stated times are again elapsed real-time, except for the last column.

Table 9 .
Computational effort of the GBD with increasing grid size.Shown are the decomposition, propagation, superposition, and the resulting total time.The number of sampling points in the target plane is 101 × 101.All stated times are again elapsed real-time, except for the last column.

Table 10 .
Parameters list for non-clipped circular Gaussian beam.

Table 11 .
The discretized NMSE and the summed relative error of the MEM and GBD respectively for different propagation distances.The number of sampling points X is 3001.

Table 12 .
Parameters list for non-clipped general astigmatic Gaussian Beam.

Table 13 .
The DNMSE and the summed relative error of the MEM and GBD for non-clipped general astigmatic Gaussian beam at different propagation distances.The number of sampling points X is 10201.

Table 14 .
Parameters list for circular Gaussian beam clipped by a circular aperture centered in the beam waist.

Table 15 .
The DNMSE and the summed relative error of the MEM and GBD respectively for different propagation distances.The number of sampling points is 3001.

Table 16 .
Parameters list of the wavefront aberration.