Memory efficient constrained optimization of scanning-beam lithography

: This article describes a memory efficient method for solving large-scale optimization problems that arise when planning scanning-beam lithography processes. These processes require the identification of an exposure pattern that minimizes the difference between a desired and predicted output image, subject to constraints. The number of free variables is equal to the number of pixels, which can be on the order of millions or billions in practical applications. The proposed method splits the problem domain into a number of smaller overlapping subdomains with constrained boundary conditions, which are then solved sequentially using a constrained gradient search method (L-BFGS-B). Computational time is reduced by exploiting natural sparsity in the problem and employing the fast Fourier transform for efficient gradient calculation. When it comes to the trade-off between memory usage and computational time we can make a different trade-off compared to previous methods, where the required memory is reduced by approximately the number of subdomains at the cost of more computations. In an example problem with 30 million variables, the proposed method reduces memory requirements by 67% but increases computation time by 27%. Variations of the proposed method are expected to find applications in the planning of processes such as scanning laser lithography, scanning electron beam lithography, and focused ion beam deposition, for example.


Introduction
During integrated circuit fabrication, materials are selectively added or subtracted by depositing a layer of resist material, then modifying certain areas using a lithography process [1,2].The lithography process involves exposing the resist to light through a complex mask or reticle.The exposed resist is then removed (or retained) through a development process.Many mask sets are required to produce a circuit, which can be prohibitively expensive for low-volume or prototype applications.
Although the exposure source and physical mechanism of serial maskless lithography methods are varied, many can be described as a controllable source with a known spatial distribution.For example, the spatial distribution of scanning optical and electron beam lithography is limited by the wavelength and numerical aperture.Since the dimension of the desired features are similar in scale to the spatial distribution of the beam, an optimization problem arises.In scanning beam lithography, a set of exposure locations and beam power settings must be found that optimize the resolution of developed features [27][28][29].A similar planning step is required by any method where the spatial distribution of the source cannot be neglected.
The first methods for exposure planning used rules [30,31] that were similar to those employed in projection lithography.These were improved by linear programming methods [29] that were commercialized for proximity correction [32][33][34].In 2017, quadratic cost functions were introduced to optimize feature geometry with regularization of exposed power [35].A non-linear programming approach with interior-point optimization was described in [36].Although this method was quick to converge due to the calculation of first and second order derivatives, it required the storage of an N 2 × N 2 matrix for an N × N problem, which is only suitable for small problems (300 × 300).The numerical efficiency was later improved by exploiting sparsity and approximating the first derivative [4,35]; however, the required memory was still on the order of N 2 × N 2 .An alternative to optimization is deconvolution based on the Fredholm integral [37], which reduces the largest matrix to N × N; however, this method requires an order of magnitude more iterations than gradient based methods and does not minimize a known cost function so optimality is neither guaranteed nor expected.These three methods are directly compared in [38].In summary, methods for planning serial maskless lithography processes are either limited by memory or do not guarantee optimality and require significant iterations.
Although this work focuses on maskless lithography, it should be noted that exposure optimization problems also exists in projection lithography [39][40][41][42].However, these methods are fundamentally different since the optimization variables are the binary mask and source pattern [43][44][45].Due to the large scale of mask optimization problems, considerable efforts have been made for improving the numerical efficiency; for example, through the use of basis functions [46], Lagrangian methods [47], set-based methods [48,49], and neural networks [50][51][52].
The contribution of this article is to reduce the memory requirements of gradient-based optimization methods for scanning-beam lithography planning.The proposed method splits the problem domain into a number of smaller overlapping subdomains with constrained boundary conditions, which are then solved sequentially using a constrained gradient search method (L-BFGS-B).Computational time is reduced by exploiting natural sparsity in the problem and employing the fast Fourier transform for efficient gradient calculation.Compared to previous methods [4,[35][36][37], the proposed method is slower but the required memory is reduced by approximately the number of subdomains.
The following three sections describe the forward process model, define the optimization problem, and derive the analytical gradients of the cost function.The problem space is them subdivided in Section 5, followed by memory efficient optimization in Section 6, and an example application in Section 7. A complexity analysis and an investigation of accuracy is presented in Section 8, then the article is concluded in Section 9.

Process model
This section describes a general model of scanning-beam lithography processes.Given an input exposure pattern, the model predicts the photoresist conversion fraction, which is directly related to developed features.The model assumes that the photoresist is sufficiently thin so that the beam profile remains approximately constant throughout the depth.The optical properties of the film, which are a function of the exposure state, are also assumed to be constant.Other optical effects, such as scattering and cavity formation, are ignored.
An introduction to the scanning-beam lithography process is illustrated in Fig. 1.This figure shows how a set of discrete exposures can be used to create features with a resolution comparable to the beam width.However, the middle row demonstrates that simple choices of exposure energy and location result in sub-optimal developed features.The bottom row illustrates how an optimized exposure pattern yields the best possible feature resolution, and it is also a numerical example of the proposed method.7)) which models the development process.In the middle row, a two dimensional exposure problem is considered.The 'naive' exposure pattern is a scaled version of the desired feature; however, the beam kernel B can be observed to significantly expand the desired feature due to over exposure.The middle and bottom row use a realistic development function f Z (•) plotted in Fig. 2. To minimize the achievable difference between the desired and developed features, the bottom row shows an example of an optimized exposure pattern, which results in the best possible developed feature given the beam kernel and development model.
The exposure pattern W, beam kernel B, dosage D, and predicted feature Ẑ are matrices which contain values at discrete locations in a workspace.The workspace is represented by a uniformly spaced rectangular grid where i denote the grid resolution in direction p, whereas N = N 1 N 2 is the grid size (number of grid points).Each grid point is associated with a desired feature value stored in the binary matrix The grid is exposed to a scanned beam with intensity modelled by an exponential function, which represents a two-dimensional Gaussian beam with arbitrary x and y axis width and rotation, as described in [35].That is, We refer to H = [h ij ] as the bandwidth, with entries defined in terms of its inverse according to In this work, the UV scanning laser lithography system described in [35] will be used as an example.The beam parameters for this system were measured to be w x0 = 570 nm, w y0 = 560 nm and ϕ = 2.2 • [35].
The dosage D represents the absorbed energy per unit area, which is determined by the exposure pattern W and beam kernel B. The exposure W can represent any variable that is proportional to energy, e.g.beam power for a constant exposure time, or exposure time with a constant beam power.The dosages and exposure values are stored in the matrices D = [d qk ] and W = [w ij ], respectively.At a given point, the total dosage is found by summing the contributions from the entire grid where w = vec W.
As the photoresist is exposed by the beam, the increasing dosage leads to a chemical change due to photon absorption or photocatalysis.A negative photoresist becomes less sensitive to a developing agent after exposure and the developed photoresist features are similar to the dosage pattern.Conversely, a positive photoresist becomes more sensitive to the developing agent and results in subtractive features that are similar to the dosage pattern.Regardless of the photoresist polarity, the development process is dependent on the fraction of photoresist that is converted; therefore, this quantity is used to define the desired feature.To complete the process model, a function is required that maps the dosage to the predicted conversion fraction.
The simplest model of the development process is a threshold function, where the photoresist is assumed to be 100% converted above a certain dosage threshold, as depicted in the top row of Fig. 1.However, in general the conversion function can be any increasing function of dosage.In this work, the sigmoid threshold function plotted in Fig. 2 is used to model the development process [35][36][37].The sigmoid mapping function f Z (•) from the dosage to the predicted feature Ẑ where α is the steepness of the curve, and d 50% is the dosage where half of the photoresist is converted.The steepness is plotted for α = 5, 10, 20, 40 in Fig. 2. In Section 7, the proposed method is applied to the optimization of ultraviolet scanning laser lithography with a positive photoresist (AZ ECI3007, MicroChemicals, Germany).The conversion fraction of this photoresist was adequately modeled using α = 5 in Ref. [35][36][37].It is convenient to normalize the fully converted photoresist to 1 and therefore d 50% = 0.5.7) models the development process.In scanning laser lithography, this function relates the exposure energy to the fraction of converted photoresist.It is parameterized by the steepness α and the dosage d 50% where half of the photoresist is converted.The steepness is plotted for α = 5, 10, 20, 40.

