# Advances in modeling and optimization for two-photon lithography

Valeriia Sedova, Florie Ogor, Joël Rovera, Odysseas Tsilipakos, Jonas

Wiedenmann, Kevin Heggarty, Andreas Erdmann

► To cite this version:

Valeriia Sedova, Florie Ogor, Joël Rovera, Odysseas Tsilipakos, Jonas Wiedenmann, et al.. Advances in modeling and optimization for two-photon lithography. Journal of Micro/Nanopatterning, Materials, and Metrology, 2025, 24 (02), 10.1117/1.JMM.24.2.023001. hal-05194777

HAL Id: hal-05194777

https://imt-atlantique.hal.science/hal-05194777v1

Submitted on 31 Jul 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Advances in modeling and optimization
for two-photon lithography

Valeriia Sedova,ª,∗ Florie Ogor,ᵇ Joël Rovera,ᵇ Odysseas Tsilipakosᶜ,ᶜ
Jonas Wiedenmann,ᵈ Kevin Heggarty,ᵇ and Andreas Erdmannᵉª

**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:** We aim 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. By 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.

**Conclusions:** 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.

© The Authors. Published by SPIE under a Creative Commons Attribution 4.0 International License. Distribution or reproduction of this work in whole or in part requires full attribution of the original publication, including its DOI. [DOI: 10.1117/1.JMM.24.2.023001]

**Keywords:** two-photon lithography; computational modeling; two-photon polymerization; metasurfaces

Paper 24096G received Dec. 29, 2024; revised Apr. 28, 2025; accepted Apr. 29, 2025; published May 24, 2025.

<footnote>ªFraunhofer IISB, Erlangen, Germany</footnote>

<footnote>ᵇIMT Atlantique, Brest, France</footnote>

<footnote>ᶜTheoretical and Physical Chemistry Institute, National Hellenic Research Foundation (NHRF), Athens, Greece, Institute of Electronic Structure and Laser, Foundation for Research and Technology Hellas, Heraklion, Greece</footnote>

<footnote>ᵈHeidelberg Instruments Mikrotechnik GmbH, Würzburg, Germany</footnote>

<footnote>*Address all correspondence to Valeriia Sedova, valeriia.sedova@iisb.fraunhofer.de</footnote>
# 1 Introduction

Metalenses and planar optics are gaining significant attention due to their potential to miniaturize and optimize optical systems.¹ These 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 dimensions (pitch and thickness) smaller than the incident light’s wavelength.² A 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.³⁻⁵

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 an alternative approach with significant advantages. TPP, a form of microscale 3D printing, utilizes a focused laser to induce polymerization in a photosensitive material. Recent progress in materials and processes for additive manufacturing has significantly increased interest in TPP, promoting an improved understanding of its capabilities and applications.⁶,⁷ Moreover, accurate computational modeling is key to using these advancements effectively, helping to better understand TPP processes and guide fabrication-aware design. This method enables the creation of structures with high aspect ratios due to the nonlinear interaction of the laser with the material.⁸ This interaction allows for unmatched precision in the three-dimensional shaping of materials, as demonstrated by Malinauskas et al., who utilized TPP to achieve highly precise fabrication of microlenses with minimal feature sizes and high reproducibility, showcasing its versatility for optical applications.⁹ Such 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 three-dimensional 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. Fabrication of larger patterns via serial voxel-by-voxel writing can indeed be highly time-consuming, a limitation common to most 3D additive manufacturing methods compared with traditional lithography techniques that pattern large areas simultaneously. This inherent slowness significantly limits throughput and scalability for large-scale production. In addition, polymers commonly employed in TPP typically exhibit relatively low refractive indices, generally in the range of ~1.5 to 1.6, such as commercially available IP resins (e.g., IP-Dip and IP-S from nanoscribe). Although higher-index materials reaching indices up to about 2.6, such as chalcogenide glasses (e.g., AsS and AsSe), have been demonstrated experimentally, most standard TPP materials still yield weaker optical index contrast compared with these specialized materials. This limited contrast poses challenges for achieving strong confinement of light, requiring precise control over the aspect ratios of structures to overcome material limitations and achieve the desired optical performance. Moreover, optical and chemical 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. Although 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. Optical models, which describe the correct superposition of polarized field components in the image plane and the impact of the polarization of the used light on the PSF are well
established in the literature.¹³,¹⁴ DeVoe et al.'s research on voxel shapes in two-photon micro-fabrication 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.¹⁵ These 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.¹⁶⁻¹⁸ Kiefer et al. expanded on the threshold concept by introducing a threshold model with an accumulation approach, describing polymerization as the result of the cumulative integration of exposure dose over time and space. The model's sensitivity figure-of-merit (FOM) provides a practical framework for comparing photoresist performance under varying experimental conditions, offering insights into the interplay between scanning speed, laser power, and photoinitiator properties for optimizing TPP processes.¹⁹ By contrast, Johnson et al. developed a more comprehensive compact model, incorporating reaction-diffusion dynamics of photoinitiators and oxygen inhibitors, as well as the self-deactivation phenomenon. Their model simulates the temporal and spatial evolution of chemical species, demonstrating how diffusion and self-deactivation dominate voxel formation and polymerization thresholds in two-photon nanolithography.²⁰ Somers et al. provide 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.²¹ This work provides essential insights into the parameters affecting voxel shape and stability, laying the groundwork for applied 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.²² Further refining these approaches, recent studies have focused on enhancing practical control over voxel dynamics. Fourkas et al. delved 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.²³ 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.²⁴ In addition, 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.²⁵ 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. Finally, 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.²⁶

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. By contrast, the compact and full models incorporate considerations of 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 toward using TPP modeling in the development of advanced optical devices. The article 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 micro- and nanostructures by selectively solidifying photoresist materials, which can be liquid or semi-solid prior to exposure. It relies on the nonlinear absorption of two photons by a photosensitive material, typically a mixture of monomer and photoinitiator. Although most TPP resins require a photoinitiator, recent studies have shown that multi-photon polymerization can also occur in pure monomers, as demonstrated by Ladika et al.²⁸ 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 exposed regions into a stable polymer network. Although under idealized conditions polymerization is confined to the focal volume, in practice, subsequent exposures can lead to polymerization outside the nominal focal volume due to proximity and memory effects of the material. In a negative-tone resist, regions outside the focus remain unexposed and can be easily washed away 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, $n_a$, can be described by the following equation:²⁹

$$n_a \approx \frac{\delta_2 P_{\text{avg}}^2}{\tau_p f_p^2} \left( \frac{NA^2}{2\hbar c\lambda} \right)^2, \quad (1)$$

