1. INTRODUCTION

Diffractive optical elements (DOEs) are flat optical components that modulate light through surface relief patterns etched onto substrates, with feature sizes comparable to the operating wavelength. The high degree of design freedom and the ability to harness diffractive effects provide DOEs with inherent versatility and distinct advantages over traditional bulk optics, driving their growing demand in various applications such as beam shaping [1], beam splitting [2], and extended-depth-of-focus imaging [3]. When combined with algorithms for compact and versatile optical signal encoding, advanced applications become feasible, including seeing through occlusions [4], lensless microscopy [5], and optical computing in the spatial domain [6]. The primary fabrication method for diffractive optical elements (DOEs) is multi-step photolithography. Despite its precision, scalability, and reproducibility [7], this approach is limited to producing DOEs with a small number of discrete phase levels, typically 4 or 8, due to the sequential nature of the photolithography steps. This restriction constrains the efficiency and range of applications of the DOEs. To achieve higher efficiency and greater versatility, continuous-relief DOEs using grayscale direct laser writing lithography (DLWL) [8,9] have become the preferred choice. Compared with other direct writing technologies, such as electron beam lithography (EBL) [10] and femtosecond direct laser writing (FLDW) [11], DLWL is particularly well-suited for mass production, offering impressive speed and cost-effectiveness [8,9,12,13]. It enables the creation of quasi-continuous profiles in a single step, with grayscale levels ranging from 256 (8-bit) to more than 4096 (12-bit), allowing precise control over microstructure dimensions and geometry. This capability is crucial for fabricating smooth-profile, multifunctional DOEs with exceptional performance.

However, DOEs fabricated by DLWL face a significant challenge from the optical proximity effect (OPE). OPE originates from the diffraction of light within the photoresist during the lithography process, leading to unexpected deviations in fabricated feature sizes. The impact of OPE is particularly prominent in DLWL, which requires the use of direct-write lenses with lower numerical apertures to achieve larger depth of focus and larger writing field areas. Conventionally, OPE can be mitigated through either DOE optical design or fabrication optimization. DOE optical design optimization means the maximization of its figures of merit characterized by its optical performance. This optimization mainly relies on iterative algorithms such as Gerchberg–Saxton [14] and direct binary search [15]. During the design process, heuristic fabrication constraints like minimum pixel size requirements can be applied to generate surface profiles with larger feature sizes and are thus easier to fabricate. There are also some more complex constraints such as continuous phase constraints [16,17]. However, these constraints inevitably lead to under-utilization of design freedom, leading to sub-optimal performance. Also, the effectiveness of such constraints varies with specific fabrication parameters and is hard to characterize. Another strategy in mitigating OPE is to modify the fabrication parameters to compensate the impact of OPE, i.e., optical proximity correction (OPC). OPC has been widely used in the lithography field for integrated circuit processing, and a series of models have been proposed. For example, the exposure and development processes of the photoresist can be described by a series of kinetic reaction equations, known as the Dill model [18–20]. To reduce complexity and enhance feasibility, empirically trained resist models have been developed, including variable-threshold models [21] and artificial neural network-based models [22–25]. Neural network-based models have demonstrated promising results in capturing complex resist behaviors and improving prediction accuracy. However, challenges persist in the further correction of OPE due to limitations in model accuracy, calibration, generalization, and the significant amount of data required to effectively train neural networks.

In this paper, we propose an “end-to-end” fabrication-integrated design (FID) approach to enable high-quality DOE design. Our method integrates the DOE optical design with fabrication optimization in an automated manner. This integrated method eliminates the need for further correction of manufacturability issues, addressing OPE and enabling comprehensive optimization by exploiting the degrees of freedom in each process to enhance overall system performance. First, a simple simulation model that characterizes the OPE and the nonlinear response of the photoresist to exposure dose is developed. This model, which comprises a low-pass filter and a nonlinear mapping, can be fully calibrated using just a single test pattern. Subsequently, a wave optics simulator generates the corresponding optical distribution based on output of the fabrication simulator. A loss function evaluates the disparity between the target and simulated optical distributions, directing the gradient-based optimization of fabrication parameters from an arbitrary starting point. The ability of our simulators to compute gradients automatically [26] streamlines the efficient optimization of extensive fabrication parameters. The simulators we utilize are intentionally crafted to be simple yet robust, requiring minimal experimental calibration and ensuring adaptability across various applications. The effectiveness of our FID methodology is validated through a series of experiments across various DOE application scenarios. Specifically, for beam splitters, we achieved more than a 25% improvement in zeroth-order suppression and uniformity across diffraction orders. Additionally, significant enhancements were observed in the quality and fidelity of digital holograms, the performance of structured light projectors, and the effectiveness of extended depth of focus (EDOF) imaging lenses. Our results validate the effectiveness and scalability of the proposed method, demonstrating its significant potential in enhancing the design and fabrication of large-scale, complex optical devices. The implications of our work extend to various domains that rely on micro–nano optics, including telecommunications, biomedical imaging, advanced manufacturing, and optical computing, promising substantial benefits and innovations across these fields.