Optimization problem
The optimization objective is to find an exposure pattern w ij ≥ 0 that minimizes the difference between the desired feature Z and the predicted feature Ẑ.The non-negative constraint on w ij is due to the nature of exposure energy, which cannot be less than zero.It also desirable to minimize the total dosage, since this minimizes negative effects such as heating, scattering, and other background exposure processes such as reflection to and from the substrate and objective lens.These optimization objectives are summarized by the following problem The penalty parameter γ>0 balances the trade-off between feature matching and total dosage.Penalizing the L 2 norm is preferred in this application as the most significant sources of background exposure are scatter within the substrate, and reflection between the substrate and objective.Both of these sources are proportional to cumulative power, which is proportional to the L 2 norm.In this work, the optimal value of γ was found to be 10 −4 , which is recommended as a starting point in other applications.Some experimentation may be required as the regularization is likely to be affected by other factors such as photoresist properties.

Gradient calculation
Following the procedure described in [53], the fast Fourier transform (FFT) is utilized for efficient computation of the cost function and its gradient.First, Eq. ( 6) is rewritten as a convolution where T ) and w q−i,k−j = 0 if the indexed point lies outside the grid.The dosage in Eq. ( 6) is recovered if L 1 = N 1 − 1 and L 2 = N 2 − 1; however, as the beam intensity of Eq. ( 2) decays exponentially and therefore in practice lacks support outside a finite region, Eq. ( 6) can be approximated by shrinking L 1 and L 2 for faster computations.A suitable choice is where τ is a user-defined scaling parameter, λ is the maximum eigenvalue of H and ⌈•⌉ denotes the ceiling function.The authors of [53] suggest using τ = 3.7 in the context of bivariate kernel density estimation; however, to safe-guard against numerical inaccuracies, a slightly higher value of τ = 10 is used here.The gradient of the cost function in Eq. ( 8) can also be expressed as a convolution [35].Specifically, ∇f = vec G where G = [g qk ] and Before applying the FFT, adequate zero-padding is required as the bandwidth H is unconstrained (non-diagonal) [53].The matrices are and where the entry w 11 is placed at index (L 1 + 1, L 2 + 1).The matrices Z, Ẑ and D are formed with the same structure as in W. The 0-blocks are built to match the dimension P 1 × P 2 , where P 1 and P 2 are chosen as a suitable composite integers.Next, the dosage and gradient are computed where F and F −1 denote the FFT and its inverse, respectively, while • denotes the Hadamard (element-wise) product of matrices.The desired matrices D and G are found in D and G as the blocks defined by the rows 2L 1 + 1 to 2L 1 + 1 + N 1 and columns 2L 2 + 1 to 2L 2 + 1 + N 2 .