where $\delta_2$ is the two-photon absorption cross-section of the photoinitiator, $P_{\text{avg}}$ is the average laser power, NA is the numerical aperture of the focusing lens, $\tau_p$ is the pulse width, $f_p$ is the pulse repetition rate, $\hbar$ is the reduced Planck's constant, c is the speed of light in vacuum, and $\lambda$ is the laser wavelength in vacuum. 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 Fig. 1. This nonlinearity confines polymerization to a highly localized volume, enabling feature sizes smaller than the diffraction-limited focal spot achievable with single-photon absorption. After the laser initiates polymerization, a gradual transformation of the monomers into a high molecular weight polymer occurs, solidifying the material in the irradiated region. Controlling the spread of the polymerization reaction is critical and can be controlled 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, material properties, and the polarization of the incident light. 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.

<figure><img src="image_2.png" /><figcaption>Fig. 1 Point spread functions: (a) single-photon absorption (SPA) and (b) two-photon absorption (TPA).</figcaption></figure>
# 3 Model Description

Modeling the TPP process differs fundamentally from single-photon polymerization due to its spatial confinement, nonlinear polymerization dynamics, and the pulsed nature of the laser source.¹⁷ 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 the TPP process, from light-matter interaction to the final structure, guiding the fabrication of microstructures in TPP.

## 3.1 Optical Model

In 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,³⁰ a tool specifically designed for lithography simulation. The Dr.Image module within Dr. LiTHO has methods to calculate aerial and bulk images in partial coherent image projection systems, which are used in lithography scanners for semiconductor lithography. It employs the Abbe method tailored for image simulation.³¹ Dr. LiTHO employs a vector imaging algorithm, which takes all polarization effects, including the correct superposition of electromagnetic field components and a Jones pupil formulation, which incorporates polarization-dependent aberration and apodization. For the simulations in this work, an ideal pupil function corresponding to an aberration-free lens is assumed. The impact of the immersion fluid, coverslip, and photoresist on image formation is included via thin-film transfer matrices. Moreover, we employed unpolarized illumination, which inherently reduces polarization-induced voxel asymmetries. As a result, the generated voxels exhibit nearly circular cross-sections, consistent with experimental observations reported by Sun et al.³² The shape of the PSF depends on the wavelength ($\lambda$) of light, numerical aperture (NA) of the system, and the refractive index of the immersion fluid ($n_{\text{immers}}$). In addition, the shape of the 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 ($n_{\text{res}}$). In our model, these refractive indices ($n_{\text{immers}}$ and $n_{\text{res}}$) are chosen to match the experimental conditions. The illumination geometry is assumed to be a Dip-In configuration. Together, these parameters define the light distribution within the photoresist, referred to as the bulk image, shown in Fig. 2. This distribution is crucial for understanding how the photoresist behaves during exposure and plays a key role in accurately modeling the lithography process. Although the intensity distribution ($I$) within the resist fundamentally controls the two-photon absorption process,

