## Abstract

Precision optics have been widely required in many advanced technological applications. X-ray mirrors, as an example, serve as the key optical components at synchrotron radiation and free electron laser facilities. They are rectangular silicon or glass substrates where a rectangular Clear Aperture (CA) needs to be polished to sub-nanometer Root Mean Squared (RMS) to keep the imaging capability of the incoming X-ray wavefront at the diffraction limit. The convolutional polishing model requires a CA to be extended with extra data, from which the dwell time is calculated via deconvolution. However, since deconvolution is very sensitive to boundary errors and noise, the existing surface extension methods can hardly fulfill the sub-nanometer requirement. On one hand, the figure errors in a CA were improperly modeled during the extension, leading to continuity issues along the boundary. On the other hand, uncorrectable high-frequency errors and noise were also extended. In this study, we propose a novel Robust Iterative Surface Extension (RISE) method that resolves these problems with a data fitting strategy. RISE models the figure errors in a CA with orthogonal polynomials and ensures that only correctable errors are fit and extended. Combined with boundary conditions, an iterative refinement of dwell time is then proposed to compensate the errors brought by the extension and deconvolution, which drastically reduces the estimated figure error residuals in a CA while the increase of total dwell time is negligible. To our best knowledge, RISE is the first data fitting-based surface extension method and is the first to optimize dwell time based on iterative extension. An experimental verification of RISE is given by fabricating two elliptic cylinders (10 mm × 80 mm CAs) starting from a sphere with a radius of curvature around 173 m using ion beam figuring. The figure errors in the two CAs greatly improved from 204.96 nm RMS and 190.28 nm RMS to 0.62 nm RMS and 0.71 nm RMS, respectively, which proves that RISE is an effective method for sub-nanometer level X-ray mirror fabrication.

© 2021 Optical Society of America under the terms of the OSA Open Access Publishing Agreement

## 1. Introduction

With the rapid development of precision technologies, the demand of high-precision optics has greatly increased in various cutting-edge applications, such as telescopes for space exploration [1–3], optics in EUV lithography [4,5], and X-ray mirrors for synchrotron radiation and free-electron laser facilities [6–10], such as National Synchrotron Light Sources II (NSLS-II). As shown in Fig. 1(a), X-ray mirrors are used to deflect and focus the incoming X-ray beams. They are usually made of Si or fused silica in a rectangular shaped substrate. As shown in Fig. 1(b), a smaller effective region called the Clear Aperture (CA) is usually polished at the center of the mirror. Due to the grazing incidence geometry, the tangential dimension (*i.e.* length $L = 100$ mm ~ $1$ m) of a CA is much longer than its sagittal dimension (*i.e.* width $W = 5$ mm ~ $20$ mm). Depending on the mirror application, the mirror surface shape can be flat, circular cylinder, elliptic cylinder, toroid, ellipsoid, *etc.* With the rapid evolution of third and fourth generation synchrotron X-ray light sources, the requirement of precision and smoothness of X-ray mirrors has greatly increased. Single or even sub-nanometer figure errors are required to avoid destruction of the incoming wavefront and keep its imaging capability at the diffraction limit [11,12]. Figure 1(c) illustrates the grazing angle $\theta$ and the figure errors $z(x)$ along the $x$ direction, from which the wavefront error $\varphi (x)$ is calculated as $\varphi (x)=-2\mathrm {sin}\theta \cdot z(x)$. According to Maréchal’s criterion [13], the wavefront diffraction limit is achieved when $\varphi (x)$ is less than or equal to $\lambda /14$ Root Mean Square (RMS). This requires, as a first approximation, $z(x)$ to be within $\lambda /(28\times \theta )$ for a small $\theta$. Assuming that $\lambda =0.1$ nm and $\theta =3$ mrad (typical parameters for hard X-ray), $z(x)$ should be approximately 1 nm RMS. Such a high level of surface quality cannot be achieved by conventional mechanical polishing techniques.

In fact, precision optical surfaces are always manufactured with Computer Controlled Optical Surfacing (CCOS) processes, including small-tool polishing [14], bonnet polishing [15], magnetorheological finishing [16], Ion Beam Figuring (IBF) [6,9,10,17–21], *etc*. Different CCOS processes can be adopted based on the requirements of the precision and the shapes of the desired optical surfaces. IBF, as one ultra-precision CCOS technique, has been frequently used for optical surface finishing [9,10,17,18,20–23]. As shown in Fig. 1(c), IBF removes material particles from an optical surface via physical sputtering. Compared with conventional small-tool CCOS techniques, IBF has the advantages of being non-contact, with no mechanical load force, minimal surface or subsurface damage, and no tool overhang. The Optical Metrology and Fabrication Group at NSLS-II has been working on developing IBF methods and systems [9,10,18,20,21,23] since early 2015. The investigation started with One-Dimensional IBF (1D-IBF) [18,20,21,23], then a Two-Dimensional IBF (2D-IBF) system [9,10] was successfully developed.

IBF follows the convolutional polishing model,

where "$*$" represents the convolution operation, $b(x,y)$ is the Beam Removal Function (BRF), $\widehat {z}(x,y)$ is the removed material, and $t(x,y)$ is the dwell time. $\widehat {z}(x,y)$ and $b(x,y)$ are known, $t(x,y)$ is thus calculated via deconvolution [10,24,25], which is an ill-posed operation that is very sensitive to boundary errors and noise.Equation (1) implies that, to achieve the desired figure errors in a CA, $t(x,y)$ is required to be calculated on a Dwell Grid (DG) that is larger than the outline perimeter of the CA with the radius $r$ of $b(x,y)$. Therefore, we can rewrite Eq. (1) as

To fulfill the requirement given in Eq. (2), in our previous 2D-IBF solution, the entire mirror was measured using our Stitching Interferometry (SI) platform [26,27] and the dwell time optimization was performed on the DG using the Robust Iterative Fourier Transform-based dwell time Algorithm (RIFTA) [10]. As shown in Fig. 2, this strategy has been successfully applied to reduce the figure errors in a rectangular CA from 3.69 nm RMS to 0.49 nm RMS, where the figure errors to be corrected are small. However, when the figure errors are large, this strategy is no long appropriate. The reason is manifested in Fig. 3.

Figure 3 gives the SI measurements during the fabrication of an elliptic cylinder from a spherical mirror (details of the experiment are explained in Section 4.1), where the figure errors are much larger than those shown in Fig. 2. The fringes in one Field Of View (FOV) of the initial sphere are shown in Fig. 3(a) and the corresponding fringes in the same FOV after one IBF run are given in Fig. 3(b). After the IBF run, the fringes become much denser due to the large figure error differences (approximately 1000 nm) inside and outside the DG, resulting in the large boundary errors in the DG map shown in Fig. 3(c) stitched from multiple FOVs in Fig. 3(b). Since deconvolution is very sensitive to boundary errors, if $t(x,y)$ for the next IBF run is deconvolved from this DG map, the solution will be unreliable and thus can hardly be used to achieve the sub-nanometer figure error convergence in the CA.