Subdomain division
For large problems sizes when the number of variables N renders the problem computationally infeasible, the optimization can be performed in different subdomains.The total solution is then obtained by concatenation of the local solutions.
Let w * denote the optimal solution on the grid X, and let w Π * ∈ w * denote the optimal solution on the subdomain Π ∈ X.The limited support of the beam intensity (Eq.( 2)) assures that a variable w ij affects the value of f only in a limited region.That being said, there is a possibility of a chain effect propagating through the grid since nearby variables are affecting each other.However, the binary structure of Z significantly reduces this effect.Hence, by solving the problem on the subdomain, an estimate ŵΠ * of w Π * is obtained.However, solving the problem only on the domain Π is insufficient as the boundary components may be incorrect.The trade-off is to use another region Γ such that Π ∈ Γ ∈ X.Then ŵΠ * can be obtained by extracting the correct components from the solution ŵΓ * computed on Γ.Thus, if Γ is sufficiently larger than Π, the difference between ŵΠ * and w Π * is negligible.In practice, for the class of problems under consideration, it was found to be reasonable to form Γ by adding L p points to Π in direction p.A simple illustration of this subdomain division is given in Fig. 3.Although this approach is novel for the application under consideration, it should be stressed that subdomain division and local optimization do not constitute a novel idea in general.In a broader context, so called domain decomposition methods [54,55] are well-established tools for solving differential equations through subdomain division, and originate from the work by Hermann Schwarz in 1870 [54].Due to the iterative procedure of these methods, they are in essence comprising local optimization problems.Example applications include finite element solvers [56], total variation de-noising [57] and plasticity with hardening [58].Moreover, the same principle has for example been used for the ptychography problem [59], including efficient parallel GPU implementations [60,61].In [62] a similar approach is used to reduce the computational time within composite pile foundation.
It can observed that the subproblems are independent of each other.In principle, this makes the optimization approach straight-forward to solve in parallel by distributing the subproblems on different CPU cores.However, this conflicts with the aim of reducing the memory requirements since the different processes compete for the same RAM memory.In the case where the RAM memory is sufficiently large, solving the subproblems in parallel will reduce the computation time by approximately the number of subdomains.Future improvements to computation time could also be achieved by utilizing a GPU for FFT computation.