<figure><img src="image_3.png" /><figcaption>**Fig. 2** This figure presents two representations of a point spread function (PSF) visualization. Panel (a) shows a “normal PSF” plot displaying the light distribution on a standard linear scale. Panel (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.</figcaption></figure>
experimental setups typically provide parameters in terms of the average laser power (P) and exposure time (t). To bridge the experimental and theoretical aspects, we convert the measured power and exposure duration into an effective local intensity distribution using the characteristics of the optical system, specifically the PSF. This approach allows the determination of the accumulated energy dose (D) per voxel (exposed volume), calculated as

$$D = It = \frac{Pt}{A}, \qquad (2)$$

where A represents the cross-sectional area derived from the PSF. By establishing this connection between measurable experimental parameters and theoretical intensity, our modeling ensures accurate simulation of the polymerization dynamics within the resist.

## 3.2 Resist Model

Following the modeling of light focusing and interaction with the wafer stack outlined in the optical model Sec. 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 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 the 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 non-polymerized regions using a straightforward method. It applies an intensity threshold to the optical model's 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.³³ 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 2 presents 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.

### 3.2.2 Compact model

The compact model provides a semi-empirical approach to include physical/chemical details of TPP, going beyond the simplicity of the threshold model. This model uses the bulk image or PSF as described in Sec. 3.1. The consideration of specific exposure kinetics, the diffusion of a single species, and the characteristic development of the processed photopolymer reflects our understanding of photopolymer behavior. Next, we will describe the specifics of each processing step in the compact model, clarifying how this method provides a more comprehensive depiction of TPP than its threshold counterpart.

**Exposure reactions.** The Dill equations³⁴ describe single-photon exposure and are given by Eqs. (3)–(4). In these equations, α is the absorption coefficient, [PAC] represents the photosensitive component, A<sub>Dill</sub> is the bleaching absorption coefficient, B<sub>Dill</sub> is the non-bleaching absorption coefficient, and C<sub>Dill</sub> defines the photosensitivity of the resist.
$$ \alpha = A_{\text{Dill}}[\text{PAC}] + B_{\text{Dill}}, \qquad (3) $$

$$ \frac{\partial[\text{PAC}]}{\partial t} = -C_{\text{Dill}}I[\text{PAC}]. \qquad (4) $$

The TPP exposure model adapts these equations to account for the nonlinear optical phenomena inherent to TPP. The likelihood of two-photon absorption is proportional to the square of the light intensity ($I^2$). Hence, in our model, the photoresist exposure transforms the intensity distribution from the optical model by squaring the intensity values to properly simulate the two-photon absorption process. To accommodate the characteristics of the TPP photoresist, it is presumed that the absorption coefficient remains invariant throughout the exposure. This assumption negates the need for a distinction between the Dill coefficients A and B ($A_{\text{Dill}}$ and $B_{\text{Dill}}$), leading to a simplification where $A_{\text{Dill}} = 0$, due to the absence of bleaching effects. To encapsulate these concepts, the modified equations for the TPP exposure reaction can be described as follows:

$$ [\text{TPPI}] = \exp(-I^2 \text{DC}_{\text{Dill}}). \qquad (5) $$

The concentration of the two-photon photoinitiator ([TPPI]) decreases exponentially with the product of a constant absorption coefficient ($C_{\text{Dill}}$), the exposure dose ($D$), and the squared intensity distribution from the optical model ($I^2$)

**Dark phase.** In the context of TPP, the “dark phase” is a term adopted to describe the period following laser exposure, where critical polymerization processes take place. Unlike the thermal post-exposure bake (PEB) in conventional photolithography, the TPP dark phase does not involve heat-driven diffusion. Instead, this phase is characterized by the diffusion of radicals and other species, which occurs due to concentration gradients established during the laser-induced polymerization. Fickian diffusion, which traditionally describes the thermal spreading of chemical species in materials science, is employed here to model the migration of light-induced reactive radicals, generated by TPPI, within the TPP resist

$$ \frac{\partial[\text{TPPI}]}{\partial t} = D\nabla^2[\text{TPPI}]. \qquad (6) $$

Here, $D$ represents the diffusion coefficient. The Fickian diffusion equation describes the isotropic spreading of radicals, which occurs as a result of spatial concentration gradients rather than thermal energy. Although the model assumes the diffusion of a single reactive species—a simplification compared with the multiple species potentially present—it offers a significant enhancement over the threshold model by capturing the resist’s response to incident light through polymerization.

Next, we will explore the development process, which forms the microstructures by removing the unexposed and unpolymerized material and finalizing the structure.

**Development.** Development involves the interaction of the resist and a developer to reveal the three-dimensional resist profile. In the conventional positive-tone development, the rate of photoresist removal is described by the Mack model,³⁵ which establishes a relationship between the local density of exposed sites and the development rate. The model is formalized by

$$ R = \frac{R_{\text{max}}(a + 1)(1 - M)^N}{a + (1 - M)^N} + R_{\text{min}}. \qquad (7) $$

Representing $a$ in terms of another parameter $M_{\text{th}}$ gives the following equation:

$$ a = \frac{(N + 1)}{(N - 1)} (1 - M_{\text{th}})^N, \qquad (8) $$

where $M$ represents the relative polymer concentration, that is the polymer concentration divided by the maximum polymer concentration. $R_{\text{max}}$ is the development rate for the unexposed regions of the resist; $R_{\text{min}}$ is the development rate function for the fully exposed (or fully polymerized) parts of the resist. $M_{\text{th}}$ is the threshold relative polymer concentration. $N$ defines the degree of
selectivity of the developer; the higher the value of N, the more the developer behaves similar to a threshold filter between exposed and unexposed regions.

Typical TPP materials are negative-tone photoresists. These resists, unlike their positive-tone counterparts, develop the unexposed areas rather than the exposed ones. As a result, the standard model is adapted into an inverse form to better represent the specific chemical development process in TPP, as expressed in³⁶

$$R_{inv} = R_{max} - R_{orig}. \qquad (9)$$

The development rates are used as inputs to a fast marching algorithm,³⁷ which evaluates the topology of the developed pattern. This algorithm calculates the local developer arrival time (DART) at every position in the resist by moving the developer front (the boundary between the developer and resist) with a speed proportional to the local development rates. The fast marching algorithm is a particular case of the level set method when the front propagates monotonically. Spatial positions with developer arrival times less than the development time are washed away after development. The part of the resist that remains after development forms the fabricated profile. The final resist profile is extracted along the iso-contour of the DART at this specific development time, capturing the precise boundary of the developed structure. Furthermore, during development, significant shrinkage may occur due to variations in cross-linking densities within the polymerized regions, which affects the final dimensions and fidelity of the printed structures.³⁸ Shrinkage was not explicitly modeled in this study because the custom-made resist used in the experiments showed no observable shrinkage.

In summary, our compact model provides a simple yet effective approach to simulate TPP. It integrates optical and photoresist components, enhancing the threshold model by incorporating the diffusion of a single substance within the resistor. In addition, the inclusion of a development process improves the model's resemblance to real-world TPP processes. Figure 3 illustrates the progression from radical distribution to the final 3D shape of the photopolymer after development. Although the compact model offers faster simulations compared with the more comprehensive full polymerization model, it struggles to capture the details of two-photon absorption, especially when multiple substances diffuse simultaneously. To address this limitation, we will explore a more comprehensive model in the subsequent section.

### 3.2.3 Full Model of Polymerization

The full model of polymerization represents a recently developed approach for TPP modeling.³⁹ The model uses the same bulk image as the threshold model and the compact model, but it provides a more sophisticated treatment of polymerization. Instead of relying on the single-component diffusion approach of the above-described compact model, it incorporates real two-photon processes, allowing for a more accurate representation of polymerization dynamics. This includes the diffusion and kinetics of multiple species, providing a more realistic description of chemical processes during fabrication.

Next, we will explore each step of the full model in detail, starting with the exposure. By explaining the specifics of this new method, we aim to provide a comprehensive understanding of TPP fabrication, overcoming the limitations of previous models.

<figure><img src="image_4.png" /><figcaption>Fig. 3 Simulation flow through the steps of TPP using the compact model.</figcaption></figure>
**Exposure reactions.** Exposure leads to the creation of radicals ($R_*$) from the photosensitive component. According to

$$[R] = [PSC]_0 - [PSC]_0 \exp(-0.5 \times \delta_{TPA} \phi \tilde{N}(x, y, z)^2 t), \quad (10)$$

where $PSC_0$ is the initial concentration of the photosensitive material in the resist and $\delta_{TPA}$ is the two-photon absorption cross-section. Two-photon absorption systems typically employ pulsed lasers due to their smaller cross-section, resulting in very high peak intensities $\tilde{N}(x, y, z)$. The photon flux $\tilde{N}(x, y, z)$ relates to the pulse duration ($\tau$) and laser repetition rate ($R_{rep}$) by

$$\tilde{N}(x, y, z) = \frac{I(x, y, z)}{\frac{hc}{\lambda} \tau R_{rep}}, \quad (11)$$

where $h$ is the Planck's constant and $c$ is the speed of light in vacuum. Figure 4 shows the generated photo-radicals for a two-photon absorption process. The generated radicals initiate polymerization reactions in the resist.

In addition to the generation of radicals, a focused light source locally raises the temperature in the resist. The extent of this local temperature increase depends on the intensity, exposure duration, and material properties of the resist: heat capacity, density, and enthalpy of polymerization. This thermal effect is mathematically described using the heat equation:40

$$\frac{\partial T}{\partial t} = \frac{\kappa}{\rho c_p} \nabla^2 T + \frac{\alpha I[M]}{\rho c_p} + \frac{H}{\rho c_p} \frac{\partial [M]}{\partial t}. \quad (12)$$

Here, $[M]$ is the concentration of monomers, $H$ is the reaction enthalpy, $I[M]$ is the rate of heat generation due to a photochemical reaction, $\kappa$ represents the thermal conductivity of the resist, $\rho$ denotes its density, $c_p$ is the specific heat capacity, and $\alpha$ stands for the molar absorption coefficient of the monomer. The solution of the heat equation yields a temperature distribution within the resist, leading to non-uniform diffusion and reaction rates. These rates are commonly described by the Arrhenius equation.41 The temperature profile at the focal region, as depicted in Fig. 4, resembles the intensity distribution of the focused light beam.

**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 model42

$$\frac{\partial [R]}{\partial t} = D_R \nabla^2 [R] - (2k_t [R]^2 + k_q [Q][R] + k_p [M][R]), \quad (13)$$

$$\frac{\partial [Q]}{\partial t} = D_Q \nabla^2 [Q] - (k_q [Q][R] + k_q [M^*][Q]), \quad (14)$$

<figure><img src="image_5.png" /></figure>

Optical parameters:

Wavelength = 0.522 µm

NA = 1.4

Immersion refractive index = 1.518

Exposure time = 0.00001 sec

Power = 0.012 W

<caption>**Fig. 4** Intensity distribution, temperature profile, and radical distribution. The plot in the middle represents the intensity distribution within the resist, whereas the figure on the right illustrates the temperature profile after exposure. The plot on the right illustrates the spatial distribution of photo-radicals.</caption>
$$ \frac{\partial[M]}{\partial t} = D_M \nabla^2[M] - \left( k_i[M][R] + \oint k_p[M]^*(t, x')(M)(t, x')G(x, x')dx' \right), \quad (15) $$

$$ \frac{\partial[M^*]}{\partial t} = D_M^* \nabla^2[M^*] + k_i[M][R] - (2k_i[M^*]^2 + k_{tp}[M^*] + k_q[M^*][Q]), \quad (16) $$

where the term R represents the amount of radicals, M is the amount of monomers, $[M^*]$ is the amount of growing active monomer chains, Q is the amount of quencher, $D_M$ is monomer diffusion coefficient, $D_Q$ is quencher diffusion coefficient, $k_i$ is the rate constant for termination, $k_p$ is the rate constant for propagation, $k_{tp}$ is the rate of termination by primary radicals, $k_q$ is the rate constant for quenching, and $k_i$ is the rate constant for initiation. A convolution is utilized in Eq. (15) 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 Eq. (17):

$$ G(x, x') = \frac{1}{\sqrt{2\pi\sigma}} \exp\left(-\frac{(x-x')^2}{2\sigma^2}\right). \quad (17) $$

