© The Author(s), published by EDP Sciences, 2026

1 Introduction

Freeform micro-optical components are central to the next generation of compact optical systems, enabling precise control over light propagation in ways that conventional optics cannot achieve. To this end, GJEs based inverse designs represent an interesting alternative to the realization of such freeform based micro-optical systems. However, such a design strategy has not been explored in combination with 2PP based micro-fabrication technologies. By combining the best from both fields, we highlight the potential for creating novel micro-optical elements tailored for specific optical transformations. In the following sections, we provide details behind the employed numerical implementation, the followed micro-fabrication process and the methodology used to extract relevant information on the smoothness associated to the computed freeform surface solutions.

2 Numerical implementation

Our implementation relies on the work presented in [1] for which an iterative algorithm is used to obtain freefrom optical surfaces via the use of a GJE. Such a design algorithm exploits the relationship between a generating function G and its unique inverse H, also referred to as a Hamilton’s characteristic function. As in [1] we rely on an angular characteristic formulation to our problem, for which the source (x) and target (m) coordinates are both expressed in terms of stereographic values. For a system consisting of a point source with N interfaces prior to the sought freeform surface, the associated Hamilton’s characteristic function can be expressed as$$\begin{array}{ll}H(x,\mathbf{m},u_{\mathrm{ff}})& =u_1\left(x\right)\frac{2|x-\mathbf{m}{|}^2}{(1+|x{|}^2\left)\right(1+\left|\mathbf{m}{|}^2\right)}+\sum _{i=2}^{N-1}\mathrm{}\left[u_i\left(x\right)\left(n_i-1+\frac{2\left|y_{i-1}\right(x)-\mathbf{m}{|}^2}{(1+|y_{i-1}\left(x\right){|}^2\left)\right(1+\left|\mathbf{m}{|}^2\right)}\right)\right]\\ & +u_{\mathrm{ff}}\left(n_N-1+\frac{2\left|y_N\right(x)-\mathbf{m}{|}^2}{(1+|y_N\left(x\right){|}^2\left)\right(1+\left|\mathbf{m}{|}^2\right)}\right)\end{array}$$(1)

with u₁, u_i being the Euclidean distance functions between the N interfaces located before the freeform surface $u_{\mathrm{ff}}$, $\left\{y_{i-1}\right(x){\}}_{i=2}^{N-1}$ represent the set of stereographic coordinates associated to the outgoing ray-directions (I) at interfaces i = 2,…, N−1 while y_N (x) corresponds to the stereographic coordinates acting as local source ones to the last term in equation (1). Finally, u_ff represents the Euclidean distance between the Nth interface and the freeform surface. To implement the algorithm, we make use of an open source finite element method (FEM) solver [2] which represents two main advantages. First, it allows the use of an unstructured grid to solve the required numerical problems enabling higher flexibility in terms of the source domain shape. Second, it allows the use of local basis functions to describe the solutions in terms of polynomials within each element. Complementary to this, we employ a hybrid numerical scheme to obtain all gradient information required as part of the iterative routine. In this case, for all quantities in equation (1) excepting m and u_ff , we rely on automatic differentiation tools via an open source implementation of a differentiable ray-tracing framework [3]. On the other hand, for m and u_ff we make use of the polynomial preserving recovery(PPR) technique [4] which allows obtaining least squares approximations to the gradients at the domain’s degrees of freedom. Importantly, by making use of this approach, we obtain a sparse matrix operator which can be directly applied on each iteration. Once a set of solution functions m and u_ff are obtained, these are used to reconstruct a three-dimensional freeform surface using$$S=u_1\left(x\right)I_0\left(x\right)+\sum _{i=2}^{N-1}\mathrm{}\left[u_i\left(x\right)I_{i-1}\left(x\right)\right]+u_{\mathrm{ff}}{I}_N\left(x\right).$$(2)

Finally, as a means of comparison and to extract information on the smoothness associated with the obtained solutions, we have also chosen to use Forbes Q freeform polynomials [5–6], with these being commonly used to express freeform surfaces in terms of global orthonormal basis functions.

3 Fabrication details

The chosen freeform lens was fabricated by means of a commercially available 2PP system (Photonic Professional GT2 by Nanoscribe GmbH, Karlsruhe, Germany) in combination with a high numerical aperture objective (Plan-Apochromat 25×/1.40 Oil DIC, Zeiss). This objective was chosen given the obtained freeform lens dimensions, which has an elliptical shape with a major axis half-diameter equal to 400 μm. At the same time, a slicing and hatching distance of 100 nm were chosen, in order to have close fidelity to the original surface design and minimize discretization artifacts. The lens was fabricated directly on top of an AlGaInP visible laser diode (LD) with a central emission wavelength at 650 nm. These LDs consist of an emission surface protected by a window glass with a thickness of 400 μm. In order to align the freeform lens center to the emission point, an additional alignment step was performed. For this, a 10×/0.35 NA objective was used to find the lateral location of the LD’s emission center point, which is found 700 μm below the LD’s upper window interface. Consequently, two reference markers were fabricated, these being used in the subsequent freeform lens fabrication step to align the structure. In Figure 1, a schematic representation of the geometrical characteristics of an LD is provided, in addition to representative rays traced from the emission point location towards the freeform optical surface direction and indications of all involved terms directly contributing to equation (1).