Memory efficient optimization
The limited-memory Broyden-Fletcher-Goldfarb-Shanno algorithm with bound constraints (L-BFGS-B) [63] is employed to solve the optimization problem of Eq. ( 8).This relies upon the BFGS method [64], which is a well-established quasi-Newton algorithm that approximates the inverse Hessian matrix and provides significantly faster convergence.
The limited-memory BFGS (L-BFGS) [65] method is a modification of BFGS that avoids building an entire matrix, which is beneficial for high-dimensional problems.The L-BFGS method expresses the BFGS update recursively and disregards all but the m latest Hessian approximations.That way the computation of the search direction can be carried out in memory-efficient for-loops.
The L-BFGS-B method [63] is a further extension that allows for bounded input constraints of the form l ij <w ij <u ij .At each iteration k, it minimizes a quadratic model of the cost function along the path In the problem under consideration, l ij = 0 and u ij = ∞ for all inputs w ij .The optimization is initialized by setting w ij = 1 if z ij = 1, which corresponds to all grid points that lie inside the desired feature pattern.The remaining variables are initialized to 0. Judging from previous experiments [35][36][37] and by intuition, grid points outside the pattern are unlikely to be exposed at the optimum.
To speed up the optimization, the sparsity of the problem is exploited.The finite region of the desired feature combined with the limited support of the kernel function allows for an a priori identification of grid points that will, in practice, be un-exposed at the optimum.For instance, variables at grid points that lie outside the desired pattern can be disregarded.Hence, the corresponding entries are removed from the gradient before it is returned to the L-BFGS-B routine.
In the implementation, we are using the MATLAB wrapper [66] of the L-BFGS-B implementation described in [67].The proposed procedure is summarized in Algorithm 1 and illustrated in Fig. 4.
Note that the L-BFGS-B optimizer is not a crucial choice.There are several other optimizers suitable for the large-scale problems, such as gradient descent as used in [35].The constraint w ij ≥ 0 can also be achieved through a variable transformation rather than by the optimizer, for instance by setting w ij = e p ij and then optimizing with respect to p ij .A pure gradient based routine is more memory efficient than L-BFGS-B, which keeps track of the m latest pairs of input points and gradients.The experiment below uses m = 5, which adds a relatively small memory requirement compared to the matrix formations and operations involved when computing the cost function and its gradient.However, a smaller value of m can be used as well (even m = 1, which resembles a gradient descent optimizer).L-BFGS-B is used in this work as it is a well-established method with an efficient and readily available implementation, and it is appealing due to the Hessian approximation and the built in support for constrained variables.Finally, L-BFGS-B worked well in practice for a range of problems.Assuming that the problem is too large to be solved directly, the grid X is split horizontally into two equally-sized subdomains Π 1 and Π 2 (middle).To avoid undesired boundary effects in the solutions, the extended regions Γ 1 and Γ 2 are constructed, of which one is shown above (right).The yellow parts indicate variables that with great certainty will be zero at the optimum and therefore are excluded from the optimization.Hence, we are interested in the variables indicated by the beige overlay.These are then concatenated with the equivalent part of Γ 2 to form the solution over X.