The term $\sigma$ in Eq. (17) is a measure of the spatial extent of the propagation reaction that results from a single radical. The extent of this spread varies with the type of resist material used. In our current implementation, diffusion constants in Eqs. (13)–(16) are treated as fixed parameters to simplify the modeling. However, we acknowledge that, in practice, these diffusion constants may vary both spatially and temporally due to local polymerization effects and cross-linking-induced changes. The incorporation of such dynamic variations represents an important direction for future model enhancements. At 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, whereas 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, monomers, polymers, and quenchers after the dark phase.

**Development.** For the development step, we utilize the Mack rate model and fast marching algorithm as described in Sec. 3.2.2.

<figure><img src="image_6.png" /></figure>

**Optical parameters:**

Wavelength = 0.522 µm

NA = 1.4

Immersion refractive index = 1.518

Exposure time = 0.00001 sec

Power = 0.012 W

Repetition rate = 63 MHz

Pulse duration = 250 fs

**Resist parameters:**

Resist refractive index = 1.525

Substrate refractive index = 1.47

$C_{DII} = 0.95 \text{ cm}^2/\text{mJ}$

Radical diffusion ($\alpha$) = 0.054 µm

Quencher loading = 4.3e-5 mol/cm$^3$

Development time = 120 sec

$R_{max} = 0.4 \text{ µm/sec}$

$R_{min} = 0.0001 \text{ µm/sec}$

M = 0.8

N = 25

<caption>**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.</caption>
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. Unlike the compact model that simulates the diffusion of just a single species, this model captures the interplay between different species such as monomers, polymers, radicals, and quenchers. The quenching of the photopolymerization reaction, typically by oxygen, provides an important component of the model with a significant impact on the resulting voxel shape. Although this model offers a more detailed simulation, capturing the complex interactions during TPP, it results in larger computation times and more, mostly unknown, model parameters. Typical simulation times for a single voxel (excluding optical simulations) are below 1 s for the threshold model, ~5 s for the compact model, and around 35 s for the full polymerization model. For larger structures created via SLM projection methods, optimized numerical algorithms enable concurrent computation of multiple voxels, typically requiring around 10 to 12 s for simulating a 15 × 15 µm structure using the compact model.

# 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. In the following examples, the light is assumed to be incident from above, propagating downward, and the interface between the immersion fluid and the photoresist is situated in the upper region of the imaging volume.

Starting with imaging, the NA of the imaging system plays an important role in determining the feature size 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 Fig. 6, where different NA values show a clear correlation with voxel dimensions. The top row of images illustrates how the NA affects the 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

Imaging - NA

<figure><img src="image_7.png" /></figure>

**Model:** Full Model of Polymerization

**Optical parameters:**
Wavelength = 0.522 µm
Immersion refractive index = 1.518
Exposure time = 0.00001 sec
Power = 0.012 W
Repetition rate = 63 MHz
Pulse duration = 250 fs

**Resist parameters:**
Resist refractive index = 1.525
Substrate refractive index = 1.47
$C_{p00}$ = 0.95 cm²/mJ
Radical diffusion (σ) = 0.054 µm
Quencher loading = 4.3e-5 mol/cm²
Development time = 120 sec
$R_{m00}$ = 0.4 µm/sec
$R_{m00}$ = 0.0001 µm/sec
M = 0.8
N = 25

**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 the size of the voxel's 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.
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 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 7 demonstrates 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 pronounced voxel elongation. The contours in the lower right of Fig. 7 suggest an increasing (top-down) asymmetry of the voxel for a 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 Fig. 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 $\sigma$ plays a significant role in the diffusion process, which critically impacts the development of the voxel. The term $\sigma$ 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 [Eq. (17)]. Figure 9 illustrates how different values of $\sigma$ affect the voxel. With smaller $\sigma$ values, indicating limited diffusion, the resulting voxels are more compact, as visualized in the DART images. As $\sigma$ 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 $\sigma$. The left graph plots an increasing trend in both the diameter and length of voxels with rising $\sigma$ values, whereas

<figure><img src="image_8.png" /><figcaption>Fig. 7 Impact of photoresist refractive index on voxel characteristics.</figcaption></figure>
Exposure - power