In fact, a more proper way to calculate $t_{DG}(x,y)$ is to only use the figure errors in a CA and to rely on surface extension to obtain the extra data outside. Therefore, a robust surface extension method becomes the key to the success of the entire fabrication process. The existing surface extension methods include the Papoulis-Gerchberg extension [28–33], the nearest-neighbors extension [19,31], the Gaussian extension [19,31,34], and the smooth extension [19]. The Papoulis-Gerchberg extension is a frequency domain method. It assumes figure errors to be band-limited and Fourier transform is used to iteratively expand the error spectrum. The nearest-neighbors extension, Gaussian extension and smooth extension directly extend the surface in spatial domain. The nearest-neighbors extension generates an extension data point by interrogating its nearest neighbors. If and only if there are at least three valid neighboring data points around an extension point, it is assigned with the value as the average of these neighbors. The Gaussian extension extends a CA by applying a Gaussian function between a certain CA boundary point and the corresponding DG boundary point. The smooth extension guarantees a piece-wise $C^{1}$ continuity at the extension boundary using triangle-based interpolation [35]. Zhou *et. al.* also proposed a descending profile that can be multiplied to the extended regions to drag down the extended errors and thus decrease the total dwell time [19]. Although these surface extension methods have been widely adopted in CCOS of circular mirrors, their direct application to rectangular CAs in X-ray mirrors, however, can hardly achieve the sub-nanometer RMS figure error requirement, exposing the following two challenges.

First, the known figure errors in a CA should be properly modeled to guarantee the continuity along extension boundary, especially when the CA of an X-ray mirror has different shapes in its tangential and sagittal dimensions. Second, high-frequency errors and noise that cannot be corrected by a particular BRF should not be extended. Improper extension will bring unsatisfactory boundary data and unexpected errors to the extended DG, which further result in unreliable solutions of $t_{DG}(x,y)$ from the ill-posed deconvolution and fail the fabrication process.

Based on Eq. (1), material removal is analogous to a low-pass filtering process, where corrections are more effective on lower-frequency figure errors. Inspired by this concept, we propose a Robust Iterative Surface Extension (RISE) method that solves the above two problems and optimizes dwell time with a polynomial-based extension strategy. First, RISE models the figure errors in a CA using orthogonal polynomials. The polynomial orders vary according to the different tangential and sagittal shapes of a CA, eliminating the continuity issues at the boundary. Second, the polynomials ensure that only correctable errors are fit and extended. In addition, boundary conditions are added during the polynomial fitting to balance the trade-off between the desired figure errors in a CA and the total dwell time. Third, an iterative refinement of dwell time is performed to compensate for the unreliable solutions brought by the extension and the ill-posed deconvolution. In each iteration, a smart piston adjustment is employed to eliminate the possible extra removal and thus reduce the total dwell time.

To verify the effectiveness of the RISE concept, an experiment is performed where the main goal is to fabricate two 10 mm $\times$ 80 mm elliptic cylindrical surfaces starting from a spherical mirror using IBF. The initial Radius of Curvature (ROC) of the spherical mirror is around 173 m, and the target ellipse parameters are the object distance $p=30$ m, image distance $q=0.3$ m, and grazing angle $\theta =3$ mrad. The figure errors in the two CAs are reduced from 204.96 nm RMS and 190.28 nm RMS to 0.62 nm RMS and 0.71 nm RMS, respectively. Also, the surface roughness before and after the IBF runs for both the CAs remain invariant, which prove that RISE is an effective method for sub-nanometer X-ray mirror fabrication, especially when the figure errors to be corrected in a CA is large.

The rest of the paper is organized as below. Section 2 reviews the principles of the existing surface extension methods, followed by an explanation of the proposed RISE concept in Section 3. The experimental verification of RISE is then given in Section 4. Section 5 discusses the limitations and the generalization of the proposed method. Section 6 concludes the paper.

## 2. Existing surface extension methods

For a rectangular CA, as shown in Fig. 4, surface extension can be defined as extending the outline perimeter of the initial figure errors $z_{CA}(x, y)$ in the CA with the radius $r$ (*i.e* four grids in Fig. 4) of a BRF. Therefore, for a CA size of $L \times W$, the size of the extended figure error map $z_{DG}(x, y)$ in the DG is $\left (L+2r\right )\times \left (W+2r\right )$. In this section, the application of the existing surface extension methods, including the Papoulis-Gerchberg extension [28,29], the nearest-neighbors extension [19,31], the Gaussian extension [19,31], the smooth extension [19,36], and descending multiplier [19,36], to a rectangular CA is reviewed.

It is worth mentioning that, two criteria are always used to decide whether a dwell time solution $t_{DG}(x,y)$ calculated with a certain surface extension method is reasonable [19,31,36]. One criterion is the estimated figure error residuals in a CA, $e_{CA|est}(x,y)$, expressed as

The other criterion is the total dwell time that is calculated as the summation of all the entries in $t_{DG}(x,y)$. Ideally, it is expected that $e_{CA|est}(x,y)$ is minimized with the shortest total dwell time.

#### 2.1 Frequency domain Papoulis-Gerchberg extension method

The Papoulis-Gerchberg extension [28,29] iteratively extends a CA in frequency domain. It assumes that the figure errors $z_{DG}\left (x,y\right )$ are band-limited. Let $z_{DG}(x,y)$ have a spectrum

where $\mathscr {F}$ represents the Fourier transform. The spectrum is assumed band-limited and is nonzero only over some region $\Omega$ so that a spectral pupil can be defined asSince $z_{DG}$ is only known within $z_{CA}$, the corresponding aperture $G_{T}$ in spatial domain is defined by

The Papoulis-Gerchberg extension then iteratively extrapolates the unknowns in $z_{DG}$ as follows: (1) Fourier transform $z_{DG}G_{T}$; (2) multiply by the spectral pupil $G_{\Omega }$; (3) inverse Fourier transform; (4) discard the known figure error in $z_{CA}$ by multiplying by $\left (1-G_{T}\right )$; (5) add in the known figure error $z_{DG}G_{T}$; and (6) go to Step 1, and repeat. From the above procedure, the $i\mathrm {th}$ estimate of $z_{DG}$ can be written as

where $z_{DG}^{0}=z_{DG}G_{T}$ and the band-limiting operator $B_{\Omega }$ is defined as $B_{\Omega }=\mathscr {F}^{-1}G_{\Omega }\mathscr {F}$, where $\mathscr {F}^{-1}$ represents the inverse Fourier transform. Although the convergence of Eq. (7) has been proved [29], the convergence rate is slow and many iterations are required to obtain a good estimation of $z_{DG}$. Otherwise, the edges of $z_{CA}$ cannot be smoothly extended. More importantly, when the tangential and sagittal dimensions of a CA have different shapes, it is hard to determine a proper spectral pupil $G_{\Omega }$, in which case either middle to high-frequency errors will be introduced or even more iterations will be required.#### 2.2 Spatial domain methods