Figure 1 Schematic representation for the laser diode system considered in this case. The emission region is modeled as a point emitter, positioned at the origin. 300 μm above this point, a window glass with a thickness of 400 μm is found. All of the interfaces u(x) and the associated intermediate ray-directions I(x) are directly included in the Hamilton’s characteristic function H shown in equation (1).

4 Results

In this section we present two freeform surface design examples obtained following the methodology introduced in Section 2. For this, two different characteristic target functions have been chosen. The first one consists of an uniform target irradiance distribution defined over a circular aperture, with a diameter of 5 cm, at a target distance of 4.5 cm from the origin. The second scenario consists of a modulated cosine-squared function, given by $\cos (\frac{2\pi x_t}{\beta }{)}^2$, where x_t denotes the x cartesian coordinates at the target plane and β the desired fringe periodicity. Both distribution are defined over the same spatial extent. As in [1], we make use of a standard transformation to obtain intensity functions from both irradiance distributions. For the source, we use the model$$I({\theta }_x,{\theta }_y)=I_0{\exp }^{-2\left(\right({\theta }_x/{\alpha }_x{)}^2+({\theta }_y/{\alpha }_y{)}^2)}$$(3)

with θ_x and θ_y being perpendicular emission angles associated with a Gaussian-like far-field intensity profile, while α_x and α_y represent half angles of divergence along both x and y axes. These were computed from the full-width at half-maximum angles provided in the manufacturer’s datasheet (28^∘ and 9^∘, respectively) using the relationship $\alpha =\frac{{\theta }_{\mathrm{fw}hm}}{\sqrt{2\ln 2}}$. Following, in Figure 2 we compare the Euclidean distance function solutions u_ff obtained for both target functions. Additionally, we contrast the extracted convergence curves in both cases, using the same definition for J_I as given in [1].

Figure 2 (a) and (b) provide the obtained uff function values after the iterative evaluations as described in Section 2. (c) Pointwise function difference Δuff between maps shown in (a) and (b) which encodes the local difference between both surface solutions. (d) provides convergence curves using the quantity JI , according to the definition provided in [1].

From a visual comparison between the functions shown in (a) and (b), a clear similarity can be observed. Interestingly, the overall global function shapes and values extend over the same range. To better contrast both u_ff functions, on (c) we display the pointwise function difference Δu_ff which reveals the local distinctions between both functions, with Δu_ff not exceeding values beyond 1 μm. Consequently, we performed non-sequential ray-tracing analysis with the reconstructed three-dimensional freeform lens models via equation (2). These results are presented in Figure 3 with, the original target distributions, the obtained mapping functions m and the irradiance information extracted from the non-sequential analysis, being shown. As described in Section 2, our numerical implementation is based on standard local polynomial basis functions. This results in C₀-only continuous solutions, with potential detrimental effects induced in the reconstructed three dimensional surfaces. To further evaluate this, we make use of an alternative surface representation based on Forbes Q freeform polynomials as briefly introduced in Section 2. We employ least-squares approximations defined directly on the original surface domains, by projecting the solutions into basis functions characterized by indices n and m, with all terms satisfying 2m + n ≤ T being included in the expansion. From the obtained projections, we calculate residue function maps ΔS_z by subtracting the polynomial fits from the FEM solutions. These residue maps, which are shown in Figures 4a and 4b were obtained with a T value of 54, resulting in a total of 1592 basis terms. From the chosen T, we estimate the minimum radial feature period being captured by the polynomial fit, by considering the shortest axis radius, yielding an approximate value of $\Delta r\approx \frac{200\mu \mathrm{m}}{T/2}=7.40\mu \mathrm{m}$. At this minimum radial period, a clear distinction in terms of the residue surface height can be observed.

Figure 3 (a), (b) and (c) show the target irradiance, obtained mapping functions m and extracted non-sequential irradiance distribution, respectively, associated with the freeform surface generated for the cosine-squared case (uff shown in Fig. 2a). Similarly, in (d), (e) and (f) the results associated with the uniform target case are presented.

Figure 4 (a) and (b) present fit residue maps ΔSz obtained after using least squares projections based on Forbes Q freeform polynomials to represent both FEM solutions. (c) presents power spectral density (PSD) cross-sections extracted from both residue maps shown in (a) and (b). (d) provides a comparison between ΔSz RMS curves as a function of the expansion index T.

For the uniform target case, a peak residue of 67 nm was obtained, with the large majority of the surface height structure being captured by the smooth polynomial fit. Differently, for the modulated-cosine target, the obtained fit residue contains relevant spatial features with height contributions in the order of 160 nm. In other words, the unresolved surface structure consist of spatial features bounded by an extension of approximately Δr with height amplitude contributions not larger than a value of λ/4.