<figure><img src="image_9.png" /><figcaption>Fig. 8 Impact of power on voxel characteristics.</figcaption></figure>

**Model:**

Full Model of Polymerization

**Optical parameters:**

Wavelength = 0.522 µm

NA = 1.4

Immersion refractive index = 1.518

Exposure time = 0.00001 sec

Repetition rate = 63 MHz

Pulse duration = 250 fs

**Resist parameters:**

Resist refractive index = 1.525

Substrate refractive index = 1.47

$C_{DII} = 0.95 \text{ cm}^3/\text{mJ}$

Radical diffusion ($\sigma$) = 0.054 µm

Quencher loading = 4.3e-5 mol/cm³

Development time = 120 sec

$R_{max}$ = 0.4 µm/sec

$R_{min}$ = 0.0001 µm/sec

M = 0.8

N = 25

DarkPhase - Sigma ($\sigma$)

<figure><img src="image_10.png" /><figcaption>Fig. 9 Impact of $\sigma$ on voxel characteristics.</figcaption></figure>

**Model:**

Full Model of Polymerization

**Optical parameters:**

Wavelength = 0.522 µm

NA = 1.4

Immersion refractive index = 1.518

Exposure time = 0.00001 sec

Power = 0.012 W

Repetition rate = 63 MHz

Pulse duration = 250 fs

**Resist parameters:**

Resist refractive index = 1.525

Substrate refractive index = 1.47

$C_{DII} = 0.95 \text{ cm}^3/\text{mJ}$

Quencher loading = 4.3e-5 mol/cm³

Development time = 120 sec

$R_{max}$ = 0.4 µm/sec

$R_{min}$ = 0.0001 µm/sec

M = 0.8

N = 25

the right graph provides a contour representation of the voxel shape change in response to different $\sigma$ values. Here, larger $\sigma$ 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 multiple 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.

## 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 (Sec. 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.¹⁵,⁴³

In summary, although 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 material 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 diffusion⁴⁴ 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 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 (Sec. 3). Calibration results [Fig. 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.

<figure><img src="image_11.png" /><figcaption>(b)</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>(c)</figcaption></figure>

<figure><img src="image_13.png" /></figure>

**Optical parameters:**

Wavelength = 0.522 µm

NA = 1.4

Immersion refractive index = 1.518

**Resist parameters:**

Resist refractive index = 1.525

Substrate refractive index = 1.47

$C_{sub}$ = 0.95 cm²/mJ

Radical diffusion (a) = 0.054 µm

Quencher loading = 4.3e-5 mol/cm³

Development time = 120 sec

$R_{max}$ = 0.4 µm/sec

$R_{min}$ = 0.0001 µm/sec

M = 0.8

N = 25

**Fig. 10** 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.
The full model of polymerization includes a critical radical control parameter that governs the interplay between diffusion kinetics and exposure time. As shown in Fig. 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. Although 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 projection lens. Although the compact model and full model of polymerization provide significant insights into voxel size prediction, there remain areas where the fit between the model and experimental data can be improved. One possible factor influencing the model’s accuracy is the filling factor of the Gaussian illumination beam. Initial discrepancies between the model predictions and experimental data could be attributed to the assumption that the pupil of the projection lens is illuminated by constant intensity over the complete pupil plane. In reality, the pupil is illuminated by a Gaussian-shaped intensity distribution. The filling factor is a number that defines the ratio between the objective aperture radius and the beam waist radius. A filling factor of 1 corresponds to a fully filled, flat beam profile, whereas values less than 1 indicate underfilling, leading to a different illumination profile and a lower effective NA. When the pupil is underfilled, the effective NA is reduced, leading to inaccuracies in voxel size predictions, especially in high-power settings.

We tested the effect of different filling factors on voxel dimensions, as shown in Fig. 12, where voxel sizes are compared for filling factors of 1 and 0.82. These simulations demonstrate that the pupil filling can be used to control the shape of the voxel and adapt it for certain applications. To improve the model’s accuracy, the filling factor was incorporated. This adjustment significantly improved the calibration results, as shown in Fig. 13, which compares voxel sizes for different laser powers. The incorporation of the filling factor corrected the earlier overestimation of voxel diameters, resulting in a much closer fit to experimental data. The filling factor adjustment allowed for better predictions across a range of exposure times and power levels. More advanced techniques for manipulating voxel shape by adjusting the intensity and phase

<figure><img src="image_14.png" /><figcaption>**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)</figcaption></figure>
<figure><img src="image_15.png" /><figcaption>**Fig. 12** Comparison of bulk images at different filling factors. Left: bulk image at a filling factor of 1.0. Middle: bulk image at a filling factor of 0.9. Right: bulk image at a filling factor of 0.8, showing the gradual changes in intensity distribution with a decreasing filling factor.</figcaption></figure>

<figure><img src="image_16.png" /></figure>

**Optical parameters:**

Wavelength = 0.522 µm

NA = 1.4

Immersion refractive index = 1.518

Filling factor = 0.85 and 0.82

**Resist parameters:**

Resist refractive index = 1.525

Substrate refractive index = 1.47

$C_{\text{DII}}$ = 0.95 cm²/mJ

Radical diffusion (a) = 0.054 µm

Quencher loading = 4.3e-5 mol/cm³

Development time = 120 sec

$R_{\text{max}}$ = 0.4 µm/sec

$R_{\text{min}}$ = 0.0001 µm/sec

M = 0.8

N = 25

<caption>**Fig. 13** Left graph shows how voxel size changes with laser power in the compact model (filling factor = 0.85). Right graph shows how voxel size changes with laser power in the full model of polymerization (filling factor = 0.82).</caption>

distribution of the illumination in the objective lens have been demonstrated, for example, by Schmidt et al.,⁴⁵ Waller et al.,⁴⁶ Tičkūnas et al.,⁴⁷ and Sun et al.⁴⁸

In this study, we evaluated three modeling approaches for two-photon lithography, with the full model providing the most accurate predictions by incorporating diffusion, radical dynamics, and exposure time effects. Adjusting the pupil filling factor improved the accuracy of models, aligning their predictions more closely with experimental data.

### 4.2.2 Application of the model to alternative exposure strategies

In this section, our focus shifts to adapting the compact model to interpret and simulate the results from a specialized two-photon lithography setup used by IMTA. This setup utilizes a spatial light modulator (SLM) to precisely control light patterns for photopolymerization processes.

The IMTA's experimental setup⁴⁴ consists of a sequence of components designed to manipulate laser light to achieve precise polymerization (Fig. 14). Starting with a laser diode as the light source, the light is homogenized and shaped using a combination of optical elements including polarizing cubes, collimation and tube lenses, and a diffractive optical element. The light is modulated by an SLM. The microscope objective projects the desired patterns onto the photoresist. This setup permits the creation of not only individual voxels but also structured patterns.