2. DESIGN METHODOLOGY

A. Principle of Fabrication-Integrated Design

Figure 1 demonstrates our FID approach, which incorporates a differentiable forward simulation model and employs gradient descent optimization enabled by automatic differentiation. The forward model contains a physics-informed fabrication model and a diffraction optics model. The term “physics-informed” indicates that the fabrication model, while simplified compared to existing precise modeling methods, effectively captures the essential physics to provide an accurate representation of the DLWL fabrication process and is straightforward to calibrate. This simplified process is characterized by low-pass filtering, reflecting the focal spread of the writing exposure lens, and a nonlinear mapping process, representing the response to exposure levels after development. Using this forward model ?, the simulated diffraction pattern II can be predicted from any given dose distribution DD. The differentiability of the model allows for direct calculation of the gradient with respect to a loss function during each simulation, enabling gradient descent-based optimization from any initial DD to iteratively enhance the optical performance as assessed by the predefined loss function. Central to this methodology is the forward model ?, composed of two cascaded simulators: the fabrication simulation model ? and the scalar diffraction-based optical simulation model ? for optical propagation. The specific characteristics of these models are thoroughly explained in the following subsections, establishing the foundation for our method. The forward model is formulated by the sequential combination of ? and ? and is thus referred to as the fabrication-integrated design model, represented mathematically as

This model ? directly predicts the diffraction pattern or light intensity (|⋅|2, the squared magnitude of the complex amplitude) on the target surface, denoted as II, resulting from the use of a DOE fabricated with specific fabrication parameters DD. In other words, the FID model establishes the direct relation between the optical fields II and the fabrication parameters DD, expressed as IDI=?[D]. Based on this model, the dose distribution DD can be thus directly optimized to achieve the intended optical outcomes. The optimization problem is formulated as follows:

where Loss(⋅) represents a chosen loss function, such as the mean square error (MSE) between the desired optical intensity distribution and the simulated one. IItarget denotes the optical design goals of the DOE, characterized by the desired intensity distribution on the target surface, and ? is a learnable scaling factor for adaptive normalization. During the optimization process, dose distribution D values are represented as floating-point numbers, with 8-bit quantization applied only after optimization. This introduces a quantization error, which we found to be negligible compared to fabrication errors. Experimentally, fabrication errors, even under simplified conditions, were larger than this 1/256 quantization error. However, as fabrication technologies improve, quantization errors could become a dominant source of error, necessitating post-design simulations to assess their impact. The FID optimization strategy integrates fabrication and optical considerations, facilitating the creation of microstructures that are both fabrication-efficient and performance-optimized, by addressing fabrication constraints in a unified optimization step.

Fig. 1. Workflow of the fabrication-integrated design (FID) using differentiable forward modeling. Step 1: a direct laser writing lithography system projects exposure light onto a photoresist layer, creating an initial dose distribution DD, which serves as the model input. Step 2: the differentiable forward model predicts the output through two components—the lithography manufacturing model, which simulates fabrication morphology using range limitation, low-pass filtering, and nonlinear mapping, and the wave optics propagation model, which calculates the intensity distribution on the target plane. Step 3: the predicted diffraction pattern I(?,?) is compared to the target pattern Itarget(?,?) using a loss function. Step 4: automatic differentiation computes the gradient of the loss function D?Loss/?D, enabling iterative updates to the dose distribution DD′ to optimize the pattern generation.