**Nearest-neighbors extension**. The nearest-neighbors extension simply calculates the error for an extension point in the DG using the known figure errors at its $N$ nearest neighbors. As shown in Fig. 4, the figure error for a point $D(x_{i}, y_{i})$ in the DG is calculated as

Since not all the extension points in the DG have valid neighbors that already have values assigned, the nearest-neighbors extension is thus required to be performed iteratively [31]. In each iteration, if and only if there are at least three valid neighboring data points around an extension point, it is assigned with the average of the figure errors at these neighbors using Eq. (8). The iteration is performed until all the points in the DG are assigned with a figure error. The nearest-neighbors extension is simple and easy to implement. However, it always leaves obvious artifacts along the boundary of $e_{CA|est}(x,y)$, since the local averaging operations cannot model $z_{CA}$ and its boundary continuity properly. Also, the result is sensitive to noise when $N$ is small.

**Gaussian extension**. The Gaussian extension is commonly used in fabrication of circular-shaped mirrors [19,31,34], since the removal functions are always axially symmetric and Gaussian shaped. As shown in Fig. 4, to obtain the extended value $z_{DG}\left (x_{i}, y_{i}\right )$ at the point $D(x_{i}, y_{i})$ in the DG, its closest CA boundary point $C(x_{p}, y_{p})$ is first identified. The figure error at $D(x_{i}, y_{i})$ is then calculated as

*i.e*point $C(x_{p}, y_{p})$ in Eq. (9)). The continuity at the entire boundary of a CA thus cannot be guaranteed. Especially, when the CA boundary point is noise, the extension is even worse and the dwell time cannot be reliably calculated, $e_{CA|est}(x,y)$ will thus be large.

**Smooth extension method**. The smooth extension [36] uses a Delaunay triangulation of the points in a CA to perform interpolation [35]. Piece-wise $C^{1}$ continuity is satisfied at each cubic surface of a triangle. Therefore, the smooth extension results in a continuous surface near the boundary of a CA so that the dwell time is always more accurate than those calculated using the other methods. However, this method triangulates the entire CA, and the uncorrectable high-frequency errors and noise are also interpolated, which prevents $e_{CA|est}$ from reaching $\leq 1$ nm RMS. Also, this method greatly increases the total dwell time if the figure errors along the boundary of a CA are large [36]. A Gaussian function given in Eq. (9) can potentially be multiplied to the extension region to lower down the extended figure errors. However, the Gaussian extension falls so fast that a large $e_{CA|est}(x,y)$ is produced [36]. To obtain smaller $e_{CA|est}(x,y)$ with relatively short dwell time, a slower descending multiplier [19] is proposed as

## 3. Robust iterative surface extension method

Based on the description above, the factors that affect the estimated figure error residuals in a CA, $e_{CA|est}(x,y)$, and the total dwell time are summarized as the followings.

The deepest impact to $e_{CA|est}(x,y)$ comes from the inappropriate modeling of the figure errors, especially the continuity along the boundary of a CA is not well handled, leading to incorrect dwell time solutions and large $e_{CA|est}(x,y)$. Even worse, if a CA initially has large figure errors along the boundary, the extension will amplify these errors, resulting in longer total dwell time than necessary. In addition, high-frequency errors and noise influence the magnification of $e_{CA|est}(x,y)$ if they are extended. These errors mainly come from the noise of metrology instruments or the footprints of previous polishing tools, and they are not correctable by a particular BRF. Added with the fact that deconvolution itself in dwell time calculation is ill-posed and thus is very sensitive to boundary errors and noise, it is always hard to obtain a reliable dwell time solution. To reduce the effects of these factors, the RISE method is proposed. As schematically illustrated in Fig. 5, RISE includes two phases. In the initial extension phase, $z_{CA}$ is modeled and extended with orthogonal polynomials, and the initial $t^{0}_{DG}(x,y)$ is calculated by RIFTA [10]. However, $t^{0}_{DG}(x,y)$ obtained from the initial extension phase may still be sub-optimal due to the errors introduced during the polynomial fitting and extension processes. Therefore, an iterative refinement phase is added to progressively refine the dwell time solution until a desired $e_{CA|est}(x,y)$ is achieved. In the rest of this section, the principle of RISE is explained in detail.

#### 3.1 Surface extension as polynomial fitting with boundary conditions

Surface extension is in fact analogous to extrapolation, in which it is required to estimate the unknown $z_{DG}$ based on known distribution in $z_{CA}$. The extended $z_{DG}$ is then followed by a deconvolution process to calculate $t_{DG}$. To achieve small $e_{CA|est}$ in a CA, $z_{CA}$ should be well modeled so that the continuity at the CA boundary is guaranteed and the extended $z_{DG}$ is smooth.

In fact, this requirement can be fulfilled if $z_{CA}$ can be described in an analytical form. In wavefront aberration analysis, orthogonal polynomials, including the Zernike polynomials [37] for circular aperture, and the Legendre polynomials [38] and Chebyshev polynomials [39] for rectangular apertures, have been used as an excellent fitting basis to analytically decompose and interpret the test data. In this study, both Legendre polynomials and Chebyshev polynomials can be used to model the known $z_{CA}$ and the extension process. A set of 2D polynomials that are orthogonal over a unit square can be described as