The triplet-triplet annihilation (TTA) photoresist used in this process is sensitive to continuous wave (CW) lasers and is pivotal for achieving high plot rates while maintaining 3D microfabrication capabilities. Although pulsed lasers are traditionally used for two-photon lithography, several authors have demonstrated 3D volumetric fabrication with CW laser sources instead.⁴⁹ The TTA process is characterized by its lower light intensity requirement and its compatibility with non-coherent light sources, as described by Limberg et al.⁵⁰
<figure><img src="image_17.png" /><figcaption>Fig. 14 Schematic of the two-photon lithography setup utilized by IMTA France.⁴⁴</figcaption></figure>

Using experimental data from IMTA, which included various structure widths and levels of defocus, we employed the compact model to simulate the lithography process. The optimization focused on the shape of the voxels, with adjustments made to both optical and resist parameters in the model. These adjustments were aimed at aligning the simulation results with the experimental profiles (Fig. 15). The results obtained demonstrate the model's ability to accurately describe pattern formation, even for more complex layouts.

Figure 15 provides a visual comparison between the experimental data and the simulated profiles obtained with the compact model, showcasing how the model can be tuned to replicate the actual fabricated structures. Although these are 2.5D structures, the goal moving forward is to extend this modeling capability to fully 3D structures. The integration of the compact model with IMTA's setup marks a significant step towards achieving this goal, providing a predictive tool that can adapt to various fabrication scenarios.

## 4.3 First Applications of Models Toward Metasurfaces

The classical design of metasurfaces neglects fabrication effects and assumes an ideal building block with sharp edges and corners. Optical and chemical effects during lithographic fabrication lead to rounded features and proximity effects, revealing the complex aspects of the manufacturing process. The manufacturable design of functional metasurfaces requires predictive modeling of the optical and photochemical effects in two-photon lithography.

FORTH has developed initial metasurface designs [Fig. 16(a)],⁵¹ which present idealized structures intended to manipulate incident light with high precision. However, the transition from concept to fabrication necessitates a design evolution to ensure mechanical stability and manufacturability [Fig. 16(b)]. Designs were modified, replacing sharp edges with smoother, more rounded shapes such as cones and half-ellipsoids. These shapes mimic the output of the real fabrication process.

We used the full model of polymerization to simulate the steps involved in standard two-photon fabrication. To construct these pillars, the simulation involves stacking multiple voxels,
<figure><img src="image_18.png" /><figcaption>**Fig. 15** Overlay of simulated and experimental photopolymerization results for different focus settings. The intensity distribution mask applied in the photopolymerization process is shown at the top. Below, the corresponding microstructures created at defocus levels from +10 to -10 µm are presented, with simulated data from the compact model superimposed (yellow areas) on the experimentally measured gray-tone SEM images. The separate plots of experimental (gray-tone SEM images) and simulated (yellow-marked areas) data on the right show a zoomed-in view of the area highlighted by the blue box, providing a clearer comparison of experimental and simulated data. This visual comparison highlights the compact model's capability to replicate the actual fabrication process across different focal plane adjustments for 2.5D structures.</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>**Fig. 16** Supercell of an optical metasurface for beam steering applications. The design in panel (a) exploits vertical waveguide segments (pillars with circular cross-sections), which supply transmission phase levels of 120 deg, 240 deg, and 360 deg. The extent of the supercell is 1350 nm, and it is designed to deflect a normally incident plane wave with a wavelength of 750 nm into an angle of 33.75 deg (1st diffraction order). (b) The design is modified to avoid sharp edges, allow for better mechanical stability, and make it friendlier to fabrication with two-photon lithography. (c) Pillars simulated by full polymerization model with the help of Dr. LiTHO lithography simulator. To generate the highest pillars, several voxels were placed on top of each other. Designs (a), (b): FORTH.</figcaption></figure>

with the calibration necessary to achieve the correct voxel size and shape. Through an optimization algorithm, we adjusted the shape and distance of the voxels to ensure that the simulated structures aligned with the desired specifications [Fig. 16(c)]. The demonstrated modeling strategy can bridge the gap between the theoretical designs of metasurfaces by FORTH and their physical realizations. The models provide a clearer expectation of the outcomes from lithographic systems, offering a more realistic framework within which to refine the designs for actual fabrication conditions. This approach illustrates a practical application of our models. It provides a first step to align the expected metasurface performance with achievable manufacturing processes.
<figure><img src="image_20.png" /><figcaption>**Fig. 17** Schematic overview of the double exposure photolithography process for generating a 3D freestanding T-shaped structure. The top-left panel shows the target 3D geometry. The next two panels display the top-view transmission patterns of the first and second masks, respectively. These masks are exposed sequentially at different focal depths. The lower row shows x-z cross-sections of the simulated intensity distributions for each individual exposure and their summation. Red-dashed lines indicate the focal planes for the first and second exposures. The final panel (top-right) illustrates the resulting 3D structure formed in the photoresist. The process highlights how spatially separated exposures in the z-direction can be combined to enable volumetric patterning beyond traditional 2D lithography. Optical and resist parameters used in the simulation are listed on the bottom left.</figcaption></figure>

Expanding the utility of the compact model previously outlined, we explore its application in simulating double exposure processes as conducted by IMTA.⁵² Through this technique, they craft complex metasurfaces by layering exposures with varying masks. The compact model, originally developed for single-exposure simulations, has been adapted to handle double-exposure methodology. The top-left panel of Fig. 17 depicts the target 3D geometry for a T-shaped structure, serving as a first simulation example of this double exposure approach. The middle-top panels display the two mask transmission patterns utilized in the double-exposure lithography process. These two masks are projected, one after the other, at distinct focal depths to form the “T”-shaped structure in the resist. The lower row shows x-z cross-sections of the simulated intensity distributions for each individual mask inside the photoresist and their summation, whereas the final top-right panel presents the resulting 3D structure formed in the photoresist. Optical and resist parameters used in the simulation are listed on the bottom-left side of the figure. This example highlights how the double-exposure method can achieve volumetric patterning beyond conventional single-layer lithography.

## 5 Conclusion and Outlook