Numerical example
This section demonstrates the proposed method on a medium scale lithography test pattern shown in Fig. 5(a).The test pattern represents a 100 µm × 100 µm area with 5486 × 5488 pixels and a resolution of 18.2 nm; that is, there are approximately 30 million variables to optimize.The beam parameters that represent the laser and optical system were experimentally identified in previous work as w x0 = 570 nm, w y0 = 560 nm and ϕ = 2.2 • [35].The photoresist fully exposed state is defined to be 1 and d 50% = 0.5.The steepness of the photoresist exposure curve (Eq.( 7)) is α = 5, and the penalty parameter γ = 10 −4 .
The test pattern shown in Fig. 5(a) examines the limits of a lithography process and is ideal for evaluating two important aspects of the optimization procedure.First, the optimization should result in consistent lines and features that do not vary along the length, or at different locations in the pattern, unless there are nearby features that require consideration.In other words, there should be no visible artefacts or asymmetry in geometric features.Secondly, the optimized features should degrade predictably as the objective becomes unachievable.The minimum feature size in scanning beam lithography is determined by the steepness of the photoresist exposure curve and the beam width.In this example where α = 5, the minimum size of an exposed feature is approximately equal to the beam width, which is shown for reference in the close-up views in Fig. 5.
The optimization procedure was completed for two cases, one where the entire pattern was optimized in one step, and another where the area is split into a 2 × 2 subdomain array.Both optimizations resulted in numerically identical solutions but the 2 × 2 subdomain required 67% less memory at the expense of 27% longer optimization time.The longer optimization time is due to additional computations required by the overlap area.The overlap area is also the reason why memory usuage is not reduced by 75%.The numerical results are summarized in Table 1.
Since the optimization results were numerically identical, only the results for the 2 × 2 solution are reported in Fig. 5.The top row of Fig. 5 shows the initial guess of the exposure pattern (Fig. 5(a)), which results in the dosage (Fig. 5(b)) and the resulting feature (Fig. 5(c)), which is grossly over exposed.The pixels of the initial guess (Fig. 5(a)) are set to 1 where the desired output feature is non-zero; therefore, (Fig. 5(a)) represents both the initial guess and the desired output feature.The middle row shows the optimized exposure pattern (Fig. 5(d)) and the resulting dosage (Fig. 5(e)) and output feature (Fig. 5(f)).In the close up view of the corner features and the resulting dosage (e) and output feature (f).In the close up view of the corner features (g), the optimization result is observed to degrade predictably as the line width and spacing approach the beam-width, which is illustrated by a white bar.In the close up views (h) and (i), the simulated exposure is observed to produce parallel lines and geometric features with no visible artefacts and degrades predictably as the feature size approaches the beam-width.Form Γ i such that Π i ∈ Γ i Extract ŵΠ i * ∈ ŵΓ i * 8: end for 9: Concatenate the solutions { ŵΠ i * } to form ŵ * (Fig. 5(g)), the optimization result can be observed to degrade predictably as the line width and spacing approach the beam width, which is illustrated by a white bar.In the close up views (Fig. 5(h) and 5(i)), the simulated exposure is observed to produce parallel lines and geometric features with predictable degradation as the feature size approaches the beam width.In summary, the proposed subdomain method optimizes an exposure pattern with approximately 30 million variables subject to positivity constraints in 30 minutes on a modest PC.Memory requirements were reduced by nearly three-quarters compared to a direct solution using the same underlying method.This example demonstrates the utility of the proposed method for increasing the scale of scanning beam lithography optimization, or significantly decreasing the memory requirements.