However, Eq. (2) introduces a significant optimization challenge, particularly for centimeter-scale DOEs, where the pixel count in the pixelated fabrication field DD can reach hundreds of millions. Thus, it is essential to emphasize the critical role of making the forward model differentiable, which enables automatic differentiation. With a differentiable forward model, the calculation of the gradient with respect to the fabrication parameter DD, D?Loss/?D, becomes not only feasible but also highly efficient. This computational ease is facilitated by the utilization of automatic differentiation libraries, such as TensorFlow or PyTorch. Based on the gradient, the process of updating via gradient descent for each optimization iteration is simply encapsulated by

where ? represents the adjustable learning rate. To initiate the optimization process, an arbitrary starting dose distribution DD0, such as randomly distributed values within a valid range, can be used to facilitate a broad exploration of the solution space. The optimization procedure also benefits greatly from the application of various stochastic gradient descent (SGD) techniques, such as adaptive moment estimation (Adam) [27] to avoid local minimum stagnation. Additionally, the optimization process can be effortlessly carried out on a graphics processing unit (GPU) with improved speed. The flexibility to employ a diverse array of loss functions further enables the execution of intricate or multi-objective optimization strategies, offering a robust framework for tackling the complexities inherent in optimizing large-scale DOEs.

B. Fabrication Simulator

The writing resolution of DLWL systems is typically at the micrometer level, constrained by the point spread function (PSF) of the writing lens. In DLWL systems, the PSF can be represented as a two-dimensional (2D) function describing the spatial intensity distribution of the exposure light at each writing pixel on the focal plane of the writing lens. Given that the depth of focus of DLWL systems generally exceeds the photoresist thickness, the PSF can be considered uniform throughout the entire depth of the photoresist. Based on the 2D PSF, the modulation transfer function (MTF) is defined as the magnitude of the 2D Fourier transform of the PSF. The MTF quantifies the transfer of high-frequency information from the exposure pattern to the photoresist. An ideal DLWL system exhibits a high and flat MTF across all spatial frequencies, corresponding to a narrow PSF, indicating minimal blurring in the transfer of fine details from the exposure pattern to the photoresist.

However, DLWL systems for grayscale lithography typically employ low-NA writing lenses to achieve an extended writing depth of focus (DOF) [8] and improved fabrication productivity. This comes at the cost of a compromised spot size or PSF due to the diffraction limit. Consequently, these low-NA writing lenses can cause significant optical proximity effects (OPEs), which can be modeled through either PSF or MTF. Due to its ease of measurement, we utilize the MTF in our fabrication simulator.

The input of the fabrication simulator is denoted as DD, which represents the grayscale 2D dose distribution used for DLWL. It specifies normalized grayscale intensity values for each writing pixel. The pixel size is selected to be smaller than the resolution of the writing lens in the DLWL system. Our experiments utilized the Heidelberg DWL 66+ system equipped with a 5-mm focal-length writing lens, facilitating an 8-bit grayscale depth with dose values ? spanning 0 to 255 and an adjustable pixel pitch (0.5, 1, or 2 µm).

To maintain the dose value within the predefined range, the physics-informed fabrication model incorporates a hard limiter ?(⋅) at the initial stage, which is defined by

Subsequent to this limitation, we introduce a low-pass filtering mechanism that simulates the blurring effect induced by OPE. This degradation is characterized by a real-valued frequency-domain modulation transfer function (MTF), represented as ?(??,??). Here, ?? and ?? represent the spatial frequencies of the exposure pattern along the ? and ? axes, respectively, following the standard optical projection modeling methods [28]. For simplification, a rotationally symmetric MTF is employed for the DLWL lens system, a plausible assumption given the rotationally symmetric nature of the writing lens. This allows for the utilization of a single-variable MTF function, ?(?), where ? is the radial spatial frequency computed as ??2+??2. Following the low-pass filtering stage, a nonlinear mapping function ?(⋅) is applied to model a contrast curve of the photosensitive polymer. This curve captures the sensitivity of photosensitive polymer and characterizes its response to light exposure and the development process, establishing the relationship between a given dose and the resulting depth on photoresist. Further information regarding the MTF and contrast curve is provided in Section 2 and Supplement 1. Finally, we can thus define the physics-informed fabrication model ?(⋅), which is a nonlinear operator that estimates the fabricated surface profile hh for any given input dose distribution DD:

where ?[⋅] and ?−1[⋅] are the Fourier transform and inverse Fourier transform, respectively. ? denotes the MTF, and ° represents element-wise multiplication. Based on the expression of the physics-informed model, we can utilize two fast Fourier transforms (FFTs) to accelerate the computation of the fabrication model. Moreover, given that this model is constructed from a series of elementary operations, its implementation as a differentiable model is straightforward using automatic differentiation libraries.

C. Optical Model

The optical model of DOE has been well studied. In the case of paraxial regime and 2? phase modulation for thin DOEs, the thin element approximation [29] alongside the scalar wave propagation model [30] is usually used with desirable accuracy. The optical wave propagation operator ?(⋅) is modeled by scalar diffraction theory using the angular spectrum method.

The optical model takes the predicted surface profile hh of the DOE, obtained from the fabrication model under specific fabrication parameters, as its input. The output of the model is the complex amplitude of the light field at the target plane. Based on the thin element approximation, the modulation of the light phase introduced by the DOE is given by h2?(?−1)h/?, where ? is the refractive index of the material at operating wavelength, and ? is the operating wavelength. On the basis of scalar diffraction theory, the expression of the optical model can then be represented as

where ?? represents the transfer function characterizing free-space propagation with a distance of ?.

Utilizing the angular spectrum method, ?? is expressed as

Thus, the evaluation of the optical model also involves only twice the FFT algorithm, which facilitates a differential and efficient calculation.

D. Model Calibration and Fitting

The model parameters are essential to achieve accurate results and can be obtained through several calibration procedures. The optical model is predefined and requires no calibration. However, in the fabrication model, both the contrast curve and the MTF must be calibrated with experimental data. Figure 2 illustrates the entire calibration process.

Fig. 2. Calibration and fitting steps for the FID model. Two specific patterns are used to calibrate low and high-frequency modes. In (a), large-area block patterns with varying dose distributions are fabricated and measured using a 3D profilometer, yielding a dose-height relationship unaffected by optical proximity effects. This relationship is polynomial-fitted to generate the contrast curve C(D), which corrects for the nonlinear photoresist response, enabling accurate dose-to-height mapping. In (b), sinusoidal dose patterns with varying spatial frequencies are used to determine the modulation transfer function (MTF). The MTF is calculated as using MTF(?)=?max′−?min′?max−?min, where ?max′ and ?min′ are the dose levels adjusted through the inverse mapping ?−1(?) of the contrast curve, and ? is the spatial frequency of the sinusoidal pattern.

The calibration of the contrast curve is straightforward and can be achieved by using large block patterns with varied uniform doses. The dimensions of each block significantly exceed the PSF associated with the writing lens, thereby mitigating the impact of OPE. Consequently, the contrast curve can be accurately determined through direct measurements obtained using a confocal microscope with interferometric measurement mode.

To accurately represent the contrast curve, a piecewise polynomial model incorporating two turning points is selected based on the characteristics observed in the measured data. The formulation of this fitting model is outlined as follows:

where ?(?) represents the contrast as a function of dose ?. The parameters ?low and ?high denote the contrast values at the low and high dose limits, respectively. ?low and ?high specify the dose thresholds below and above which the contrast remains constant. Within the range ?low<?<?high, the contrast curve is modeled as a polynomial function of ?, where ?0, ?1, ?2, etc., are the polynomial coefficients. The order of the polynomial can be modified to achieve low-error and robust fitting outcomes. The calibration results and robustness tests are included Supplement 1.

The MTF calibration process builds on the contrast curve established during the initial calibration step. As the MTF calibration process is also influenced by the nonlinear response of photoresist, this nonlinearity should be compensated to facilitate MTF estimation. Therefore, our strategy involves counteracting the nonlinear effect of photoresist through a digital inverse mapping, using a pre-calibrated contrast curve as a prerequisite for MTF calibration.

Sinusoidal patterns with varying periods are used as a series of patterns for MTF calibration, which facilitate the comprehensive calibration of MTF. The dose level distribution for one of these patterns is described by the equation

where ?max and ?min denote the maximum and minimum dose levels used, respectively. The selection of periods, ?, with their wide range and sampling density, is specifically chosen based on the resolution of the writing lens to effectively capture the rapid changes in the MTF curve. Typically, the MTF curve in most imaging systems displays a smooth profile, with a sharp decline beyond cutoff spatial frequency. This characteristic allows for a limited selection of periods in our calibration efforts. Following fabrication, we measure the features using a high-resolution confocal 3D microscopy, chosen for its ability to resolve details much finer than the minimum pattern period. The measured height profile then undergoes inverse mapping ?−1(?) against the contrast curve to compensate the nonlinear response of the photoresist. This step allows for an accurate estimation of the actual absorbed dose level. By comparing these measurements with the target dose levels, we can calculate the MTF values at specific spatial frequencies in the following manner:

where the spatial frequency ?pattern in unit of µµm−1 is the reciprocal of the period of the sinusoidal pattern. ?max′ and ?min′ represent the extremities of the dose distribution post inverse mapping, while ?max and ?min indicate the extremities of the original dose distribution. Therefore, the MTF value at a specific spatial frequency ?pattern can represent the reduction factor of the amplitude of the prescribed dose pattern caused by OPE. The MTF values are computed at specific frequencies corresponding to the employed pattern periods. To accurately model the MTF curve, we adopt a two-term exponential function:

incorporating fitting parameters ?, ?, ?, and ? to effectively capture the characteristic shape and features of the MTF curve. More details related to model calibration and fitting can be found in Supplement 1. While the MTF model assumes rotational symmetry, we observed that some systems exhibit minor directional variations in their MTF, while others demonstrate more pronounced asymmetry. For systems with minimal asymmetry, the rotationally symmetric model is sufficient. However, for systems where directional variations are significant, a rotationally symmetric model may not adequately describe the behavior of the MTF. To address this, we propose an enhanced 2D MTF model. This approach reconstructs the MTF in two dimensions by estimating MTF values at multiple angles and applying interpolation techniques such as spline fitting and mirror symmetry. This enhanced model provides improved accuracy for systems with pronounced directional variations. For practical applications, we recommend a dual-model framework: the rotationally symmetric model is suitable for systems with negligible asymmetry, while the enhanced 2D MTF model is better suited for systems where directional effects significantly impact performance. Details of the experiments validating the presence of directional variations and the methodology for constructing the 2D MTF are provided in Supplement 1.

3. EXPERIMENTAL RESULTS

A. Diffractive Beam Splitters

To validate the effectiveness of the proposed approach, we designed, fabricated, and characterized 1×3 and 1×5 diffractive beam splitters for wavelengths of 650 and 1550 nm. Diffractive beam splitters are a type of DOEs that can split a single incoming beam of light into multiple beams with specific intensity and phase distributions. To compare the performance of different design methods, we evaluated three approaches: the direct method without OPC; the conventional two-step design process, which includes phase retrieval followed by OPC; and our proposed comprehensive FID method. These approaches are referred to as “naïve,” “conventional,” and “this method” in the figures, respectively. The experimental results, including comparisons of both predicted and measured height profiles as well as diffraction efficiencies, are presented in Fig. 3.

Fig. 3. (a)–(d) Diffraction efficiency comparison for 1×3 and 1×5 beam splitters using the naive approach (direct height-to-dose mapping), conventional OPC, and the proposed method (“this method”). Each sub-figure shows (left to right) diffraction efficiency across methods and performance metrics. Beam splitters are designed for 650 and 1550 nm: 1×3 splitters at (a) 650 nm and (b) 1550 nm; 1×5 splitters at (c) 650 nm and (d) 1550 nm. (e)–(h) Measured profiles (colored, thick lines), simulated profiles (shaded areas with gray outlines), and target designs (gray thick lines for naive and conventional cases) obtained via atomic force microscopy. Beam splitters are designed for 650 and 1550 nm: 1×3 splitters at 650 nm (e) and 1550 nm (f); 1×5 splitters at 650 nm (g) and 1550 nm (h).

Figures 3(a)–3(d) show both the measured and simulation profiles for these beam splitters. The simulation profiles are represented as shaded regions, while the measured profiles—obtained using confocal 3D microscopy—are depicted as thick lines, enabling a detailed visual comparison. In traditional DOE designs, manufacturability constraints require pixelation of the pattern, where each pixel represents a constant height. While this results in discrete steps that may appear to have limited levels, the profiles actually incorporate 256 distinct heights optimized for diffraction efficiency. However, the pixelated nature can lead to significant variations within contours, as seen in Fig. 3 for both “naïve” and “conventional” cases. These variations and discrete edges can cause challenges in achieving alignment with fabrication capabilities even with OPC, contributing to deviations between designed and fabricated profiles. In contrast, our FID method overcomes this limitation by using pixelated exposure patterns with a pixel size significantly smaller than the feature size of the structure, resulting in quasi-continuous height profiles through the fabrication model. While it does not explicitly avoid challenging features such as large-aspect ratio or sharp structures, it utilizes design flexibility to promote fabrication-friendly microstructures. In essence, this approach enhances the overall manufacturability of DOE designs by creating structures that are inherently easier to fabricate.

Figures 3(e)–3(h) displays the measured diffraction efficiencies and non-uniformities across each set of beam splitters sets for the three approaches. This figure also elucidates the total diffraction efficiency and relative uniformity among operational orders as performance metrics. These metrics are decided by the distribution of light intensity across both operational and all diffraction orders, quantitatively defined as follows:

Here, “efficiency” quantifies the proportion of light energy directed into desired diffraction orders, while “non-uniformity” assesses the consistency of intensity distribution across these operational orders. Considering both metrics, the FID method demonstrates superior performance.

B. Structured Light Generation

Unlike diffractive beam splitters that typically require the design and fabrication of one-dimensional profiles, DOEs for structured light generation necessitate the creation of two-dimensional profiles. To illustrate the practicality of our method, we fabricated 15×15 and 3×3 dot projectors. Figure 4 presents a photograph of the diffraction pattern generated by the 15×15 dot projector, demonstrating the capability of the device to produce extensive and uniform patterns across a full angular range of nearly 90°. Similarly, we designed and fabricated three DOEs for each type using different approaches: the naive method (without OPC), the conventional two-step method (phase profile design followed by OPC), and finally the proposed (FID) method. The comparison clearly reveals that our (proposed) method not only achieves the best uniformity and efficiency, but also significantly enhances the practical applicability of DOE fabricated by DLWL for structured light generation. Furthermore, we also demonstrated the ability to correct Fourier-space distortion for wide-angle DOEs. For far-field wide-angle DOEs, the angular distribution of the target intensity experiences a distortion due to the nonlinear relationship between the diffraction angles and the diffraction orders [30]. Such distortion can be compensated through the pre-distortion method (Please refer to Supplement 1 for calculation details). Figure 4(e) shows the projection results designed by the FID method with pre-distortion.

Fig. 4. 15×15 dot projector with a full angular range of approximately 90°. (a) Optical setup. Captured diffraction pattern for samples designed and fabricated using (b) naïve method (direct height-to-dose mapping), (c) conventional OPC, and (d) proposed FID method. (e) Proposed FID method (with Fourier-space distortion correction).

The efficiency measurements for the 3×3 dot projectors are depicted in Fig. 5 with three different approaches. The FID technique improves diffraction efficiency by 19.1% and 10.7% over the naive method and conventional OPC, respectively, for 650 nm devices and by 13.6% and 9.7% for 1550 nm devices. Besides, the FID method drastically reduces dot array non-uniformity to 14.3% at 1550 nm and 7.6% at 650 nm, representing a remarkable improvement over both the conventional OPC and the naive method. Supplement 1includes photographs of the dot arrays captured by the camera and surface profiles of the devices, presenting visual evidence of the improved uniformity and efficiency achieved through our proposed technique.

Fig. 5. Measured diffraction efficiencies for 3×3 beam splitters operated at (a)–(c) 650 nm and (d)–(e) 1550 nm. Operating diffraction orders are marked as red dots. Methods compared include the naive approach (direct height-to-dose mapping), the conventional OPC, and the proposed FID method.

C. Far-field Digital Hologram

To validate the scalability and efficacy of our proposed method, we first targeted far-field digital holograms operating under Fraunhofer diffraction conditions. These holograms were designed with a pixel resolution of 800×800 but fabricated as 3000×3000 arrays after periodic extension, corresponding to a physical size of 5mm×5mm at a wavelength of 650 nm. Two illustrative target patterns were selected: a grayscale image of the mathematician Gauss and a binary image of the Tsinghua University logo.

Figure 6 presents the resulting diffraction patterns. The naive design method, depicted in Figs. 6(b) and 6(e), produced suboptimal image fidelity due to significant fabrication errors. By contrast, the conventional two-step approach (phase retrieval with subsequent OPC) improved fidelity as shown in Figs. 6(c) and 6(f). The proposed FID method, illustrated in Figs. 6(d) and 6(g), achieved the highest fidelity, effectively reducing fabrication errors and closely matching the target images.

The convergence of the FID method is visualized in Fig. 8(a). Using the Adam optimizer, satisfactory quality was achieved within 4000 iterations, with the entire process completed in under 3 min on a consumer-grade NVIDIA GeForce RTX 2060 GPU. This efficiency underscores the practicality of the method for far-field digital holography.

D. Large-scale Digital Hologram

To assess the scalability of the FID method for large-area, high-resolution DOE designs, simulations and experiments were conducted for designs ranging from tens to hundreds of millions of pixels. A key feature of the method is its scalability, primarily constrained by GPU memory. Memory requirements increase linearly with design size, allowing GPUs like the NVIDIA GeForce RTX 4090 (24 GB memory) to handle designs up to 6000×6000 pixels (36 million pixels). High-performance GPUs, such as the NVIDIA H100 (80 GB memory), support designs exceeding 10000×10000 pixels (100 million pixels).

For experimental validation, a large-scale hologram operating under near-field (Fresnel diffraction) conditions was designed with 6000×6000 pixels and subsequently fabricated. This choice of optical model demonstrates the adaptability of the method to varying diffraction regimes. Figure 7(a) illustrates the optical setup, while Fig. 7(b) presents the fabricated hologram. In Fig. 7(c), human images from the CelebA dataset [31] and animal images from the AFHQ dataset [32] are used as target patterns for comparison. The conventional methods produced patterns with significant deviations from the target, including reduced sharpness and increased distortion. In contrast, the FID method achieved high fidelity, with experimental diffraction patterns closely matching the target figures, confirming its accuracy and scalability.

Optimization tests for larger holograms, such as 10000×10000 pixels, demonstrated completion in tens of minutes on an NVIDIA H100 GPU. Figures 8(b) and 8(c) illustrate the usage of GPU memory and the optimization time on various design scales, further demonstrating the efficiency of the method in resource utilization.

Fig. 6. Far-field holographic pattern generation. (a) Optical setup. Diffraction patterns generated by devices using (b), (e) naive method (direct height-to-dose mapping); (c), (f) conventional OPC; and (d), (g) proposed FID method.

Fig. 7. Large-scale holographic pattern generation. (a) Optical setup. (b) A large-scale hologram with 6000×6000 pixels (36 million pixels) compares conventional methods (left), experimental results (middle), and the target figure (right).

Fig. 8. Quantitative evaluation of (a) convergence curve of the FID method, (b) GPU memory usage, and (c) optimization time across various design scales for large-scale optimization.

The robustness and adaptability of the FID method were evaluated using two distinct laser writing systems: a laser scanning system, which sequentially scans the substrate with a focused laser beam, and a projection-based system, which projects entire patterns using a digital micromirror device (DMD). For each system, different objective lenses or resolution modes were employed to further assess adaptability. Both systems successfully fabricated the same 6000×6000 pixel hologram, achieving consistent diffraction patterns (see Figures S27 and S28 in Supplement 1). A single calibration step sufficed to adapt the method to different machines, underscoring its practicality for diverse fabrication setups. Comprehensive details are available in Supplement 1, further highlighting the method’s potential for real-world applications.

E. Extended Depth-of-Focus Imaging

Advances in DOE have revolutionized imaging systems, allowing engineering of point spread functions to enhance functionalities such as superior detection [13], extended depth-of-field (EDOF) [3,33,34], and imaging through obstructions [4]. Central to these advances is the use of customized planar DOE lenses, which underscore the importance of meticulous DOE design and optimization for augmented imaging and extended focus capabilities.

In this context, our work demonstrates the application of DOEs fabricated by DLWL for EDOF imaging, highlighting the wide applicability of the proposed method. Figures 9(a) and 9(b) show infrared images obtained using EDOF DOE lenses designed and fabricated via the conventional OPC method and our proposed FID method, respectively. The scene, illuminated by an LED at approximately 850 nm, spans a range from 50 to 1200 mm. Figure 9(c) provides a hand-drawn schematic of the setup of the scene, helping readers better understand the arrangement of objects within the imaging environment. The experimental results reveal that our method significantly outperforms the conventional approach, particularly in aspects such as contrast and sharpness. Our methodology results in finer details in the image, such as shadows and specific areas on a curtain situated approximately 1200 mm from the imaging sensor, which are especially evident in the dashed rectangular regions.

Fig. 9. EDOF imaging results using DOEs fabricated by (a) conventional OPC method and (b) proposed FID method. (c) Hand-drawn schematic of the scene setup. Confocal microscope images of DOEs at 5×, 2 0×, and 50× magnifications: (d)–(f) correspond to the conventional OPC method, while (g)–(i) correspond to the proposed FID method.

The EDOF lens design ensures that the point spread function remains relatively constant across a range of focal lengths, enabling good imaging quality over an extended depth range. This is achieved by sampling uniformly within the focal length range [fmin,fmax], corresponding to the working distance range [?min,?max], as determined by Newton’s imaging equation. For a given plane, the light intensity distribution is denoted as ??, and the loss function is defined as

where ?0 is the diffraction-limited Airy pattern distribution, approximated by a two-dimensional Gaussian. This multi-objective loss function enables lens optimization to achieve PSFs closer to ideal imaging across multiple planes by iteratively modifying the exposure dose distribution.

This study adopts a rotationally symmetric phase profile as the core design for EDoF imaging, which differs from the commonly used cubic phase or its variants [35]. While cubic phases are popular due to their simplicity and separability in Cartesian coordinates, rotationally symmetric designs, such as logarithmic aspheres [36] and axicons [37], also provide effective EDoF imaging functionality. Similar examples leveraging rotational symmetry can be found in [3,38], demonstrating high-performance EDoF lens optimization under rotational symmetry constraints. To simplify the optimization process, rotational symmetry is enforced as a constraint in our method, allowing the parameterization of the two-dimensional dose distribution into a one-dimensional profile, as illustrated in Supplement 1 Figure S36. This approach significantly reduces the number of optimization variables without compromising performance, resulting in a lens profile that is inherently rotationally symmetric.

Confocal microscope images of the fabricated DOEs are presented for both methods at 5×, 20×, and 50× magnifications: Figs. 9(d)–9(f) correspond to the conventional OPC method, while Figs. 9(g)–9(i) correspond to the proposed FID method.

4. DISCUSSION AND CONCLUSION

This study introduces a fabrication-integrated design method that combines physics-informed modeling with scalable computational techniques, advancing the design and manufacturing of DOEs. Validated through experiments on hundred-million-pixel designs, computational scaling analyses, and fabrication variability evaluations, the FID approach demonstrates high precision and adaptability across different laser direct writing lithography systems, maintaining excellent optical fidelity and supporting centimeter-scale DOE designs for applications in holography, beam shaping, and structured light generation. While streamlining the design process and improving both performance and efficiency, it also eliminates the need for complex calibration procedures, such as extensive scanning electron microscope measurements or fabrication parameter estimations. Even under simplified fabrication models, the FID method exhibits remarkable adaptability and practical value for high-precision grayscale lithography, achieving significant performance gains without added procedural complexity.

Despite these advancements, challenges remain. The simplified lithography model used in this study, while computationally efficient, lacks the precision of advanced neural network-based methods, which often require significant retraining for new scenarios. Neural network-based approaches offer higher model fidelity but are constrained by their computational demands and limited interpretability. Figure S29 in Supplement 1 highlights these differences, demonstrating how the proposed FID method offers a balance between computational efficiency and practical integration advantages, particularly in scalable DOE design. Additionally, constraints in budget and software accessibility prevent direct comparisons with commercial OPC tools, limiting the scope of evaluation. Furthermore, fabrication errors inherent in grayscale direct laser writing, combined with variability in development times and material properties, can affect DOE profile accuracy despite standardized procedures designed to mitigate these issues.

Future work will focus on optimizing fabrication techniques. Enhancements to the advanced lithography model are expected to improve its precision while retaining computational efficiency. Incorporating advanced optical feedback mechanisms, such as real-time optical monitoring during the lithographic development stage, could further enhance reliability in fabrication. These advancements will enable the production of high-resolution DOEs tailored for specialized applications, establishing a robust foundation for next-generation optical element manufacturing.