In this comprehensive study, we have embarked on a detailed journey to enhance the modeling and fabrication of optical metasurfaces through two-photon lithography. We have introduced three 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 15 mW for a 10-µs 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 of reliably predicting voxel lengths.

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 versus 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 Fig. 16 and of the formation of a T-shaped pattern by a double exposure demonstrate the first steps toward process-aware design and fabrication of optical metasurfaces.

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. We will also investigate temporal effects, such as crosslinking-induced changes in diffusion during sequential voxel exposures, to improve the realism of our models. 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.

## Disclosures

The authors declare there are no financial interests, commercial affiliations, or other potential conflicts of interest that have influenced the objectivity of this research or the writing of this paper.

## Code and Data Availability

Fraunhofer IISB, HIMT, and IMTA proprietary information will not be made available, but the manuscript content is consistent with JM3 technical content guidelines.

## Acknowledgments

This project has received funding from the European Union's Horizon Europe research and innovation program (Grant 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 et al., "Multifunctional metaoptics based on bilayer metasurfaces," *Light: Sci. Appl.* **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)," *Proc. SPIE* **13023**, 1302306 (2024).

4. J. Zeng et al., "Broadband and high-efficiency multi-tasking silicon-based geometric-phase metasurfaces: a review," *Photonics* **9**(9), 606 (2022).

5. S. Zahra et al., "Electromagnetic metasurfaces and reconfigurable metasurfaces: a review," *Front. Phys.* **8**, 593411 (2020).

6. G. Zyla et al., "Two-photon polymerization as a potential manufacturing tool for biomimetic engineering of complex structures found in nature," *J. Opt. Microsyst.* **2**(3), 031203 (2022).

7. E. Skliutas et al., "X-photon laser direct write 3D nanolithography," *Virt. Phys. Prototyp.* **18**(1), e2228324 (2023).

8. V.-C. Su et al., "Advances in optical metasurfaces: fabrication and applications [invited]," *Opt. Express* **26**(10), 13148–13182 (2018).

9. M. Malinauskas et al., "A femtosecond laser-induced two-photon photopolymerization technique for structuring microlenses," *J. Opt.* **12**(3), 035204 (2010).

10. C. Arnoux et al., "Understanding and overcoming proximity effects in multi-spot two-photon direct laser writing," *Addit. Manuf.* **49**, 102491 (2022).

11. L. A. Pérez Covarrubias et al., "Proximity effect in parallelized microfabrication using two-photon polymerization," *Proc. SPIE* **11349**, 1134900 (2020).

12. S. K. Saha et al., "Effect of proximity of features on the damage threshold during submicron additive manufacturing via two-photon polymerization," *J. Micro Nano-Manuf.* **5**, 031002 (2017).

13. A. S. van de Nes et al., "Calculation of the vectorial field distribution in a stratified focal region of a high numerical aperture imaging system," *Opt. Express* **12**, 1281–1293 (2004).

14. P. Török et al., "Electromagnetic diffraction of light focused through a planar interface between materials of mismatched refractive indices: an integral representation," *J. Opt. Soc. Am. A* **12**, 325–332 (1995).

15. R. J. DeVoe et al., "Voxel shapes in two-photon microfabrication," *Proc. SPIE* **4797**, 310–316 (2003).

16. S. Kuebler and M. Rumi, "Nonlinear optics, applications | three-dimensional microfabrication," in *Encyclopedia of Modern Optics*, R. D. Guenther, Ed., pp. 189–206, Elsevier, Oxford (2005).

17. C. N. LaFratta et al., "Multiphoton fabrication," *Angew. Chem.* **46**, 6238–6258 (2007).

18. K. Lee et al., "Recent developments in the use of two-photon polymerization in precise 2D and 3D micro-fabrication," *Polym. Adv. Technol.* **17**, 72–82 (2006).

19. P. Kiefer et al., "Sensitive photoresists for rapid multiphoton 3D laser micro- and nanoprinting," *Adv. Opt. Mater.* **8**(19), 2000895 (2020).

20. J. E. Johnson, Y. Chen, and X. Xu, "Model for polymerization and self-deactivation in two-photon nanolithography," *Opt. Express* **30**, 26824–26840 (2022).

21. P. Somers et al., "The physics of 3D printing with light," *Nat. Rev. Phys.* **6**, 99–113 (2023).

22. J. Serbin et al., "Femtosecond laser-induced two-photon polymerization of inorganic-organic hybrid materials for applications in photonics," *Opt. Lett.* **28**(5), 301–303 (2003).

23. 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, pp. 45–61, William Andrew Publishing, Oxford (2016).

24. J. B. Mueller, J. Fischer, and M. Wegener, "Chapter 3.1 - Reaction mechanisms and in situ process diagnostics," in *Three-Dimensional Microfabrication Using Two-photon Polymerization*, T. Baldacchini, Ed., Micro and Nano Technologies, 2nd ed., pp. 175–196, William Andrew Publishing (2020).

25. R. Pingali and S. K. Saha, "Printability prediction in projection two-photon lithography via machine learning based surrogate modeling of photopolymerization," *J. Micro Nano-Manuf.* **10**, 031005 (2023).

26. J.-F. Xing et al., "Improving spatial resolution of two-photon microfabrication by using photoinitiator with high initiating efficiency," *Appl. Phys. Lett.* **90**(13), 131106 (2007).

27. V. Sedova et al., "Advances in modeling and optimization for two-photon lithography," *Proc. SPIE* **13023**, 1302309 (2024).

28. D. Ladika et al., "X-photon 3D lithography by fs-oscillators: wavelength-independent and photoinitiator-free," *Light: Adv. Manuf.* **5**(4), 567–579 (2024).

29. A. Diaspro, *Confocal and Two-Photon Microscopy: Foundations, Applications, and Advances*, Wiley-Liss, New York (2002).

30. T. Fühner et al., "Dr.LITHO: a development and research lithography simulator," *Proc. SPIE* **6520**, 65203F (2007).

31. A. Erdmann, *Optical and EUV Lithography: A Modeling Perspective*, SPIE Press, Bellingham, Washington (2021).
32. H.-B. Sun et al., "Experimental investigation of single voxels for laser nanofabrication via two-photon photo-polymerization," *Appl. Phys. Lett.* **83**, 819-821 (2003).

33. D. Fuard, M. Besacier, and P. Schiavone, "Assessment of different simplified resist models," *Proc SPIE* **4691**, 1266-1277 (2002).

34. F. H. Dill et al., "Characterization of positive photoresist," *IEEE Trans. Electron Devices* **22**, 445 (1975).

35. C. A. Mack, "New kinetic model for resist dissolution," *J. Electrochem. Soc.* **139**(4), L35-L37 (1992).

36. W. Gao et al., "Application of an inverse Mack model for negative tone development simulation," *Proc. SPIE* **7973**, 79732W (2011).