Complexity and accuracy analysis
As the memory requirement is directly proportional to the grid size, it scales with O(N).The complexity of the cost function and gradient computations are dictated by the convolutional operations, which are of order O(N log N) when using the FFT.The approximative FFT computations should improve this scaling, as the reduced convolutional limits in Eq. ( 10) are independent of N. As the L-BFGS-B routine has a computational cost of approximately O(N) per iteration [63], the overall complexity is expected to be (at most) of order O(N log N).
To confirm the above estimates, an empirical study was performed on the example problem described in Section 7. The results of the analysis are plotted in Fig. 6.In Fig. 6(a), the computation time of the cost function and gradient is plotted against the number of variables, averaged over 10 trials.The plot includes results for exact FFT computations (τ = ∞) and the approximate version with τ = 10.Reference lines corresponding to scaling rates of O(N log N) and O(N) are also plotted.The results show that computation time scales with O(N) rather than O(N log N), especially when τ = 10.Furthermore, the approximation reduces the computation time with a factor of 4-6 on this interval (increasing with N).
The time complexity of the complete optimization procedure (here comprising 50 iterations) is plotted against N in Fig. 6(b).The results also show a scaling of rate O(N), which confirms that the L-BFGS-B routine does not degrade the overall complexity.In addition to L-BFGS-B with memory limit m = 5, this plot includes two versions of gradient descent (GD) obtained by setting m = 1: one using the same bound constraint treatment as in L-BFGS-B, and the other using the variable transformation w ij = e p ij mentioned in Section 6.It is seen that the time consumption per iteration is lower in these cases (roughly 40% reduction).Also, Fig. 6(c) shows the evaluation of the cost function for the largest number of variables (N = 10 8 ) used in the experiments of Fig. 6(b).The bound constrained versions converge faster initially.However, the results are very similar for the memory limits 5 and 1 -as noted in [64], the optimal choice of memory limit is highly problem dependent.A study was also performed on the effect of subdomain division on the accuracy of results and the computation time.Figure 7 plots RMS difference between the desired and predicted feature, and the computation time for different subdomain divisions.The results are reported relative to the full grid solution with exact FFT computations (1 × 1, τ = ∞).It is seen that the accuracy is not affected by the FFT approximation, or by the subdomain division (maximum error less than 0.02%).Although there is a computational overhead from the subdomain division, it is small in comparison to the time we gain from the approximation.These studies were conducted on a Scientific Linux 6.10 server with an AMD Opteron (Bulldozer) 6282SE processor and 128 GB RAM.