Additionally, power spectral density (PSD) analysis was performed from the extracted ΔS_z residue maps. These results are presented in Figure 4c as a comparison between cross-sectional views of the obtained PSDs. Both curves reveal that spatial features located at frequencies larger than 1 μm⁻¹ are highly suppressed. Complementary to this, in Figure 4d we compare root-mean-square (RMS) estimates for ΔS_z as functions of the expansion index T. In this case, clear RMS reduction trends can be observed, with values of 3.38 and 7.75 nm achieved at T = 106, for the uniform and cosine-squared targets, respectively. Finally, the freeform lens corresponding to the cosine-squared target pattern was fabricated following the methodology described in Section 3. The obtained freeform lens was characterized using a white-light interferometer(WLI, Zygo, Nexview NX2, 10X, Δxy: 750 nm) for which the recorded intensity map is presented in Figure 5a.

Figure 5 (a) Intensity record by employed white-light interferometer system. (b) Surface residue fit ΔSz from extracted topography via WLI system. Only features with a magnitude bounded by λ are presented. (c) Recorded irradiance distribution from freeform lens fabricated directly on top of LD system in accordance to the process described in Section 3.

As performed for the design surface case, we make use of Forbes Q freeform polynomials (T = 54) to fit the measured surface height profile. In (b), a portion of the obtained height fit residue map ΔS_z with contributions bounded by a value of λ is presented. From the full extracted ΔS_z, a total residue RMS of 110 nm, in addition to a peak-to-valley value of 4.67 μm were estimated, which compare to values of 37 and 263 nm for the designed surface case at the same T value, respectively. These parameter differences can be attributed to the utilized fabrication parameters and development process, with material shrinkage representing a significant portion of these contribution. Following, we measured the irradiance distribution generated by the manufactured freeform lens. This was recorded with a CMOS camera in close proximity to the LD sample, for which the obtained irradiance map is presented in Figure 5c. From the recorded map, a distorted fringe pattern in comparison to the target one can be recognized, with higher intensity values obtained towards the domain’s boundary. Such a distortion is caused by the employed α_x and α_y values which were directly taken from the manufacturer datasheet and were not extracted from intensity measurements as part of the followed process. Nevertheless, despite all mentioned factors, the pattern presented in (c) still preserves most of the fringe-like structure from the intended target distribution. From the Euclidean distance function difference Δu_ff shown in Figure 2c with an RMS value of 345 nm and the results presented above in Figure 5, it can be concluded that the generated freeform surface design is well resolved by the employed 2PP fabrication process. Based on these results, we believe that the presented design and fabrication pipeline represents potential advantages for the realization of compact micro-optical freeform lens systems, benefiting from the flexibility offered by both, the employed GJE based design technique and the two photon polymerization based fabrication process.

5 Conclusions

In this work we have demonstrated the feasibility of combining novel inverse design algorithms with the micro fabrication capabilities offered by two photon polymerization. The presented integrated engineering pipeline benefits from the developments in both fields and serves as an important bridge between advanced design techniques and cutting edge micro fabrication technologies, enabling in this form, the further development of complex freeform surface based micro-optical systems.

Funding

This work was funded by the German Research Foundation (DFG) – project number: 537519988.

Conflicts of interest

The authors declare no conflicts of interest.

Data availability statement

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

Author contribution statement

All authors contributed equally to this work.

References

All Figures

Figure 1 Schematic representation for the laser diode system considered in this case. The emission region is modeled as a point emitter, positioned at the origin. 300 μm above this point, a window glass with a thickness of 400 μm is found. All of the interfaces u(x) and the associated intermediate ray-directions I(x) are directly included in the Hamilton’s characteristic function H shown in equation (1).

In the text

	Figure 2 (a) and (b) provide the obtained uff function values after the iterative evaluations as described in Section 2. (c) Pointwise function difference Δuff between maps shown in (a) and (b) which encodes the local difference between both surface solutions. (d) provides convergence curves using the quantity JI , according to the definition provided in [1].
In the text

	Figure 3 (a), (b) and (c) show the target irradiance, obtained mapping functions m and extracted non-sequential irradiance distribution, respectively, associated with the freeform surface generated for the cosine-squared case (uff shown in Fig. 2a). Similarly, in (d), (e) and (f) the results associated with the uniform target case are presented.
In the text

	Figure 4 (a) and (b) present fit residue maps ΔSz obtained after using least squares projections based on Forbes Q freeform polynomials to represent both FEM solutions. (c) presents power spectral density (PSD) cross-sections extracted from both residue maps shown in (a) and (b). (d) provides a comparison between ΔSz RMS curves as a function of the expansion index T.
In the text

	Figure 5 (a) Intensity record by employed white-light interferometer system. (b) Surface residue fit ΔSz from extracted topography via WLI system. Only features with a magnitude bounded by λ are presented. (c) Recorded irradiance distribution from freeform lens fabricated directly on top of LD system in accordance to the process described in Section 3.
In the text