where $j$ is a polynomial ordering index starting with $j=1$, and $m$ and $n$ are the orders of the corresponding 1D polynomials $P_{m}(x)$ and $P_{n}(y)$ in $x$ and $y$ directions, respectively [38,39]. Their orthonormality is expressed by where $K$ is a normalization factor, and $\delta _{jj'}$ is the Kronecker delta.To fit $z_{DG}(x,y)$ using Eq. (11), its $x$ and $y$ coordinates are normalized to $[-1, 1]$ by the dimensions of the DG $(L+2r, W+2r)$ (see Fig. 4). We define a boundary condition on a DG as the outermost perimeter of $z_{DG}(x,y)$, *i.e* $x,y=$ -1 or 1, and denote it as $\partial z_{DG}(x,y)$. In the fitting process, besides the known $z_{CA}(x,y)$, $\partial z_{DG}(x,y)$ is also included to tune the fitting coefficients and the extended data.

The simplest way is to set $\partial z_{DG}(x,y)=0$, however, it is so strict that may cause sudden jumps in the extension region, which will in turn affect the dwell time calculation and increase $e_{CA|est}$. Therefore, we initialize $\partial z_{DG}(x,y)$ as the boundary of $z_{DG}(x,y)$ after the smooth extension, *i.e* $\partial z_{DG}(x,y) = \partial z_{DG|SE}(x,y)$, where $z_{DG|SE}(x,y)$ represents the figure errors in the DG extended using the smooth extension method, since it guarantees the $C^{1}$ surface continuity. As mentioned in Section 2, however, the smooth extension may provide unsatisfactory extension and increase the total dwell time when the figure errors at the boundary of $z_{CA}(x,y)$ are large. We thus introduce a factor $k, k\in \left [0,1\right ]$ that is multiplied to $\partial z_{DG|SE}(x,y)$ to flexibly adjust the figure errors in the extension region and consequently tune the dwell time. Also, an iterative scheme that will be explained in Section 3.3 is added to further refine the dwell time solution.

Based on the above description, the optimal polynomial coefficients $\widehat {c}_{j}$ used in the polynomial fitting are solved from the following least-squares problem

The initial extension is thus obtained as $z^{0}_{DG}=\sum _{j}\widehat {c}_{j}Q_{j}$, after which the initial dwell time $t^{0}_{DG}$ is calculated using RIFTA [10]. It can be expected from Eq. (14) that the shortest total dwell time is obtained at $k=0$ while the minimum $e_{CA|est}(x,y)$ is achieved at $k=1$. In practice, $k$ can be adjusted based on the desired $e_{CA|est}(x,y)$. One simple way is to initialize $k$ as $k=0$ and progressively increase it until the desired $e_{CA|est}(x,y)$ is achieved.

#### 3.2 Determination of the polynomial orders

The chosen polynomials in Eq. (11) are separable in the Cartesian coordinates $x$ and $y$ of a point in the unit square. Therefore, the polynomial orders in the $x$ and $y$ directions can be determined individually. This is especially useful when the figure errors in a CA have different shapes in its tangential and sagittal dimensions.

Ideally, the polynomial order that can perfectly describe a certain figure error map is equal to the number of pixels in the map. However, due to the low-pass convolutional polishing model in Eq. (1), the over-fit high-frequency errors will affect the reliability of the calculated dwell time. Therefore, it is always necessary to only fit and extend the lower-frequency errors. The residuals between $z_{CA}$ and $\left [z^{0}_{DG}\right ]_{CA}$ can be used to guide the selection of the appropriate polynomial orders. The polynomial orders in $x$ and $y$ directions are increased from zero until the RMS of the residual is less than or equal to a threshold $\gamma$. A good candidate value for $\gamma$ is the random error of a particular metrology instrument. In this study, $\gamma = 0.3$ nm is employed, which is equal to the RMS of the random errors of the SI platform.

#### 3.3 Iterative refinement of dwell time

Even with the proposed polynomial-based extension, the initial $t^{0}_{DG}$ may still be sub-optimal due to the errors introduced during the extension and the unreliable, ill-posed deconvolution. Therefore, an iterative refinement scheme is employed to allow the dwell time $t^{i}_{DG}$ to progressively approach a better solution.

The basic idea of the iterative refinement is that, in the $i$th ($i\geq 1$) iteration, the polynomial-based extension is applied to the residual $e_{CA|est}$ to obtain $z^i_{DG}$. The resulted $\Delta t_{DG}$ is then added to the final dwell time solution $t^{i}_{DG}$ as

where $h$ is added to control the convergence rate of $t^{i}_{DG}$. The fastest rate is achieved at $h=1$, however, overshoot may occur and convergence would not be reached. On the contrary, if $h$ is too small, the convergence can always be reached but the convergence rate will be slow. In this study, $h = 0.5 \sim 0.8$ is found to be a proper range that allows RISE to reach the convergence in less than ten iterations.It is worth noting that, the extended $z^{i}_{DG}$ is not directly used to calculate $\Delta t_{DG}$. Instead, a piston adjustment is performed on $z^{i}_{DG}$ and only the resulted $\Delta z_{DG}$ is used in RIFTA to calculate $\Delta t_{DG}$ and update $t^{i}_{DG}$ as

where in which $\mathrm {min}(\cdot )$ represents the minimum value in "$\cdot$", and $\widehat {z}^{i-1}_{DG} = t^{i-1}*b$ is the estimated material removed in the DG for the $(i-1)$th iteration. This piston adjustment eliminates the extra removal added in each iteration. It relaxes the non-negativity requirement of $\Delta t_{DG}$ but only guarantees $t^{i}_{DG}$ to be non-negative. To our best knowledge, such an algorithmic level iterative solution has not been attempted in the literature of precision optical fabrication. As will be shown in Section 4.1, it help us achieve the desired sub-nanometer level estimation before IBF, which has the potential to reduce the number of iterations required in the real fabrication process and thus increase the overall efficiency.#### 3.4 Flow of the RISE method

The flow chart of the RISE method is schematically illustrated in Fig. 5. The initial extension phase includes Steps $1 \sim 4$. In Steps $1 \sim 3$, the polynomial-based extension described in Section 3.1 is performed to extend $z_{CA}$ and obtain the initial $z^{0}_{DG}$, from which the initial dwell time $t^{0}_{DG}$ is calculated by RIFTA in Step 4.

Afterwards, an iterative refinement of $t^{0}_{DG}$ is performed in Steps $5 \sim 15$ as follows.

- • Step 5: Initialize $\Delta t_{DG}=0$ and the iteration number $i=1$.
- • Steps 6 and 7: In Step 6, the estimated material removal in the DG, $\widehat {z}^{i-1}_{DG}$, which will be used in Step 12 for piston adjustment, is calculated from the convolution between $t^{i}_{DG}$ and $b$. Based on $\widehat {z}^{i-1}_{DG}$, the residual $e_{CA|est}$ is updated using Eq. (3) in Step 7.
- • Step 8: If $e_{CA|est}$ satisfies the convergence criterion $\mathrm {RMS}\left (e_{CA|est}\right )\leq \epsilon$, the iteration ends. Otherwise, it goes to Step 9. The convergence threshold is set as $\epsilon = \gamma = 0.3$ nm.
- • Steps $9\sim 12$: The iteration number $i$ increases by 1 and the polynomial-based extension is used to extend $e_{CA|est}$. Different from the initial extension, in Step 9, an additional polynomial fitting is performed on $e_{CA|est}$ to obtain $e_{CA|est\_fit}$. The boundary condition is then generated by the smooth extension on $e_{CA|est\_fit}$ instead of $e_{CA|est}$ in Step 10. The reason for adding this additional fitting is that $e_{CA|est}$ is already dominated by high-frequency errors. If the boundary condition is obtained by directly extending $e_{CA|est}$ using the smooth extension, it will be affected by the high-frequency components. Therefore, to only extend the remaining lower-frequency errors, $e_{CA|est}$ is fit with polynomials first and the boundary condition is obtained from the extension of $e_{CA|est\_fit}$. Afterwards, in Step 11, $e_{CA|est\_fit}$ is replaced by $e_{CA|est}$ and the polynomial-based extension is performed to obtain $z^{i}_{DG}$ in Step 12. It is worth mentioning that, the same polynomial orders obtained from the initial extension can be applied to both of the polynomial fittings if the same metrology instrument is used, since the error level remains invariant.
- • Steps $13 \sim 15$: The piston adjustment and dwell time update processes described in Section 3.3 are performed.

## 4. Experiment

To verify the effectiveness of the RISE concept, as shown in Fig. 6, an experiment was performed where the main goal is the fabrication of two 10 mm $\times$ 80 mm elliptic cylinders (*i.e* CA1 and CA2) starting from a pitch-polished silicon spherical mirror with ROC around 173 m using our 2D-IBF system [9,10]. The ion source parameters are beam voltage = 600 V, beam current = 10 mA, accelerator voltage = -90 V, and accelerator current = 2 mA. The target ellipse parameters for both the CAs, as shown in Fig. 7(a), are $p = 30$ m, $q = 0.3$ m, and $\theta =3$ mrad. The metrology data are provided by the SI platform [26,27]. The lateral resolution of SI is set as 0.09 mm/pixel and the IBF machining interval is 0.3 mm. The figure errors from the initial sphere to the target elliptic cylinders in CA1 and CA2 are 204.96 nm RMS and 190.28 nm RMS, respectively.

For CA1, the dwell time is initially calculated using the metrology data of the corresponding DG without using RISE. However, the figure errors in CA1 cannot be improved to sub-nanometer level after two IBF runs. Therefore, RISE is then applied to correct the figure error residuals in CA1. For CA2, only the CA data is used and the RISE method is applied in all the IBF runs. Chebyshev polynomials are used in RISE for the experiment, since they are computationally more efficient than Legendre polynomials. The same BRF shown in Fig. 7(b) is used in each IBF run. The peak removal rate of the BRF is 7.5 nm/s and its full width at half maximum is 4.2 mm.

The rest of this section is organized as follows. Section 4.1 studies and compares the performances of the existing surface extension methods by applying them to calculate the dwell time and the estimated residual for CA2. The experimental results in CA1 and CA2 are explained in Sections 4.2 to 4.4. Lastly, to quantify the performance of the fabricated mirrors, we have simulated their diffraction limited capabilities in Section 4.6 using the Synchrotron Radiation Workshop (SRW) [40,41].

#### 4.1 Existing surface extension methods vs. RISE

To compare the performances of the existing surface extension with the proposed RISE method, they are applied to calculate the dwell time and estimated figure error residuals for CA2. Figure 8 demonstrates the extended figure error maps, $z_{DG2}$, the dwell time maps, $t_{DG2}$, and the estimated figure error residual maps in CA2, $e_{CA2|est}$, calculated with the Papoulis-Gerchberg extension (Fig. 8(a), the nearest-neighbors extension (Fig. 8(b)), the Gaussian extension (Fig. 8(c)), the smooth extension (Fig. 8(d)), the smooth extension multiplied by descending multiplier (Fig. 8(e)), and RISE (Fig. 8(f)).

It is obvious that the boundary of CA2 is not smoothly extended with the Papoulis-Gerchberg extension, the nearest-neighbors extension, and the Gaussian extension. The Papoulis-Gerchberg extension is performed for 500 iterations. However, as shown in Fig. 8(a), the boundary of CA2 is still not well handled, since the errors along the boundary have higher frequencies than those inside and the band-limited expansion is more effective on lower frequencies. As shown in Fig. 8(b), the nearest-neighbors extension smooths the boundary by the local averaging operations, thus it arrives at a lower $e_{CA2|est}$. The Gaussian extension shown in Fig. 8(c) achieves the shortest dwell time due to the descending profile multiplied to the extension region of CA2. Nonetheless, the Gaussian function descends so fast that the error frequencies along the boundary of CA2 is high and thus cannot be corrected.

The smooth extension shown in Fig. 8(d) guarantees the piece-wise $C^{1}$ continuity so that the boundary of CA2 is more appropriately handled. The RMS of $e_{CA2|est}$ reaches 1.5 nm with a slightly longer total dwell time (*i.e* 9.3 mins) than the nearest-neighbors extension. The descending multiplier $m(l)$ given in Eq. (10) is then multiplied to $z_{DG2}$ extended by the smooth extension. As shown in Fig. 8(e), the errors along the boundary of CA2 is dragged down and the total dwell time becomes smaller, while the RMS of $e_{CA2|est}$ increases.

Therefore, it is found that the smoothness at the boundary of CA2 is the most important factor of achieving a lower $e_{CA2|est}$. On the other hand, a descending profile along the boundary, such as the Gaussian extension or the descending multiplier approaches, helps reduce the dwell time but the RMS of $e_{CA2|est}$ will increase in return, since the high-frequency errors along the boundary may not be correctable.

The proposed RISE method is then applied to CA2. Figure 8(f) demonstrates the results of the initial extension followed by the iterative refinement. The factor $k=1$ (see Steps 2 and 10 in Fig. 5) is used for RISE, since we focus more on a small $e_{CA2|est}$. The iteration criterion is set as $\epsilon =0.3$ nm, which is equal to the random error of the SI platform. The polynomial orders in $x$ and $y$ directions obtained with the method explained in Section 3.2 are 6 and 10, respectively. After the initial extension, the RMS of $e_{CA2|est}$ is 0.88 nm, which is further reduced to 0.26 nm after the iterative refinement. Also, thanks to the piston adjustment in Step 12 (see Fig. 5), only 1.5 mins additional dwell time is added to the initial dwell time solution. It is also worth mentioning that, as shown by the magnified parts of the dwell time maps in Figs. 8(a) to 8(e), the high-frequency errors and noise are left in the dwell time maps, since they are also extended with the existing surface extension methods. These errors are eliminated by RISE as shown in Fig. 8(f).

#### 4.2 Experimental results in CA1

For comparison purpose, in CA1, RISE is not applied at first. Instead, the metrology data in the DG is used directly to calculate the dwell time. Figure 9(a) shows the estimated figure error residuals $e_{CA1|est}$ and the real figure error residuals $e_{CA1|real}$ in CA1 for the first two IBF runs.

The estimated figure error residuals $e_{CA1|est}$ for the 1st IBF run are estimated as 0.64 nm RMS. It achieves the sub-nanometer level since the boundary has not been affected by the IBF process. After the 1st IBF run, the real figure errors in CA1 are reduced from 204.96 nm RMS to 3.66 nm RMS. The convergence rate is 98.2%. The same dwell time calculation is performed again for the 2nd IBF run, however, due to the unreliable measurement and the large errors at the boundary of the DG (see Section 1 and Fig. 3) after the 1st IBF run, the RMS of $e_{CA1|est}$ is estimated as 2.88 nm, and $e_{CA1|real}$ is barely reduced from 3.66 nm RMS to 3.42 nm RMS.

To further reduce the figure errors in CA1, the proposed RISE method is then applied to calculate the dwell time for the 3rd IBF run. As shown in Fig. 9(b), $e_{CA1|est}$ estimated by RISE is 0.59 nm RMS, and the final $e_{CA1|real}$ after the 3rd IBF run achieves 0.62 nm RMS, which closely duplicates the RISE estimation.

#### 4.3 Experimental results in CA2

In CA2, RISE is applied to calculate the dwell time for all the three IBF runs shown in Fig. 10. The final $e_{CA2|real}$ reaches 0.71 nm RMS after three IBF runs, which is similar to the result obtained in CA1. It can be found that $e_{CA2|est}$ is getting bigger after each IBF run due to the footprints of the BRF left on CA2. These footprints cannot be corrected unless a BRF with smaller size is used or a smoothing process is followed.

It is also worth mentioning that, the main contribution to the figure errors in the final $e_{CA2|real}$ after the 3rd IBF run is the region of higher error residuals at its bottom-right corner, which also causes that $e_{CA2|real}$ is slightly higher than $e_{CA1|real}$. These errors are a result of the fault of the horizontal translation stage during the 1st IBF run. If this region is excluded, $e_{CA2|real}$ can be improved to 0.67 nm RMS, which is higher than $e_{CA1|real}$ by only 0.05 nm RMS. This small difference is negligible, since it is below the measurement repeatability of SI.

#### 4.4 Estimation vs. real experimental results

Achieving the desired ultra-low figure errors indeed relies on two important factors: accurate dwell time solutions and precise dwell time implementations. The RISE method is proposed to produce robust dwell time solutions, from which the figure error residuals can be reliably estimated. In practice, the estimated residuals should be as small as possible to minimize the errors produced during the calculation stage. The precise implementation of dwell time, however, is attributed to the determinism of the hardware configuration of a polishing system. The higher determinism a system achieves, the real polishing result will be closer to the estimation. In our 2D-IBF system, the determinism is affected by the heating from the ion beam, the motion limits of the translation stage, and the footprint of the BRF.

As shown in Figs. 9(a) and 10, the real figure errors does not converge to the estimation after the 1st IBF runs, due to the heating time of around 74 mins (see Fig. 8(f)). As the figure errors became smaller and thus the required total dwell time decreases, the determinism improves. As shown in Fig. 9(b), the figure error residuals for the 3rd IBF run in CA1 closely duplicate the estimation. The influence of the determinism caused by the translation stage is demonstrated in Fig. 10, where a small positioning error during the 1st IBF run in CA2 brings the artifacts to the final figure error residual map. Also, the BRF footprints are obviously left in the vertical direction on the polished surfaces, which further degrades the level of surface figure accuracy. The reason for the generation of these vertical footprints is that we applied the velocity-based feeding mode only in the horizontal direction. For the vertical direction, a constant step size was fed with the position-based mode.

#### 4.5 Surface roughness

Figures 11(a) and 11(b) shows the Power Spectral Density (PSD) curves and their corresponding integrated PSD RMS curves along the tangential direction in CA2 measured in different spatial frequencies using different metrology instruments. The SI measurement results demonstrates that the IBF process using the proposed RISE method is effective for sub-nanometer level figure error correction.

In addition, the middle-frequency surface roughness before and after the IBF process was measured by a Zygo New View white-light microscope interferometer with 20$\times$ magnification. The spatial resolution of the measurements is 0.55 $\mu$m and the FOV is 0.35 mm $\times$ 0.26 mm. For the high-frequency surface roughness, since our Atomic Force Microscopy (AFM) instrument was not available until the IBF process was completed, a region near CA2 that had not been processed by the ion beam was selected as a reference. The spatial resolution of the AFM measurements is 31.3 nm and the FOV is 7.98 $\mu$m $\times$ 7.98 $\mu$m. We can see that the surface roughness, at the explored resolution levels, is not affected by the IBF process and still remain at a very acceptable level for X-ray mirrors.

#### 4.6 Simulation using SRW

To demonstrate the focusing capability of these two IBF-processed elliptical cylinders, we use the SRW simulation engine [40,41] to perform a virtual X-ray focusing experiment. The layout of the simulation is shown in Fig. 12. The X-ray source is simulated as a coherent Gaussian beam with a 0.1 $\mu$m RMS $\times$ 0.1 $\mu$m RMS beam waist. The fabricated elliptic cylindrical mirror ($p = 30$ m, $q = 0.3$ m, and $\theta =3$ mrad) is used as a horizontal focusing mirror placed at 30 m. Slit apertures are used to illuminate only the mirror’s CA, and a watch point (detector) is placed at 30.3 m where the nominal focal spot is located.

Three 2D figure error maps of the polished CAs from the fabrication process are individually used as the figure errors of the horizontal focusing mirror in the SRW simulation. The first map is the figure errors after the first two IBF runs in CA1 without using RISE (Fig. 9(a)). The second one is the figure error map after the third run CA1 obtained with RISE (Fig. 9(b)). The third one is the figure error map after the third run in CA2 using RISE (Fig. 10).

The focusing ability of these three mirror surfaces are studied with the X-ray wavelength $\lambda$ varying from 0.02 nm to 0.2 nm. We calculated the Strehl ratio of the focal spots and our results are shown in Fig. 13.

As expected, the Strehl ratios [42] of these three mirror surfaces increase with the X-ray wavelength. For $\lambda =$ 0.02 nm, all the Strehl ratios are smaller than 0.8 in Fig. 13(a). Based on Maréchal’s criterion, the Strehl ratios should be close to 0.8 with the two $0.6 \sim 0.7$ nm RMS figure errors around $\lambda =$ 0.05 nm, which shows good agreement with our simulation in Fig. 13(b). Figure 13(c) shows that the Strehl ratio of CA1 without using RISE method is higher than 0.8 for $\lambda =$ 0.2 nm. These simulations clearly indicate that using the RISE method, we can successfully allow fabricating diffraction limited X-ray mirrors for $\lambda$ as short as 0.05 nm.

## 5. Discussion

**Limitation of the 2D-IBF system**. More IBF runs can be potentially performed to further improve the figure accuracy in both the CAs shown in Figs. 9 and 10, however, the BRF used in the experiment is so powerful (7 nm/s peak removal rate) that the total dwell time at each individual dwell point is less than 0.05 s. The acceleration of our translation stage can hardly match this rapid motion requirement. Also, to correct the BRF footprints, a smaller BRF is required. We are designing and fabricating 1 mm and 2.5 mm pinholes to constrain the size and power of the BRF, with which we are going to push the figure errors to less than 0.5 nm RMS level in future work.

**Generalization of the RISE concept**. It is also worth mentioning that, although RISE is explained and tested on polishing rectangular surfaces using IBF, it is generally applicable to any CCOS techniques and arbitrary surface shapes. For example, Zernike polynomials can be used in the polynomial fitting steps in RISE for circular surfaces. For any other special shaped surfaces, RISE can be applied based on the known measurement data and either a rectangular or a circular polynomial basis can be employed.

## 6. Conclusion

In this study, a Robust Iterative Surface Extension (RISE) method was proposed to obtain reliable dwell time solutions in sub-nanometer X-ray mirror fabrication. RISE employs a polynomial-based extension strategy to guarantee the continuity at the extension boundary as well as ensure that only the correctable figure errors are extended. An iterative scheme is then applied to progressively refine the dwell time solution so that the desired figure errors can be achieved. The RISE concept has been applied to the fabrication of two X-ray diffraction limited elliptic cylindrical surfaces on the same spherical mirror using ion beam figuring. Both clear apertures reached sub-nanometer figure errors, which proves the effectiveness of the proposed RISE method.

## Funding

Brookhaven National Laboratory (BNL LDRD 17-016); Office of Science (DE-SC0012704).

## Acknowledgment

The authors would like to sincerely thank Richie Nethery from Kaufman & Robinson for his strong support in helping us tune and stabilize the ion source. Special thanks go to Dr. Konstantine Kaznatcheev from NSLS-II for his strong support and close collaboration in performing the AFM measurement for the experiment.

## Disclosures

The authors declare no conflicts of interest.

## Data availability

Data underlying the results presented in this paper are not publicly available at this time but may be obtained from the authors upon reasonable request.

## References

**1. **R. A. Jones, “Computer-Controlled Polishing Of Telescope Mirror Segments,” in Advanced Technology Optical Telescopes I, vol. 0332L. D. Barr and G. Burbidge, eds., International Society for Optics and Photonics (SPIE, 1982), pp. 352–356.

**2. **J. Fanson, R. Bernstein, G. Angeli, D. Ashby, B. Bigelow, G. Brossus, A. Bouchez, W. Burgett, A. Contos, R. Demers, F. Figueroa, B. Fischer, F. Groark, R. Laskin, R. Millan-Gabet, M. Pi, and N. Wheeler, “Overview and status of the giant magellan telescope project,” in Ground-based and Airborne Telescopes VIII, vol. 11445 (International Society for Optics and Photonics, 2020), p. 114451F.

**3. **M. Ghigo, G. Vecchi, S. Basso, O. Citterio, M. Civitani, E. Mattaini, G. Pareschi, and G. Sironi, “Ion figuring of large prototype mirror segments for the E-ELT,” in Advances in Optical and Mechanical Technologies for Telescopes and Instrumentation, vol. 9151R. Navarro, C. R. Cunningham, and A. A. Barto, eds., International Society for Optics and Photonics (SPIE, 2014), pp. 225–236.

**4. **M. Weiser, “Ion beam figuring for lithography optics,” Nucl. Instrum. Methods Phys. Res., Sect. B **267**(8-9), 1390–1393 (2009). [CrossRef]

**5. **L. Wischmeier, P. Graeupner, P. Kuerz, W. Kaiser, J. van Schoot, J. Mallmann, J. de Pee, and J. Stoeldraijer, “High-na euv lithography optics becomes reality,” in Extreme Ultraviolet (EUV) Lithography XI, vol. 11323 (International Society for Optics and Photonics, 2020), p. 1132308.

**6. **A. Schindler, T. Haensel, A. Nickel, H.-J. Thomas, H. Lammert, and F. Siewert, “Finishing procedure for high-performance synchrotron optics,” in Optical Manufacturing and Testing V, vol. 5180 (International Society for Optics and Photonics, 2003), pp. 64–72.

**7. **H. Thiess, H. Lasser, and F. Siewert, “Fabrication of x-ray mirrors for synchrotron applications,” Nucl. Instrum. Methods Phys. Res., Sect. A **616**(2-3), 157–161 (2010). [CrossRef]

**8. **L. Peverini, I. Kozhevnikov, A. Rommeveaux, P. Vaerenbergh, L. Claustre, S. Guillet, J.-Y. Massonnat, E. Ziegler, and J. Susini, “Ion beam profiling of aspherical x-ray mirrors,” Nucl. Instrum. Methods Phys. Res., Sect. A **616**(2-3), 115–118 (2010). [CrossRef]

**9. **T. Wang, L. Huang, Y. Zhu, M. Vescovi, D. Khune, H. Kang, H. Choi, D. W. Kim, K. Tayabaly, N. Bouet, and M. Idir, “Development of a position–velocity–time-modulated two-dimensional ion beam figuring system for synchrotron x-ray mirror fabrication,” Appl. Opt. **59**(11), 3306–3314 (2020). [CrossRef]

**10. **T. Wang, L. Huang, H. Kang, H. Choi, D. W. Kim, K. Tayabaly, and M. Idir, “RIFTA: A Robust Iterative Fourier Transform-based dwell time Algorithm for ultra-precision ion beam figuring of synchrotron mirrors,” Sci. Rep. **10**(1), 8135 (2020). [CrossRef]

**11. **K. Yamauchi, H. Mimura, T. Kimura, H. Yumoto, S. Handa, S. Matsuyama, K. Arima, Y. Sano, K. Yamamura, K. Inagaki, H. Nakamori, J. Kim, K. Tamasaku, Y. Nishino, M. Yabashi, and T. Ishikawa, “Single-nanometer focusing of hard x-rays by kirkpatrick–baez mirrors,” J. Phys.: Condens. Matter **23**(39), 394206 (2011). [CrossRef]

**12. **K. Yamauchi, H. Mimura, S. Matsuyama, H. Yumoto, T. Kimura, Y. Takahashi, K. Tamasaku, and T. Ishikawa, “Focusing mirror for coherent hard x-rays,” Synchrotron Light Sources and Free-Electron Lasers: Accelerator Physics, Instrumentation and Science Applications pp. 1093–1122 (2020).

**13. **M. Born and E. Wolf, “Electromagnetic theory of propagation, interference and diffraction of light,” Principles of Optics pp. 370–458 (1980).

**14. **D. W. Kim, S.-W. Kim, and J. H. Burge, “Non-sequential optimization technique for a computer controlled optical surfacing process using multiple tool influence functions,” Opt. Express **17**(24), 21850–21866 (2009). [CrossRef]

**15. **C. Wang, Z. Wang, Q. Wang, X. Ke, B. Zhong, Y. Guo, and Q. Xu, “Improved semirigid bonnet tool for high-efficiency polishing on large aspheric optics,” Int. J. Adv. Manuf. Technol. **88**(5-8), 1607–1617 (2017). [CrossRef]

**16. **D. Golini, S. D. Jacobs, W. I. Kordonski, and P. Dumas, “Precision optics fabrication using magnetorheological finishing,” in Advanced Materials for Optics and Precision Structures: A Critical Review, vol. 10289 (International Society for Optics and Photonics, 1997), p. 102890H.

**17. **M. Demmler, M. Zeuner, F. Allenstein, T. Dunger, M. Nestler, and S. Kiontke, “Ion beam figuring (IBF) for high precision optics,” in Advanced Fabrication Technologies for Micro/Nano Optics and Photonics III, vol. 7591 (2010), p. 75910Y.

**18. **M. Idir, L. Huang, N. Bouet, K. Kaznatcheev, M. Vescovi, K. Lauer, R. Conley, K. Rennie, J. Kahn, R. Nethery, and L. Zhou, “A one-dimensional ion beam figuring system for x-ray mirror fabrication,” Rev. Sci. Instrum. **86**(10), 105120 (2015). [CrossRef]

**19. **L. Zhou, Y. Dai, X. Xie, X. Peng, F. Shi, and S. Li, “Methods to extend surface error map in dwell time algorithm,” in EUSPEN’s 16th International Conference & Exhibition, (2016).

**20. **L. Zhou, L. Huang, N. Bouet, K. Kaznatcheev, M. Vescovi, Y. Dai, S. Li, and M. Idir, “New figuring model based on surface slope profile for grazing-incidence reflective optics,” J. Synchrotron Radiat. **23**(5), 1087–1090 (2016). [CrossRef]

**21. **T. Wang, L. Huang, M. Vescovi, D. Kuhne, K. Tayabaly, N. Bouet, and M. Idir, “Study on an effective one-dimensional ion-beam figuring method,” Opt. Express **27**(11), 15368–15381 (2019). [CrossRef]

**22. **R. Rückriem, M. Zeuner, and R. Köhler, “P2. 2-fabrication of high-sensitivity pyroelectric sensors by ion beam etching,” Tagungsband pp. 551–554 (2016).

**23. **L. Zhou, M. Idir, N. Bouet, K. Kaznatcheev, L. Huang, M. Vescovi, Y. Dai, and S. Li, “One-dimensional ion-beam figuring for grazing-incidence reflective optics,” J. Synchrotron Radiat. **23**(1), 182–186 (2016). [CrossRef]

**24. **C. L. Carnal, C. M. Egert, and K. W. Hylton, “Advanced matrix-based algorithm for ion-beam milling of optical components,” in Current Developments in Optical Design and Optical Engineering II, vol. 1752 (1992), pp. 54–63.

**25. **C. Jiao, S. Li, and X. Xie, “Algorithm for ion beam figuring of low-gradient mirrors,” Appl. Opt. **48**(21), 4090–4096 (2009). [CrossRef]

**26. **L. Huang, T. Wang, K. Tayabaly, D. Kuhne, W. Xu, W. Xu, M. Vescovi, and M. Idir, “Stitching interferometry for synchrotron mirror metrology at National Synchrotron Light Source II (NSLS-II),” Opt. Lasers Eng. **124**, 105795 (2020). [CrossRef]

**27. **L. Huang, T. Wang, J. Nicolas, A. Vivo, F. Polack, M. Thomasset, C. Zuo, K. Tayabaly, D. W. Kim, and M. Idir, “Two-dimensional stitching interferometry for self-calibration of high-order additive systematic errors,” Opt. Express **27**(19), 26940–26956 (2019). [CrossRef]

**28. **R. Gerchberg, “Super-resolution through error energy reduction,” Optica Acta: Int. J. Opt. **21**(9), 709–720 (1974). [CrossRef]

**29. **R. J. Marks, “Gerchberg’s extrapolation algorithm in two dimensions,” Appl. Opt. **20**(10), 1815–1820 (1981). [CrossRef]

**30. **S. Wilson and J. McNeil, “Neutral ion beam figuring of large optical surfaces,” in Current Developments in Optical Engineering II, vol. 818 (1987), pp. 320–325.

**31. **Y. Liu, H. Cheng, Z. Dong, and H.-Y. Tam, “Edge effect of optical surfacing process with different data extension algorithms,” Front. Optoelectron. **7**(1), 77–83 (2014). [CrossRef]

**32. **T. Caixue, Y. Hao, L. Zijian, Z. Yuanhang, and W. Shenglin, “High precision edge extrapolation technique in continuous phase plate magnetorheological polishing,” Infrared and Laser Eng. **48**(4), 442001 (2019). [CrossRef]

**33. **X. Qian, Z. Ma, Y. Yao, and L. Shen, “Investigation of dwell time based on lucy-richardson algorithm and gercherg surface continuation algorithm,” in AOPC 2020: Optics Ultra Precision Manufacturing and Testing, vol. 11568 (International Society for Optics and Photonics, 2020), p. 1156814.

**34. **L. Yun, D. Guping, and X. Tingwen, “Smoothly extending algorithm for the surface error of optics,” Infrared and Laser Eng. **42**, 408–412 (2013).

**35. **I. Amidror, “Scattered data interpolation methods for electronic imaging systems: a survey,” J. Electron. Imaging **11**(2), 157–176 (2002). [CrossRef]

**36. **B. Yang, X. Xie, F. Li, and L. Zhou, “Edge effect correction using ion beam figuring,” Appl. Opt. **56**(32), 8950–8958 (2017). [CrossRef]

**37. **J. Wang and D. E. Silva, “Wave-front interpretation with zernike polynomials,” Appl. Opt. **19**(9), 1510–1518 (1980). [CrossRef]

**38. **V. N. Mahajan, “Orthonormal aberration polynomials for anamorphic optical imaging systems with rectangular pupils,” Appl. Opt. **49**(36), 6924–6929 (2010). [CrossRef]

**39. **F. Liu, B. M. Robinson, P. Reardon, and J. M. Geary, “Analyzing optics test data on rectangular apertures using 2-d chebyshev polynomials,” Opt. Eng. **50**(4), 043609 (2011). [CrossRef]

**40. **M. S. Rakitin, O. Chubar, P. Moeller, R. Nagler, and D. L. Bruhwiler, “Sirepo: a web-based interface for physical optics simulations - its deployment and use at NSLS-II,” in Advances in Computational Methods for X-Ray Optics IV, vol. 10388O. Chubar and K. Sawhney, eds., International Society for Optics and Photonics (SPIE, 2017), pp. 158–168.

**41. **R. Celestre, O. Chubar, T. Roth, M. S. del Rio, and R. Barrett, “Recent developments in x-ray lens modelling with SRW,” in Advances in Computational Methods for X-Ray Optics V, vol. 11493 (International Society for Optics and Photonics, 2020), p. 114930J.

**42. **V. N. Mahajan, “Strehl ratio for primary aberrations in terms of their aberration variance,” J. Opt. Soc. Am. **73**(6), 860–861 (1983). [CrossRef]