Conclusion
The article extends previous work on scanning beam lithography processes by proposing a memory efficient optimization procedure.The benefits of reduced memory requirements include lower minimum hardware specifications, or direct solution of larger scale problems, or a finer pixel resolution.
Instead of optimizing for the full grid directly, the problem is divided into smaller subdomains that are partly overlapping to avoid inaccurate boundary effects.The sub-problems are solved one at a time using the L-BFGS-B gradient search algorithm, and are then concatenated into a global solution.Furthermore, efficient computations are obtained through local application of the fast Fourier transform and utilization of natural sparsity.The overlaps between the subdomains imply a computational overhead resulting in a longer run-time; the upside is that the memory requirements are reduced by a factor almost equal to the number of subdomains, which opens up for problem sizes much larger than otherwise possible.
Empirical experiments demonstrate the performance of the proposed method; for instance, by splitting a problem with 30 million variables into 4 subdomains, it is solved with 67% less memory and a time increase of merely 27%.This trade-off between reduced memory and longer computation time is expected to be desirable in large-scale applications that tend to be memory limited.Moreover, it is shown that the time complexity scales linearly with the problem size and that the subdomain division has a negligible impact on the accuracy of the solution.
This work focuses on applications in scanning laser lithography.Future work will consider modified process models and cost functions that are suited to other point-wise lithography and direct fabrication processes such as scanning electron beam lithography, focused ion beam deposition, and plasma etching.

Fig. 1 .
Fig. 1.Illustration of the scanning-beam lithography process.The top row shows a one-dimensional interpretation, where a discrete exposure pattern W represents the beam exposure time or power setting (vertical axis) at discrete locations (horizontal axis).The resulting dosage D (exposure energy per unit area) is the sum of beam kernels B, scaled and shifted by the exposure pattern W. The output feature Ẑ is predicted by the threshold function f Z (•) (Eq.(7)) which models the development process.In the middle row, a two dimensional exposure problem is considered.The 'naive' exposure pattern is a scaled version of the desired feature; however, the beam kernel B can be observed to significantly expand the desired feature due to over exposure.The middle and bottom row use a realistic development function f Z (•) plotted in Fig.2.To minimize the achievable difference between the desired and developed features, the bottom row shows an example of an optimized exposure pattern, which results in the best possible developed feature given the beam kernel and development model.

Fig. 3 .
Fig. 3. Illustration of subdomain division.The problem is solved on the Γ-regions, and then the Π-regions are concatenated.

Fig. 4 .
Fig. 4. Illustration of Algorithm 1.The desired feature is the C-shaped black region (left).Assuming that the problem is too large to be solved directly, the grid X is split horizontally into two equally-sized subdomains Π 1 and Π 2 (middle).To avoid undesired boundary effects in the solutions, the extended regions Γ 1 and Γ 2 are constructed, of which one is shown above (right).The yellow parts indicate variables that with great certainty will be zero at the optimum and therefore are excluded from the optimization.Hence, we are interested in the variables indicated by the beige overlay.These are then concatenated with the equivalent part of Γ 2 to form the solution over X.

Fig. 5 .
Fig. 5. Example application of the proposed method to a lithography test pattern.The top row shows the desired output feature (a) which is also used as the initial guess of the exposure pattern.The exposure (a) results in the dosage (b) and the resulting feature (c), which is grossly over exposed.The middle row shows the optimized exposure pattern (d)and the resulting dosage (e) and output feature (f).In the close up view of the corner features (g), the optimization result is observed to degrade predictably as the line width and spacing approach the beam-width, which is illustrated by a white bar.In the close up views (h) and (i), the simulated exposure is observed to produce parallel lines and geometric features with no visible artefacts and degrades predictably as the feature size approaches the beam-width.

5 : 6 :
Identify sparse entries in Γ i and keep them constant 0 Compute ŵΓ i * = argmin w f (w) on Γ i using L-BFGS-B 7:

Fig. 7 .
Fig. 7.The impact of subdomain division on accuracy and computation time.The results reported are relative to the full grid with exact FFT computations (1 × 1, τ = ∞).

Table 1 . Optimization time and memory usage for the example in Fig. 5. When the area is split into 4 overlapping subdomains, the memory usage is reduced by 67% at the expense of a 27% increase in optimization time. The optimization was implemented in MATLAB on a Windows 10 PC with an i5-6300HQ processor and 32 GB RAM.
Exposure optimization 1: Input: Grid X, desired feature Z 2: Divide X into n ≥ 1 subdomains {Π i } 3: for i = 1 . . .n do 4: