Advances in modeling and optimization for two-photon lithography

Advances in Modeling and Optimization for Two-Photon Lithography Valeriia Sedovaa, Florie Ogorb, Jo¨ el Roverab, Odysseas Tsilipakosc, Jonas Wiedenmannd, Kevin Heggartyb, Andreas Erdmanna aFraunhofer IISB (IISB), Erlangen, Germany bIMT Atlantique (IMTA), Brest, France cTheoretical and Physical Chemistry Institute (TPCI), National Hellenic Research Foundation, Athens, Greece dHeidelberg Instruments Mikrotechnik GmbH (HIMT), W¨ urzburg, Germany Abstract. Background: Two-photon polymerization (TPP) is one of the most promising methods for the fabrication of metasurfaces due to its ability to create complex, high-resolution nanostructures. However, fabricating these structures using TPP is complex, and predicting the results of the fabrication process is challenging. Aim: This paper aims to address these challenges by demonstrating how different modeling techniques can be employed to support the fabrication-aware design of metasurfaces. Approach: We introduce and explore three modeling techniques: a simple Threshold model, a Compact model, and a Full model of Polymerization. Each model offers different levels of complexity and accuracy. We assess how well each model performs and what limitations they have, using practical examples to show how they can guide the fabrication process. Results: Our comparison highlights the advantages and limitations of each modeling approach. The basic Threshold model provides a general overview, focusing solely on the optical aspects and using a simple threshold for the photoresist. Thus, it lacks detailed descriptions of resist behavior. The Compact model is semi-empirical, focusing on simplified chemical dynamics of a single species while including essential photochemical processes. In contrast, the Full model of Polymerization is the most advanced, offering a detailed description of various species involved in the process, such as monomers, polymers, and quenchers. Although it provides the most accurate predictions, it is also the most complex and computationally demanding. Conclusion: By comparing these modeling approaches, we show that the decision on the most appropriate model depends on the specific requirements of the metasurface fabrication process. This analysis helps researchers and engineers determine the most suitable modeling approach for their work. Keywords: metasurfaces, two-photon lithography, computational modeling, two-photon polymerization. *Valeriia Sedova, [email protected] 1 Introduction Metalenses and planar optics are gaining significant attention due to their potential to miniaturize and optimize optical systems.1These technologies enable reduced assembly sizes, novel functionalities, and enhanced performance of optical components. This is achieved through metasurfaces—arrays of sub-wavelength structures that modulate incident light at a very fine scale. The creation of a metasurface involves engineering meta-atoms on a surface to create structures with 1 dimensions (pitch and thickness) smaller than the incident light’s wavelength.2A critical challenge in metasurface technology lies in the fabrication of these components, which is essential for both laboratory research and high-volume production. The field is experiencing rapid advancements in fabrication techniques, each with its unique advantages and limitations. These are significantly influenced by application requirements, material selection, and the available production infrastructure.3–5 In contrast to traditional lithographic and pattern transfer methods such as nanoimprint lithography, deep ultraviolet (DUV) projection lithography, and electron beam lithography, Two-Photon Polymerization (TPP) presents a novel approach with significant advantages. TPP, a form of microscale 3D printing, utilizes a focused laser to induce polymerization in a photosensitive material. This method enables the creation of structures with high aspect ratios due to the nonlinear interaction of the laser with the material.6This interaction allows for unmatched precision in the three-dimensional shaping of materials.7Such capabilities make TPP particularly valuable for applications that require intricate control over material structure. TPP is distinguished by its ability to directly ’write’ sub-wavelength structures into the material, offering significant versatility in threedimensional design. The process, which has seen substantial growth in recent years, enables the construction of complex geometries that were previously challenging or even impossible to achieve with conventional lithographic methods. However, this method is not without trade-offs. The fabrication of larger patterns by the serial writing of many small voxels is highly time-consuming. Additionally, the polymers commonly used in TPP typically have low refractive indices, which result in weak optical index contrast. This weak contrast can hinder the strong confinement of light. To address this, precise control over the aspect ratios of structures is required to overcome these material limitations and achieve the desired optical performance. Moreover, optical and chemical 2 phenomena in TPP can lead to proximity effects and deviations from the intended design. These effects make it harder to reproduce the target structures accurately. Careful adjustments to process parameters are needed to minimize these issues. Given these considerations, modeling becomes a crucial tool in supporting the optimization of complex fabrication processes. By enabling fabrication-aware design, modeling can help predict and control the outcomes of TPP, thereby improving the efficiency and reliability of producing high-quality metasurfaces. Furthermore, transferring modeling techniques from semiconductor fabrication to the fabrication of three-dimensional metasurfaces presents an opportunity to leverage established methodologies for new applications. While TPP is gaining traction in research and industry, one of the key challenges is modeling the process. Few mathematical models exist that can predict the diameter, length, and overall shape of the voxel, the fundamental building block in TPP. Several notable contributions have been made toward modeling the TPP process. DeVoe et al.’s research on voxel shapes in two-photon microfabrication highlights the complexities involved. Their work demonstrates the nonlinear relationship between laser dose and voxel shape, where higher doses result in highly asymmetric voxels and lower doses yield nearly spherical ones.8These experimental observations underscore the complexity of voxel formation, demanding a deeper understanding of the underlying polymerization kinetics to explain how voxels grow over different time scales. Although threshold behavior has been proposed by various researchers working on the TPP process, the precise mechanism behind thresholding remains unclear due to experimental limitations.9–11 Expanding on foundational studies, Somers et al. offer a comprehensive overview of the physics governing TPP, covering the fundamental interactions of light and material properties that drive three-dimensional printing at the nanoscale.12 This work provides essential insights into the parameters affecting voxel shape and stability, laying the groundwork for applied 3 modeling techniques. Serbin et al., for example, developed a model predicting changes in radical concentration during polymerization, both spatially and temporally, though this approach did not account for factors such as molecular diffusion, polymerization kinetics, or temperature variations that can significantly affect voxel formation and resolution.13 Further refining these approaches, recent studies have focused on enhancing practical control over voxel dynamics. Fourkas et al. delve into critical aspects of TPP, including the nonlinear behavior of photoresist materials and threshold effects that impact voxel geometry and quality, emphasizing the importance of laser exposure conditions and photoresist properties in effectively achieving high-resolution structures.14 Mueller et al. extend these findings by exploring specific reaction mechanisms, demonstrating how laser parameters and photoresist composition can be optimized to control polymerization onset and voxel stability under diverse conditions.15 Additionally, Pingali and Saha introduce a machine learning-based surrogate model, addressing the computational challenges in traditional finite element modeling by predicting printability in projection TPP.16 This model enables the rapid exploration of parameter spaces for optimal printing conditions, marking a significant advancement in scaling TPP processes without sacrificing print fidelity. Lastly, Xing et al. incorporated radical kinetics into their time-integrated model, offering valuable insights into polymerization dynamics, although this model simplifies TPP as a steady-state process, excluding key thermal effects that are crucial in controlling voxel size and shape during prolonged laser exposure.17 In light of these efforts, this paper compares three models of varying complexity that aim to enhance the understanding and control of the TPP process. These models include the Threshold model, the Compact model, and the Full model of Polymerization. The Threshold model focuses solely on the optical aspects of TPP, determining the minimum conditions required for initiating the polymerization process. In contrast, the Compact and Full models incorporate considerations of 4 the photoresist’s behavior, detailing the chemical reactions and physical changes that occur during and after exposure to light. Following the detailed model descriptions, we transition to model applications. This section showcases how modifications of the models’ parameters influence the geometry of the resulting voxels—the fundamental units of structure in TPP. Using experimental data provided by our partners, we demonstrate the calibration of our models against measured voxel dimensions achieved under various exposure conditions. This validation shows how different models can handle the experimental data and their potential utility in optimizing the TPP process. Subsequently, we explore the adaptation of these models for different TPP setups, illustrating their ability to adapt to various fabrication scenarios. An initial exploration into the application of these models for creating metasurfaces is then presented, marking a first step towards using TPP modeling in the development of advanced optical devices. The manuscript concludes with a summary of our findings and an outlook on future research directions. 2 Two Photon Polymerization for 3D fabrication Before going into the detailed modeling of TPP, we will describe the fundamental mechanisms behind this process, which enable the fabrication of complex 3D structures. TPP is a technique for creating microand nanostructures by selectively solidifying a liquid resin. It relies on the nonlinear absorption of two photons by a photosensitive material, typically a mixture of monomer and photoinitiator. When an ultrashort laser pulse is tightly focused inside the resin, two photons are absorbed almost simultaneously at the focal point, generating radicals. These radicals initiate a chain polymerization reaction, converting the liquid resin into a solid polymer. The polymerization occurs only within the focal volume of the laser, allowing for precise spatial control. In a negative-tone resist, regions outside the focus remain unexposed and can be easily washed away 5 Fig 1 Point Spread Functions: (A) Single-Photon Absorption (SPA) and (B) Two-Photon Absorption (TPA). after processing. The localized polymerization ensures high resolution and enables the creation of intricate 3D geometries that would be difficult to achieve using traditional lithography techniques. The number of photons absorbed per molecule per pulse, na, can be described by the following equation:18 na≈δ2P2 avg τpf2 pNA2 2ℏcλ2 ,(1) where δ2is the two-photon absorption cross-section of the photoinitiator, Pavg is the average laser power, NA is the numerical aperture of the focusing lens, τpis the pulse width, fpis the pulse repetition rate, ℏis the reduced Planck’s constant, cis the speed of light, and λis the laser wavelength. This equation highlights the quadratic dependence of photon absorption on the input laser power, in contrast to single-photon processes that exhibit a linear dependence as can be seen in Figure 1. This nonlinearity confines polymerization close to the focal point, providing the ability to fabricate structures with resolutions beyond the diffraction limit of the focusing optics. After the laser initiates polymerization, a gradual transformation of the monomers into a highmolecular-weight polymer occurs, solidifying the material in the irradiated region. Controlling 6 the spread of the polymerization reaction is critical and can be controlled by using quenchers or oxygen to limit radical diffusion. Through precise scanning of the laser focus, complex 3D structures can be constructed voxel by voxel (volumetric pixel). The size and shape of these voxels depend on multiple factors, including laser power, exposure time and the material’s properties. This voxel-based approach is the fundamental building block of the TPP process and understanding its evolution is key to achieving high resolution and control over the final structure. 3 Model description Modeling TPP process differs fundamentally from single-photon polymerization due to its spatial confinement, nonlinear polymerization dynamics, and the pulsed nature of the laser source.10 These factors are crucial for accurately predicting the behavior of photopolymers. Here, we present three distinct models that range from simple to complex, capturing various aspects of TPP with different degrees of detail. Typically, these models consist of two parts: the optical part that simulates the image formation and the light interaction with the resist, and the resist part that delves into the material’s response. This approach allows for a thorough understanding of TPP process, from light-matter interaction to the final structure, guiding the fabrication of microstructures in TPP. 3.1 Optical model In the optical modeling, our goal is to compute what is known as the bulk image or 3D representation of light intensity within the photoresist, which is obtained for a specific optical setup. Here we will focus on the image of a point object - the so called - point spread function (PSF). The PSF is computed using Dr. LiTHO,19 a tool specifically designed for lithography simulation. The Dr.Image module within Dr. LiTHO has methods to calculate aerial and bulk images in partial co7 Fig 2 This figure presents two representations of a Point Spread Function (PSF) visualization. (A) shows a ”Normal PSF” plot displaying the light distribution on a standard linear scale. (B) shows a ”Log PSF” plot, presenting the same distribution but with a logarithmic scale for intensity values, enhancing the visibility of details across a wide range of intensities. herent image projection systems, which are used in lithography scanners for semiconductor lithography. It employs the Abbe-method tailored for image simulation.20 The shape of the PSF depends on the wavelength (λ) of light, numerical aperture (NA) of the system, and the refractive index of the immersion fluid (nimmers). Additionally, the shape of the Point Spread Function (PSF) is significantly determined by the defocus level (defocus) and system reduction (reduction), alongside the properties of the wafer stack, particularly the photoresist’s thickness and refractive index (nres). Together, these parameters define the light distribution within the photoresist, referred to as the bulk image, shown in Figure 2. This distribution is crucial for understanding how the photoresist behaves during exposure and plays a key role in accurately modeling the lithography process. 3.2 Resist model Following the modeling of light focusing and interaction with the wafer stack outlined in the Optical Model Section 3.1, we now turn our attention to the Resist Model. This section delves into the mathematical descriptions of the photochemical reactions and the subsequent alterations in the 8 photoresist due to light exposure. The Resist Model is critical as it captures the changes in the photoresist composition triggered by the absorbed intensity, a key process that defines the success of TPP. This section introduces three resist models, ranging from basic Threshold model to complex descriptions of kinetics and diffusion processes, which are essential for accurately simulating TPP. 3.2.1 Threshold model The Threshold model simplifies the process of TPP by distinguishing polymerized from nonpolymerized regions using a straightforward method. It applies an intensity threshold to the optical model’s point spread function (PSF) to determine the volume where polymerization will occur. This model does not describe details of the resist’s behavior, focusing exclusively on the optical interactions. In this approach, polymerization begins where the squared intensity within the PSF exceeds a certain threshold, reflecting the quadratic response characteristic of TPP processes. The threshold is approximately inversely proportional to the exposure dose.21 At low exposure doses, only the center of the PSF—with squared intensity exceeding higher threshold values—is polymerized. As the exposure dose increases, polymerization extends over larger volumes corresponding to regions with lower squared intensity. This simplification permits rapid calculations of voxel formations and can be particularly effective for initial approximations or in scenarios where detailed material dynamics are less critical. Figure 2presents the output of the Threshold model, showing the PSF with a defined threshold contour and the resulting 3D shape that represents the boundary of polymerization. This visual representation demonstrates the immediate outcome of applying the threshold, marking the limits of the polymerized area within the photoresist. 9 Dark Phase. The Dark Phase begins after exposure reactions, occurring without light, and continues until resist development. Besides polymerization, termination and quenching reactions also occur during this phase. The polymerization process is driven by a nonlocal photopolymerization diffusion model:29 ∂[R] ∂t =DR∇2[R]−(2kt[R]2+kq[Q][R] + kp[M][R]),(10) ∂[Q] ∂t =DQ∇2[Q]−(kq[Q][R] + kq[M∗][Q]),(11) ∂[M] ∂t =DM∇2[M]−(ki[M][R] + Ikp[M]∗(t, x′)[M](t, x′)G(x, x′)dx′),(12) ∂[M∗] ∂t =D∗ M∇2[M∗] + ki[M][R]−(2kt[M∗]2+ktp[M∗] + kq[M∗][Q]),(13) where the term Rrepresents the amount of radicals, Mthe amount of monomers, [M∗]- the amount of growing active monomer chains, Qthe amount of quencher, DMis monomer diffusion coefficient, DQquencher diffusion coefficient, ktrate constant for termination, kprate constant for propagation, ktp - rate of termination by primary radicals, kqrate constant for quenching, kirate constant for initiation. A convolution is utilized in Equation 12 to model the effect of initiation of propagation reactions at a spatially distant position. The function G(x, x′)is a Gaussian function applied in the convolution product in Equation 14: G(x, x′) = 1 √2πσ exp −(x−x′)2 2σ2.(14) The term σin Equation 14 is a measure of the spatial extent of propagation reaction that results from a single radical. The extent of this spread varies with the type of resist material used. At 16 Fig 5 Spatial distribution of species during the dark phase of polymerization. (A) Radical distribution, R, indicating the concentration of radicals. (B) Monomer distribution, M, representing the remaining monomers. (C) Polymer distribution, M∗, showing the formation of polymerized regions. (D) Quencher distribution, Q, showing the spatial variation of the quenching agent. the beginning of the dark phase, there is a relatively high concentration of radicals in the resist. However, as the number of radicals decreases, the formation of new polymer chains starts to slow down, while the reactions that stop and stabilize the chains start to take over. This is why quenching is a key part of the reaction — it helps keep the size of the final features under control. Figure 5 presents the concentration of: radicals, monomer, polymer and quencher after the dark phase. Development. For the development step, we utilize the Mack rate model and fast marching algorithm as described in subsection 3.2.2. In summary, the Full model of Polymerization introduces a comprehensive approach to simulate TPP by closely tracking the progression of multiple reactive components within the resist. 17 Unlike the Compact model that simulates diffusion of just a single species, this model captures the interplay between different species like monomers, polymers, radicals, and quencher. The quenching of the photopolymerization reaction, typically by oxygen, provides an important component of the model with significant impact on the resulting voxel shape. Although this model offers a more detailed simulation, capturing the complex interactions during TPP, it results larger computation times and more, mostly unknown, model parameters. 4 Model application 4.1 Impact of parameters on voxel shape The shape of the voxel is a key metric of the TPP process, as it directly influences the architecture of the resulting metasurface. In this section, we explore the impact of critical parameters on the voxel shape. Starting with imaging, the NA of the imaging system plays an important role in determining the resolution and quality of the voxel shape, which subsequently impacts the resulting metasurface structure. As the NA increases, it allows for a tighter focus of the laser beam, which can lead to a more defined and reduced voxel size. This effect is evident in Figure 6, where different NA values show a clear correlation with voxel dimensions. The top row of images illustrates how the NA affects the point spread function (PSF), which in turn influences the intensity distribution within the resist. A higher NA typically results in a more concentrated light distribution, leading to smaller voxels as shown in the DART images in the second row. The bottom left graph underscores this relationship by plotting voxel dimensions against NA values, demonstrating that both the diameter and length of voxels decrease as the NA increases. The bottom right plot shows the voxel shape 18 Fig 6 Impact of NA on voxel shape. The top row depicts the intensity distribution of light for different NA values, illustrating how the light focus narrows as NA increases. The middle row shows the corresponding DART images of voxels, revealing the progressive sharpening and reduction in size of the voxels higher NA values. The bottom left graph illustrates the dependency of voxel dimensions on NA. The bottom right graph presents the voxel shape’s dependency on NA values. dependency on NA values, offering a visual representation of how varying the NA can alter the voxel profile, with lower values creating more elongated shapes. The refractive index of the photoresist is another important parameter. The refractive index determines how light propagates through the resist, which in turn affects the voxel’s formation, influencing the aspect ratio of the voxel—crucial for defining the metasurface’s functionalities. Figure 7demonstrates that the refractive index has a significant impact on the voxel length and the aspect ratios of the voxel. A higher refractive index leads to an extended focal depth, which also tends to stretch the voxel along the z-axis, effectively increasing its length. The DART images visually confirm these changes, where an increasing refractive index correlates with a more 19 Fig 7 Impact of Photoresist Refractive Index on voxel characteristics. pronounced voxel elongation. The contours in the lower right of Figure 7suggest in increasing (top-down) asymmetry of the voxel for larger refractive index of the photoresist. This asymmetry is caused by spherical wavefront aberrations introduced by the interface between the immersion fluid and the photoresist. In practice such asymmetries are avoided by an index match between photoresist and immersion fluid. In the exposure step of TPP, the power of the laser source is a critical factor that influences the polymerization sensitivity and thus the voxel dimensions. As depicted in Figure 8, we can see how varying the laser power affects the voxel’s shape. With low power values, the voxel appears smaller due to the limited energy available to initiate the polymerization process. As the power increases, we observe a substantial increase in both the diameter and the length of the voxel, evidenced by the widening and elongation in the DART images. In the dark phase of TPP, the parameter σplays a significant role in the diffusion process, which 20 Fig 8 Impact of power on voxel characteristics. critically impacts the development of the voxel. The term σrepresents the root mean square of the diffusion length and dictates how far the polymerizing species can travel within the photoresist during the time between exposure and development (Equation 14). Figure 9illustrates how different values of σaffect the voxel. With smaller σvalues, indicating limited diffusion, the resulting voxels are more compact, as visualized in the DART images. As σincreases, we observe a broadening effect; the polymerizing species diffuse more widely, leading to a growth in the voxel’s size. The graphs underscore the dependency of voxel dimensions on σ. The left graph plots an increasing trend in both diameter and length of voxels with rising σvalues, while the right graph provides a contour representation of the voxel shape change in response to different σvalues. Here, larger σ values result in a more pronounced spread in the voxel shape. In conclusion, the voxel formation in TPP is significantly governed by the interplay of mul21 Fig 9 Impact of σon voxel characteristics. tiple parameters across different phases of the process. From the initial imaging, where NA sets the stage for resolution, to the exposure phase, where laser power and exposure time dictate the extent of polymerization, and finally to the dark phase and development phase, where diffusion and developer selectivity refine the final structure—each parameter has a profound impact on the outcome. By carefully examining and optimizing these parameters, we can achieve better control over the TPP process. This allows for the creation of metasurfaces with complex geometries and functionalities tailored to specific applications. The insights gained from the dependency of voxel dimensions on factors such as power,NA,σ, and N, along with other contributing parameters, provide us with the ability to tailor the TPP process and produce optimized metasurfaces. 22 4.2 Model calibration and verification 4.2.1 Model parameter calibration using experimental voxel data The predictivity and flexibility of the models were evaluated through calibration against experimental data, using voxel measurements as benchmarks. The experimental setup employs a femtosecond (fs) laser system for high-resolution TPP. The collected data illustrates the relationship between laser power, exposure time, and voxel dimensions, providing a foundation for model calibration. In the following, we calibrate our three models with the described experimental data for voxels. The aim is to not only validate our models but also to get insights that could refine our understanding and enhance the fidelity of the TPP simulation. Threshold model. The Threshold model predicts the regions of polymerization based on an intensity threshold (Section 3). As shown in Fig. 10 A, the model is effective at estimating the shape of smaller voxels under low power but fails to accurately predict larger voxel dimensions at higher powers. This limitation is due to the model’s inability to consider diffusion and kinetics, which become important at higher powers. Similar observations are reported in the literature.8,30 In summary, while the Threshold model serves as a useful tool for rapid, preliminary simulations, particularly at lower power settings, it is insufficient for practically relevant high exposure powers. Compact model. The Compact model improves upon the Threshold model by accounting for diffusion effects during polymerization. Figure 10 B shows that this model provides better predictions of voxel diameters across varying power levels. However, the model struggles to predict voxel lengths accurately, suggesting that additional factors, such as anisotropic diffusion or other mate23 Fig 10 The graph illustrates the relationship between laser power and voxel dimensions for three models: (A) Threshold Model, (B) Compact Model, and (C) Full Model of Polymerization. Experimental and simulated data are shown for both voxel length (blue) and diameter (red) across all models. rial properties, may need to be included. Despite this, the Compact model represents a significant improvement over the Threshold model in matching experimental data for voxel diameters. Future improvements of the model could incorporate anisotropic diffusion31 to better simulate the directional dynamics of the polymerization process. Alternatively, a more comprehensive model that includes a wider range of physical and chemical interactions during TPP could help to overcome to the current limitations. Such developments would significantly enhance the model’s predictive capabilities, e.g., a quenching effect, providing a more accurate tool for optimizing the TPP process. Full model of Polymerization. The Full model of Polymerization integrates temperature effects, multi-species diffusion, and dark-phase reactions, providing a comprehensive view of the photopolymerization process (Section 3). Calibration results (Figure 10 C) demonstrate improved alignment with experimental data for both voxel diameter and length under various power conditions. However, some discrepancies remain in predicting voxel lengths, especially at higher powers. The Full model of Polymerization includes a critical radical control parameter that governs the 24 Fig 11 Calibration graphs for the Full model of Polymerization showing the relationship between laser power and voxel diameter at four distinct exposure times (1000 µs, 100 µs, 10 µs, 1 µs) interplay between diffusion kinetics and exposure time. As shown in Figure 11, this parameter can be adjusted to describe data, which were obtained by different combinations of exposure power and exposure time. This showcases the model’s flexibility in adapting to the kinetic and diffusion dynamics for different exposure dynamics. This model stands as the most advanced among the three we have discussed, particularly in its representation of voxel length and exposure dynamics. While there is room for improvement, especially in the more accurate prediction of length, the current calibration results represent the most accurate trend of voxel dimensions in response to power variations. Moving forward, the application of anisotropic diffusion (for example for radicals) within this model may refine its predictive capabilities further. Nonetheless, the Full model of Polymerization has already shown its potential to describe the TPP process, offering a substantial improvement in the simulation of photopolymerization dynamics. Another important (optical) effect, which was not discussed so far is the illumination of the 25 Fig 17 Illustration of double exposure photolithography process. Top row: Target ”T” shape design and the respective mask transmission patterns for the first and second exposures. Bottom row: Cross-sectional simulation results showing the summation of exposures and the resulting 3D structure of the ”T” shape as achieved through the double exposure method. distinct modeling approaches—each with varying complexity and capability. The Threshold model, our initial step, is efficient for simulations at lower powers. However, its simplicity, focusing only on the optical model without considering resist kinetics or diffusion, limits its applicability under higher power conditions where these factors become critical. It struggles to accurately predict voxel diameters when the laser power exceeds 15mW for a 10-microsecond exposure time, indicating its insufficiency for higher power applications. Experimental data and literature confirm that achieving a universal quantitative predictability of TPP processes requires more advanced solutions, such as the Compact model or the Full model of Polymerization. The Compact model emerged as an advancement, offering improved predictions of voxel dimensions by integrating diffusion effects during polymerization. This model has shown considerable promise, aligning more closely with experimental results, particularly for voxel diameters. However, it falls short in reliably predicting voxel lengths. 32 The Full model of Polymerization, the most comprehensive of the three, includes additional physical and chemical effects to the simulations by incorporating a quencher—a feature not present in the other models. This addition enables the model to better replicate the termination of polymerization, a critical phase in two-photon lithography. Furthermore, this model’s capability to adapt to various exposure times, thanks to a unique radical control parameter, sets it apart from its predecessors. A key advantage of the Full model is its ability to capture the diffusion of multiple species, including radicals, quenchers, and monomers, within the photoresist, providing a more detailed understanding of the photopolymerization process. A significant advancement was the inclusion of the filling factor in our calculations. Initially, models showed discrepancies in fitting the experimental data, particularly due to the assumption of a flat beam profile. However, upon adjusting the filling factor to account for a Gaussian illumination beam profile (e.g., a filling factor of 0.82 vs. 1.0), the models exhibited far better agreement with experimental measurements across different power levels and exposure times. This adjustment corrected the underfilling of the pupil, which reduced the effective NA and caused inaccuracies in earlier predictions. These simulations also demonstrate that the pupil filling can be used to control the shape of the voxel, enabling its adaptation for specific applications where precise voxel geometry is required. We conducted a thorough analysis to discern the influence of various parameters on the voxel’s final shape. By systematically examining each step of the process, from exposure reactions to the development phase, we established how particular adjustments can modify the shape and size of the voxel. The examples of a simulated supercell of a metasurface in Figure 16 and of the formation of a T-shaped pattern by a double exposure demonstrate first steps towards process-aware design and fabrication of optical metasurfaces. 33 In future work we will enhance the computational speed and accuracy of our simulations, exploring anisotropic diffusion, and incorporating additional physical and chemical interactions that occur during TPP. The goal is to transform these models into more predictive, versatile, and universally applicable tools for metasurface fabrication. Moreover, we anticipate the expansion of our models’ applications, exploring their adaptability to novel TPP techniques and materials. As metasurfaces continue to advance, it is imperative to understand and predict the used fabrication process and their outcome. We are also exploring the integration of neural networks to provide faster and more robust computations. Acknowledgments This project has received funding from the European Union’s Horizon Europe research and innovation program under grant agreement no 101091644. Funded by the European Union. Views and opinions expressed are however those of the author(s) only and do not necessarily reflect those of the European Union or Horizon Europe. Neither the European Union nor the granting authority can be held responsible for them. References 1 D. C. Zografopoulos and O. Tsilipakos, “Recent advances in strongly resonant and gradient all-dielectric metasurfaces,” Mater. Adv. 4, 11–34 (2023). 2 Y. Zhou, I. I. Kravchenko, H. Wang, et al., “Multifunctional metaoptics based on bilayer metasurfaces,” Light: Science & Applications 8, 80 (2019). 3 B. Robben, C. Beckerleg, and L. Penninck, “Local full-wave methods for accurate modelling of large area meta-surfaces (invited paper),” in Proceedings of the SPIE,13023(13023-13) (2024). 34 4 J. Zeng, Y. Dong, J. Zhang, et al., “Broadband and high-efficiency multi-tasking silicon-based geometric-phase metasurfaces: A review,” Photonics 9(9), 606 (2022). 5 S. Zahra, L. Ma, W. Wang, et al., “Electromagnetic metasurfaces and reconfigurable metasurfaces: A review,” Frontiers in Physics 8, 593411 (2020). 6 V.-C. Su, C. H. Chu, G. Sun, et al., “Advances in optical metasurfaces: fabrication and applications [invited],” Optics Express 26(10), 13148–13182 (2018). 7 M. Malinauskas, H. Gilbergs, A. ˇ Zukauskas, et al., “A femtosecond laser-induced two-photon photopolymerization technique for structuring microlenses,” Journal of Optics 12(3), 035204 (2010). 8 R. J. DeVoe, H. W. Kalweit, C. A. Leatherdale, et al., “Voxel shapes in two-photon microfabrication,” in Multiphoton Absorption and Nonlinear Transmission Processes: Materials, Theory, and Applications,4797, International Symposium on Optical Science and Technology, SPIE, (Seattle, WA, United States) (2003). 9 S. Kuebler and M. Rumi, “Nonlinear optics-applications: Three dimensional microfabrication,” in Nonlinear Optics, 189–206, Elsevier, Oxford (2004). 10 C. N. LaFratta, J. T. Fourkas, T. Baldacchini, et al., “Multiphoton fabrication,” Angewandte Chemie 46, 6238–6258 (2007). 11 K. Lee, D. Yang, S. Park, et al., “Recent developments in the use of two-photon polymerization in precise 2d and 3d microfabrication,” Polymers for Advanced Technologies 17, 72–82 (2006). 12 P. Somers, A. M¨ unchinger, S. Maruo, et al., “The physics of 3d printing with light,” Nature Reviews Physics 6(2023). 35 13 J. Serbin, A. Egbert, A. Ostendorf, et al., “Femtosecond laser-induced two-photon polymerization of inorganic-organic hybrid materials for applications in photonics,” Optics Letters 28(5), 301–303 (2003). 14 J. T. Fourkas, “Chapter 1.3 - fundamentals of two-photon fabrication,” in Three-Dimensional Microfabrication Using Two-photon Polymerization, T. Baldacchini, Ed., Micro and Nano Technologies, 45–61, William Andrew Publishing, Oxford (2016). 15 J. Mueller, J. Fischer, and M. Wegener, Reaction Mechanisms and In Situ Process Diagnostics, 82–101 (2016). 16 R. Pingali and S. K. Saha, “Printability Prediction in Projection Two-Photon Lithography Via Machine Learning Based Surrogate Modeling of Photopolymerization,” Journal of Micro and Nano-Manufacturing 10, 031005 (2023). 17 J.-F. Xing, X.-G. Dong, W.-Q. Chen, et al., “Improving spatial resolution of two-photon microfabrication by using photoinitiator with high initiating efficiency,” Applied Physics Letters 90(13), 131106 (2007). 18 A. Diaspro, Confocal and Two-Photon Microscopy: Foundations, Applications, and Advances, Wiley-Liss, New York (2002). 19 T. F¨ uhner, T. Schnattinger, G. Ardelean, et al., “Dr.LiTHO: a development and research lithography simulator,” in Optical Microlithography XX, D. G. Flagello, Ed., 6520, 65203F, International Society for Optics and Photonics, SPIE (2007). 20 A. Erdmann, Optical and EUV Lithography: A Modeling Perspective, SPIE Press, Bellingham, Washington, USA (2021). 36 21 D. Fuard, M. Besacier, and P. Schiavone, “Assessment of different simplified resist models,” in Proc SPIE,4691 (2002). 22 F. H. Dill, W. P. Hornberger, P. S. Hauge, et al., “Characterization of positive photoresist,” IEEE Transactions on Electron Devices 22, 445 (1975). 23 C. A. Mack, “New kinetic model for resist dissolution,” Journal of The Electrochemical Society 139(4), L35–L37 (1992). 24 W. Gao, U. Klostermann, T. Muelders, et al., “Application of an inverse mack model for negative tone development simulation,” Proceedings of SPIE - The International Society for Optical Engineering 7973 (2011). 25 J. A. Sethian, “Fast marching level set methods for three-dimensional photolithography development,” in Proc. SPIE,2726, 262 (1996). 26 T. P. Onanuga, Process modeling of two-photon and grayscale laser direct-write lithography. Dissertation, Friedrich-Alexander-Universit¨ at Erlangen-N¨ urnberg, Erlangen, Germany (2019). Doktor-Ingenieur. 27 J. B. Mueller, J. Fischer, Y. J. Mange, et al., “In-situ local temperature measurement during three-dimensional direct laser writing,” Applied Physics Letters 103(12), 123107 (2013). 28 A. K. O’Brien and C. N. Bowman, “Modeling the effect of oxygen on photopolymerization kinetics,” Macromolecular Theory and Simulations 15(2), 176–182 (2006). 29 M. R. Gleeson, J. Guo, and J. T. Sheridan, “Recent developments in the npdd model,” in Optical Modeling and Design II, F. Wyrowski, J. T. Sheridan, J. Tervo, et al., Eds., SPIE Proceedings 8429, 84291C, SPIE (2012). 37 30 S. Wu, J. Serbin, and M. Gu, “Two-photon polymerisation for three-dimensional microfabrication,” Journal of Photochemistry and Photobiology A: Chemistry 181(1), 1–11 (2006). 31 F. Ogor, T. Le Deun, V. Sedova, et al., “Modelling and compensating proximity effects in massively parallelized multi-photon photoplotting,” in Proceedings of the SPIE 2024, SPIE (2024). 32 F. Ogor, T. Le Deun, V. Sedova, et al., “Modelling and compensating proximity effects in a massively parallelized multi-photon photoplotter,” 14 (2024). 33 M. Thiel, J. Fischer, G. von Freymann, et al., “Direct laser writing of three-dimensional submicron structures using a continuous-wave laser at 532 nm,” Applied Physics Letters 97(22), 221102 (2010). 34 D. K. Limberg, J.-H. Kang, and R. C. Hayward, “Triplet–triplet annihilation photopolymerization for high-resolution 3d printing,” Journal of the American Chemical Society 143(9), 3387–3395 (2021). 35 G. Perrakis, M. Kafesaki, and O. Tsilipakos, “Optical metasurfaces with two-photon lithography: design considerations for beam steering applications,” in Proceedings of SPIE, (1302314), SPIE (2024). 36 J. Rovera, V. Sedova, F. Ogor, et al., “Digital modelling of a massively parallelised multiphoton polymerisation plot process,” in Proceedings of the SPIE 2024, SPIE (2024). Valeriia Sedova is currently pursuing a PhD at FAU Erlangen-Nuremberg. She works in the lithography group, focusing on computational lithography. She completed her MS degree in Advanced Optical Technologies from the Friedrich-Alexander-Universit¨ at Erlangen-N¨ urnberg in Germany. She is also a member of SPIE. 38 List of Figures 1 Point spread function 2 Point spread function 3 Compact model simulation flow 4 Full model of Polymerization simulation flow 5 Spatial distribution of species 6 Impact of NA 7 Impact of refractive index 8 Impact of exposure 9 Impact of √σ 10 Calibration results of the Threshold model 11 Calibration graphs 12 Comparison of voxel profiles at different filling factors. 13 Voxel size with different filling factors 14 Schematic of the Two-Photon Lithography Setup 15 Overlay of simulated and experimental photopolymerization results 16 Supercell of an optical metasurface 17 Double exposure 39