37. J. A. Sethian, "Fast marching level set methods for three-dimensional photolithography development," *Proc. SPIE* **2726**, 262-272 (1996).

38. N. Lang et al., "Towards efficient structure prediction and pre-compensation in multi-photon lithography," *Opt. Express* **30**, 28805-28816 (2022).

39. T. P. Onanuga, "Process modeling of two-photon and grayscale laser direct-write lithography," Dissertation, Friedrich-Alexander-Universität Erlangen-Nürnberg, Erlangen, Germany, Doktor-Ingenieur (2019).

40. J. B. Mueller et al., "In-situ local temperature measurement during three-dimensional direct laser writing," *Appl. Phys. Lett.* **103**(12), 123107 (2013).

41. A. K. O'Brien and C. N. Bowman, "Modeling the effect of oxygen on photopolymerization kinetics," *Macromol. Theory Simul.* **15**(2), 176-182 (2006).

42. M. R. Glee son, J. Guo, and J. T. Sheridan, "Recent developments in the NPDD model," *Proc. SPIE* **8429**, 84291C (2012).

43. S. Wu, J. Serbin, and M. Gu, "Two-photon polymerisation for three-dimensional micro-fabrication," *J. Photochem. Photobiol. A: Chem.* **181**(1), 1-11 (2006).

44. F. Ogor et al., "Modelling and compensating proximity effects in massively parallelized multi-photon photo-plotting," *Proc. SPIE* **12995**, 1299509 (2024).

45. V. Schmidt et al., "Two-photon 3D lithography: a versatile fabrication method for complex 3D shapes and optical interconnects within the scope of innovative industrial applications," *J. Laser Micro/Nanoeng.* **2**, 16 (2007).

46. E. H. Waller, M. Renner, and G. von Freymann, "Active aberration- and point-spread-function control in direct laser writing," *Opt. Express* **20**, 24949-24956 (2012).

47. T. Tičkūnas, D. Paipulas, and V. Purlys, "Dynamic voxel size tuning for direct laser writing," *Opt. Mater. Express* **10**, 1432-1439 (2020).

48. M. Sun et al., "Fast, precise, high contrast laser writing for photonic chips with phase aberrations," *Laser Photonics Rev.* **18**, 2300702 (2023).

49. M. Thiel et al., "Direct laser writing of three-dimensional submicron structures using a continuous-wave laser at 532 nm," *Appl. Phys. Lett.* **97**(22), 221102 (2010).

50. D. K. Limberg, J.-H. Kang, and R. C. Hayward, "Triplet-triplet annihilation photopolymerization for high-resolution 3D printing," *J. Am. Chem. Soc.* **143**(9), 3387-3395 (2021).

51. G. Perrakis, M. Kafesaki, and O. Tsilipakos, "Optical metasurfaces with two-photon lithography: design considerations for beam steering applications," *Proc. SPIE* **13023**, 1302307 (2024).

52. J. Rovera et al., "Digital modelling of a massively parallelised multiphoton polymerisation plot process," *Proc. SPIE* **PC13023**, PC1302309 (2024).

**Valeriia Sedova** is a PhD student at Friedrich-Alexander University Erlangen-Nürnberg (FAU) and a researcher in the Computational Lithography and Optics Group at Fraunhofer IISB. Her research focuses on modeling and optimization of lithographic processes, including grayscale and two-photon lithography, as well as the development of inverse lithographic techniques. She received her MSc in Advanced Optical Technologies from FAU. She is the recipient of the Best Paper Award at SPIE Optical Systems Design 2024 and the Applied Photonics Award for Best Master's Thesis 2023. She is also a member of SPIE.

**Florie Ogor** studied engineering at Télécom Physique Strasbourg before pursuing a PhD in optics at IMT Atlantique. Her doctoral research focused on massively parallelized multiphoton 3D microfabrication. She is currently working as an optical engineer at Cristalens.

**Joël Rovera** is a PhD student in photonics at the Optics Department of IMT Atlantique in Brest, France. His research focuses on 3D microfabrication using two-photon polymerization (2PP), with a particular emphasis on process simulations for precompensating fabrication instructions. He received his master's degree in optics and image processing from the University of Saint-Étienne and an engineering degree in optical engineering from the Institut d'Optique Graduate School (IOGS), completed while working for a leading aerospace company.
**Odysseas Tsilipakos** obtained his PhD from the School of Electrical and Computer Engineering, Aristotle University of Thessaloniki (AUTH), in 2013. In 2016, he joined the Institute of Electronic Structure and Laser (IESL) of the Foundation for Research and Technology Hellas (FORTH) as a postdoctoral fellow. Since October 2022, he holds a researcher position (assistant Prof. level) at the Theoretical and Physical Chemistry Institute (TPCI) of the National Hellenic Research Foundation (NHRF), where he leads the activity on “Theoretical and Numerical Methods for Photonics, Optoelectronics and Metamaterials.” He is an author/co-author of 70 publications in refereed journals and 95 contributions to international conferences. He has received the Best Student Paper Award in SPIE Photonics Europe 2012. He is a senior member of IEEE and a member of Optica.

**Jonas Wiedenmann** is currently the head of two-photon polymerization (TPP) development at Heidelberg Instruments Mikrotechnik GmbH in Würzburg. He holds a master of physics degree in optoelectronics and laser from Heriot-Watt University in Scotland and his MSc in physics from the University of Würzburg, Germany. From 2014 to 2019, he was a research assistant at the Chair of Prof. Molenkamp, where he worked on the fabrication and analysis of II–VI semiconductor–superconductor hybrid nanostructures. In his current role, he draws on deep expertise in direct laser writing lithography to develop and optimize multiphoton 3D printing systems and workflows.

**Kevin Heggarty** studied natural sciences at the University of Cambridge (UK) and received his doctorate from the E.N.S.T. in Paris. He is now a professor of optics/photonics at the French "Grande Ecole" IMT Atlantique where he leads the diffractive optics group. His research interests include non-display applications of spatial light modulators, the design and fabrication of diffractive micro-optical elements and their applications in optical telecommunications and optical information processing.

**Andreas Erdmann** is the head of the Fraunhofer IISB Computational Lithography and Optics Group and teaches as a “Privatdozent” at the University of Erlangen. He has more than 25 years of experience in optical and EUV lithography. He chaired several SPIE conferences and is an organizer of the International Fraunhofer Lithography Simulation Workshop. He contributed to the development of several advanced lithography simulators including the development and research lithography simulator Dr. LiTHO. He is a fellow of SPIE.