1. Introduction

Two-dimensional (2D) image sensors are the most common detectors for light, so a leading optical engineering task is how to best extract the information from the incident optical field using an optical system and a planar image sensor. For example, a black and white camera maximizes the amount of spatial information by performing a one-to-one mapping between the direction of propagation and a pixel on the image sensor, but spectral and polarization information is lost. A line scan multispectral camera maps only one spatial coordinate to a direction on the image sensor while the other direction records the spectrum of the spatial pixel. Various other mappings are used in cameras to image color, polarization, light field, etc. Since the information capacity of a 2D image sensor is finite, getting more information about some degrees of freedom for light comes at the expense of information in other degrees of freedom. Also, in general purpose systems with trivial mapping implementations that use conventional optical components like lenses, gratings, and prisms, a single pixel on the image sensor detects a specific combination of single degrees of freedom. For example, a camera may detect the combination of a specific direction, wavelength band, and polarization. Thus, if $S$ spatial directions, $W$ wavelength bands, and $P$ polarizations need to be resolved, then $S \times W \times P$ pixels are needed, where $P$ is at most 4 to fully classify polarization through Stokes parameters [1,2].

Oftentimes there is prior knowledge about the input light field, in which case it is possible to use mappings that more efficiently utilize the pixels of the sensor [3–7]. Here, we demonstrate an inverse designed volumetric meta-optic device that can efficiently map different combinations of wavelengths, directions and polarizations into combinations of pixels on an image sensor. The compressed information can fully classify properties of the incident fields under certain approximations, namely that the wavefronts are monochromatic, locally linear in phase, and linearly polarized.

The device is based on metaoptics, which describes materials patterned with subwavelength resolution that impart customized transformations to incident light. Most research in dielectric meta-optics focuses on metasurfaces that consist of a single, approximately wavelength-thick layer of subwavelength scale antennas with the ability to completely control the phase and polarization of an incident wavefront [8–10]. The ability to impose independent phase profiles to two orthogonally polarized inputs allows a single metasurface to perform tasks that historically could only be achieved with cascaded systems of bulk components, which has improved applications that must obey strict size constraints while maintaining high efficiencies. However, further expanding the multifunctionality of metasurfaces to multiple angles and wavelengths tends to come at the cost of efficiency [11,12].

To recover the performance lost by increasing multifunctionality, the size of the metaoptics system can be increased. This idea has been explored in recent years by employing multiple metasurfaces in a cascaded system [13–18]. While these systems have a reduced size relative to systems that employ only conventional bulk optical elements, the size of the system is primarily dictated by the distance between metasurfaces. The distance between metasurface elements must be sufficiently large to preserve the assumptions that metasurface design is reliant on. Broadly speaking, placing meta-atoms directly in contact with other meta-atoms will alter their transmission in ways that are difficult to predict using only the individual response of each meta-atom. These considerations are described in detail in Refs. [19–22]. The paradigm of cascading elements that are spatially separated enough to preserve their design independence provides an intuitive way to design complex systems, but is not strictly necessary. Instead, design methodologies that account for the effects of near-field coupling and multiple scattering events (i.e. full-wave Maxwell solvers, including FDTD) are required for more advanced operation.

Recent work has demonstrated the tractability of designing 3D volumetric metaoptics using inverse-design techniques. These techniques are aided by the adjoint method for electromagnetics, which utilizes adjoint symmetry in Maxwell’s equations to efficiently compute the gradient of arbitrary figures of merit (FoMs) with respect to material permittivity [23–25]. This can be subsequently used to optimize the shape of a structure in a process referred to as topology optimization [26,27]. Full-wave simulations require significant computational resources to perform, so applications of inverse-designed metaoptics have been limited to integrated waveguide components [28] or small 3D components [29], on the order of several wavelengths per side. A recent pioneering work by Lin et al. co-optimized the scattering behavior of a multi-layer structure consisting of variable height polymer with a computational imaging reconstruction back-end to extract spatial, spectral, and polarization properties [30]. The device is placed over a $10\times 10$ grid of imaging sensors and offers remarkable noise insensitivity. The work discussed in this paper is similar, but features an increased emphasis on layered fabrication and a read-out strategy using straight-forward intensity ratios rather than a computational back-end to classify optical states.

For this work we simulate a device that can classify incident light based on the fundamental properties of propagation direction, wavelength, and polarization over a $3\times 3$ grid of imaging sensors. To our knowledge this amount of multi-functionality at both high efficiency and compactness has not been obtained to this degree in the past. The intent of the study done here is to explore the possibilities ahead for volumetric metaoptics, provide a meaningful goal for future photonic foundry processes, and pose a device that could be of high interest to the computational imaging community.

Due to the computational difficulty of this task, we target a platform in which a periodic array of the designed 3D devices can be placed above a sensor array. This allows the devices to be small during the design, then tiled to cover a desired area. The devices are designed for three functionalities: polarization splitting of two linearly polarized states, wavelength splitting of two distinct wavelengths, and a basic form of imaging in which five planewave components are focused to five different locations. This totals to nine states we seek to classify, so in practice the device can be placed above a 3$\times$3 grid of pixels such as CMOS or CCD arrays. Although the device is designed for a discrete set of states, it can be used to interpolate states on a continuum by analyzing the ratio of intensities among all the pixels.

The polarization and wavelength splitting functionalities are designed to work for all of the assumed incident angles, which include a normally incident case and non-normal incidences that are 5 degrees tilted from normal (polar angle $\theta =5^{\circ }$) with azimuthal angles $\phi \in \{0^{\circ },90^{\circ },180^{\circ },270^{\circ }\}$. The two design wavelengths are 532 nm and 620 nm. The two polarizations are linear and orthogonal, with one polarized in the xz-plane and the other in the yz-plane. The device is a 3$\mathrm {\mu}$m $\times$ 3$\mathrm {\mu}$m $\times$ 4$\mathrm {\mu}$m stack of 20 layers. Each 200 nm layer is comprised of titanium dioxide (TiO$_2$, n=2.4) and silicon dioxide (SiO$_2$, n=1.5), with a minimum feature size of 50 nm. The background material is assumed to be SiO$_2$. A schematic of the device, sensor array, and examples of the fields at the focal plane under a 620 nm xz-polarized planewave excitation angled at $(\theta,\phi ) = (5^\circ,90^\circ )$ is shown in Fig. 1. Note that since this $5^\circ$ polar angle is in SiO$_2$ with $n=1.5$, these planewave components couple to $7.5^\circ$ planewave components in air.

 

Fig. 1. The layout of the system, which features a device region placed above an array of sensors. (a) An incident planewave is input to a device comprised of SiO$_2$ and TiO$_2$ at an angle $(\theta, \phi )$. The bottom of the device is 1.5 $\mu m$ above a sensor array. The background material is assumed to be SiO$_2$. (b) The distribution of functionalities across the pixel array. The device focuses light to different pixels depending on the state of the input light. (c) The output of 620 nm, xz-polarized planewave input incident at an angle $(\theta,\phi )=(5^{\circ },90^{\circ })$.

Download Full Size | PDF

We begin by describing the inverse design process used to optimize the devices. Next, the performance of the device under the assumed input conditions is studied. Following this, we observe that the behaviors of functionalities are preserved and predictable when excited at states in between the states the device was optimized for. For example, as the input wavelength is continuously shifted from 532 nm to 620 nm, the ratio of the 620 nm pixel transmission to the 532 nm pixel transmission monotonically increases, while the imaging and polarization functions remain efficient. This occurs despite the device only being optimized in a narrow wavelength range around 532 nm and 620 nm. Similar behavior is observed in the imaging functionality when excited at input angles that were not explicitly optimized within a $5^{\circ }$ cone. The read-out from these pixels can thus be used to infer the state of the light incident on the device so long as the input is assumed to be a monochromatic planewave. We find that this behavior does not depend on the chosen mapping of functionality to specific pixels in the focal plane by analyzing the behavior of devices optimized for 20 unique pixel distributions.

2. Methods

The goal of photonic topology optimization is to find a refractive index distribution that maximizes an electromagnetic figure-of-merit (FoM). Since the device presented here is highly multi-functional, its optimization is multi-objective. Each objective is a mapping of an input to a FoM: the input is a planewave with a specific angle, polarization, and wavelength; the FoM is the power transmission through the desired pixel. The general procedure for this optimization is shown in Fig. S1 of Supplement 1. It consists of three main steps: first, the FoMs and their associated gradients are computed [24]; second, the gradients are combined with a weighted average [31]; third, the device is updated in accordance with the averaged gradient using either a density-based optimization of a continuous permittivity [32] or a level-set optimization of a discrete permittivity [33].

The individual FoM gradients are evaluated at every point in the design region using the adjoint method, which entails combining the electric fields in the design region for a “forward” and an “adjoint” simulation to compute the desired gradient. In this case, the forward case simulates the device under the assumed planewave excitation, and the adjoint case simulates a dipole (with a particular phase and amplitude based on the forward simulation) placed at the center of the desired pixel. This choice of sources optimizes the device to focus light to the location of the dipole. However, we record the performance of the device as power transmission through the desired pixel rather than intensity at a point, since power transmission better represents the signal a sensor pixel would record.

All unique forward and adjoint simulations are first simulated in Lumerical FDTD in parallel on a high-performance cluster (Caltech Resnick High Performance Computing Cluster). Next, the associated FoMs are recorded from the forward simulation results, and the FoM gradients are computed by combining the results of appropriate forward and adjoint source pairs. The individual gradients are spatially averaged in the z-direction for each layer, yielding a 2D gradient for each layer of the device. All FoM gradients are combined using a weighted average into a single gradient, still evaluated at every point in the design region. Information on this weighting procedure is described in Ref. [31]. To prevent the structure from being highly resonant, frequencies within a $\pm 5$ THz range (approximately $\pm 4$ nm) of the design frequencies are also optimized. This adds very little computational complexity since the FDTD method simulates a broadband pulsed excitation that already contains these frequencies. Finally, this interpolated gradient is used to update the permittivity of the device structure.

The optimization is done in two phases: a density-based phase, and a level-set phase. Each phase has a unique update procedure. In the density-based optimization the permittivity of the device is modelled as a grid of grayscale permittivity values between the permittivity of the two material boundaries. This permittivity representation is effectively fictitious (unless an effective index material can be reliably fabricated), and the goal is to converge to a binary device that performs well and is comprised of only two materials. We use the methods described in Ref. [32] to achieve this.

While the density-based optimization can converge to a fully binary solution, it is faster to terminate the optimization early and force each device voxel to its nearest material boundary. This thresholding step reduces the device performance, which we recover by further optimizing the device with a level-set optimization. Level-set optimization models the device boundaries as the zero-level contour of a level-set function ($\phi (x,y) = 0$), and thus benefits from describing inherently binary structures [34]. Empirically the final device performance is dependent on initial seed, hence the need for the improved density-based optimization that converges to a near-binary solution. Here we use level-set techniques to simultaneously optimize device performance and ensure the final device obeys fabrication constraints. This technique is described in detail in Ref. [33], and elaborated on in Supplement 1 Section S2. The full convergence plot for the optimization and the resulting device layers are discussed in Section S3 and S4, respectively, of Supplement 1. Details on computer hardware and simulation times are discussed in Supplement 1 Section S5.

3. Results

We quantify the performance of the device in two ways. First, we check the performance of the device under the excitation beams that were assumed when optimizing the device, which we refer to as training modes. Second, we study the device at different input angles, wavelengths, and polarization states to analyze its ability to classify states. We refer to these as validation modes.

The design methodology co-optimizes 20 different training modes featuring unique combinations of wavelength, polarization, and angle of incidence. The transmission to each pixel for each input state is shown in Fig. 2. For each state the transmission to the three correct pixels, the transmission to the six incorrect pixels, and the transmission elsewhere (oblique scattering, back-scattering, etc.) is shown. On average, the transmission to the correct pixels totals 47.7%, the transmission to incorrect pixels totals 16.8%, and the transmission elsewhere totals 35.5%. The imaging pixels tend to be more susceptible to error in the sense that more power is erroneously transmitted to them than the wavelength- or polarization-sorting pixels on average. A possible reason for this is that the different angled input states have higher correlation than the different spectral and polarization input states due to the small device aperture, giving rise to more cross-talk between imaging pixels. Undesired scattering could be explicitly minimized as part of the optimization in future tests. The most straight-forward way to do this would be to incorporate the minimization of transmission to these pixels as FoMs in the overall optimization.

 

Fig. 2. The performance of the device under excitation of training modes. The input angle, polarization state, and wavelength are represented by the y-axis. The x-axis represents the fraction of power scattered to different regions. The bold colored bars quantify transmission to the correct pixel, while the pale colored bars quantify transmission to the incorrect pixel. The grey bar quantifies the power that is back-scattered or obliquely scattered, thus not reaching any of the pixels in the focal plane.

Download Full Size | PDF

The design process optimizes the performance of the device only for the training modes, and it is unclear what the output of the device will be if the excitation parameters are continuously varied. Here we quantify the performance of the device at input states between that of the designed input states. We study the performance at wavelengths in between 532 nm and 620 nm, at a large set of angles within a $10^\circ$ input angle cone, and at arbitrary polarization states.

The full angular response of each individual pixel in a $10^{\circ }$ cone is shown in Fig. 3(a)-(i). The results along the horizontal dashed lines ($\phi =0^\circ$) and vertical dashed lines ($\phi =90^\circ$) are plotted in Fig. 3(j),(k). To obtain these plots, the transmission values are averaged across wavelength, polarization, or both depending on the functionality. For the wavelength sorting functionality the values are averaged across polarization; for the polarization sorting functionality the values are averaged across wavelength; and for the imaging functionality the values are averaged across both wavelength and polarization. The wavelength and polarization demultiplexing functionalities show a more uniform response across $(\theta, \phi )$ values than the imaging functionality, which means these functions are not highly dependent on incident angle. The imaging functionality is, by design, sensitive to the incident angle. The variation in azimuth for the polarization and spectral pixels is expected because the device is not symmetric, which is required since the desired output fields are not symmetric. However, these deviations are non-ideal, as deviations from a uniform angular response can cause uncertainties in angle classification to cascade to uncertainties in spectral and polarization properties.

 

Fig. 3. The angular sensitivity of the device. (a)-(i) plot the transmission to the different pixels in the focal plane for $\theta$ from $0^\circ$ to $10 ^\circ$ and $\phi$ from $0^\circ$ to $360^\circ$. The top-down left-right ordering of these plots are matched to the ordering of the pixels depicted in Fig. 1(b). The wavelength demultiplexing plots (a) and (g) are averaged in polarization. The polarization sorting plots (c) and (i) are averaged in wavelength. The imaging plots (b), (d)-(f), and (h) are averaged in both polarization and wavelength. (j),(k) are horizontal (j) and vertical (k) traces of the surface plots (a)-(i) along the black dashed lines. For (j) and (k), the line colors represent transmission to different pixels: the green line is the 532 nm pixel, the red line is the 620 nm pixel, the solid blue is the xz-polarization pixel, the dashed blue is the yz-polarization pixel, the solid yellow lines are the on-axis imaging pixels, and the dashed yellow lines are the off-axis imaging pixels. For (l), the green and red lines represent the transmission to the 532 and 620 nm pixels respectively. The solid blue line is the transmission to the correct polarization pixel given the polarization of the input source, and the dashed blue line is the transmission to the incorrect polarization pixel. For (l), all traces are averaged over all evaluated input angles within a $5^\circ$ cone.

Download Full Size | PDF

The red and green traces of Fig. 3(l) show that as the wavelength shifts from 532 nm to 620 nm, the transmission through the 532 nm pixel smoothly decreases, and the transmission through the 620 nm pixel smoothly increases. These transmission values are averaged over both polarization and across all simulated incident angles $\phi \in [0, 360^\circ ]$ and $\theta \in [0,5^\circ ]$. As expected, the transmission to the 532 nm pixel is maximized for 532 nm wavelengths, and likewise for the 620 nm pixel. The traces smoothly vary between these two wavelengths, indicating that the power is predictably redistributed between the pixels in these cases. The solid blue line represents the average transmission to the correct polarization pixel (e.g., xz-polarized input being focused to the xz-polarization pixel), and the dashed blue line represents the average transmission to the incorrect polarization pixel (e.g., xz-polarized input being focused to the yz-polarization pixel). These transmission values are averaged over wavelength, and are averaged over incident angle in the same way the red and green traces are. While there is a drop in efficiency at the non-optimized wavelengths, there remains a high contrast between the solid and dashed blue lines at all wavelengths.

The imaging functionality transmission values vary smoothly with respect to incident angles $(\theta,\phi )$, but this occurs in a less intuitive way than it would with a lens. Rather than the focus spot continuously shifting across the focal plane as $(\theta,\phi )$ varies, the focused spots instead brighten or dim while remaining approximately centered in the individual pixels. As an example, when $(\theta,\phi ) = (2.5^{\circ },45^{\circ })$ the light is primarily scattered to the center of the $(0^{\circ },0^{\circ })$, $(5^{\circ },0^{\circ })$, $(5^{\circ },90^{\circ })$ pixels, as well as the relevant wavelength and polarization pixels. In this specific case 532 nm light is not focused to the 620 nm pixel which is the pixel that the incident planewave is oriented towards, and as such is where light would be focused if the device were replaced with a lens. Instead, the wavelength sorting and polarization sorting functionalities are mostly preserved under all excitation angles within an acceptance cone. Visualization 1, Visualization 2, Visualization 3, and Visualization 4 show the output fields changing as these parameters are swept, with further details available in Supplement 1 Section S8.

The effect of altering the polarization state is predictable since the polarization state of any input excitation can be described by two orthogonal linearly polarized states with a relative amplitude and phase shift between the orthogonal components. As an example, a simple experiment involving a linearly polarized source input into a Wollaston prism would show that as the source is rotated, the output of the Wollaston prism fluidly shifts from one linearly polarized output beam to the other. This behavior is observed in this device in the polarizing functionality, where the ratio between the polarizer pixels transmission can be used to infer the relative power of the two orthogonal linear polarization states. However, we would hope that the device performance for the wavelength and imaging functionalities is not adversely effected by the polarization states that were not explicitly optimized for. The efficiency of these functionalities may be altered as the cross-polarized output of one beam interferes with the parallel-polarized output of the orthogonally polarized input beam. We define the amount of cross-polarization by integrating the ratio $|E_x|^2/|E_y|^2$ under yz-polarized illumination over the output plane for all incident angles. The highest amount of cross-polarization is −11.9 dB at 620 nm and −13.6 dB at 532 nm. To verify that the various functionalities are not strongly affected, we analyze the results of all functionalities at all possible input polarization states by sweeping the relative amplitude and phase of the orthogonal input components, and we observe that interference effects can alter the transmission to the various pixels by a few percent. Data and further commentary on this matter is available in Supplement 1 Section S6. Note that if a particular design seeks to further minimize cross-polarization, then cross-polarization can be explicitly minimized during the optimization. Doing so will consume some design degrees of freedom, which may detract from the efficiency of the various functionalities, but the optimization will not require more time since the number of electromagnetic simulations per iteration will be unchanged.

Based on these results the device could be used to analyze the state of incident light at wavelengths between 532 nm and 620 nm and inputs with incident angles within an approximately 5$^{\circ }$ cone. Additionally, the relative amplitude of the xz-polarized and yz-polarized component can be classified. The classification is done by measuring the ratio of intensity between pixels and using a look-up table to then approximate the state of light. In the case of classifying wavelength and imaging angle, the behavior of smoothly transitioning between states is non-trivial and is not exhibited in typical optical devices such as gratings and lenses. These components tend to spatially shift a beam in the focal plane as the wavelength is varied (grating) or the incident angle is varied (lens), whereas our device only alters the ratio of transmission to the various pixels with very little spatial shift of the beams. In the case of classifying polarization this behavior tends to occur naturally in other optical devices such as birefringent prisms because any polarization state can be decomposed into a coherent sum of two orthogonally polarized states.

While classifying the wavelength and input angle states without ambiguity is not significantly affected by the incident polarization state it does require that the incident light be a monochromatic planewave. The planewave assumption is a common assumption employed in Shack-Hartman sensors and plenoptic sensors, and the assumption of monochromaticity can be satisfied with color filters, or if the light fundamentally comes from a narrow-band source such as a laser. Thus, there are numerous applications in which this device can be useful. The device can be tiled across a sensor array to enhance the functionality of an image sensor [29,35,36]. Such an array could be used to accurately classify the properties of a laser beam, including all fundamental properties of wavelength, polarization, and incident angle within the device’s acceptance cone. The angle-dependent nature of the device is similar in principle to angle-sensitive CMOS pixels [37], and could be used for lightfield imaging since the incident wavefront angle can be computationally determined using the relative intensity of the imaging pixels [5,38–40]. If coupled with a device such as a tunable bandpass filter, then a wavelength-dependent light field image can be obtained, with the added functionality of measuring the relative intensity of the two orthogonal linear polarization components for basic polarimetry. In general, we believe this type of device will open up new degrees of freedom in computational imaging applications.

One intriguing question we investigated is whether the qualitative behavior of the device depends on the specific assignment of functionalities to pixels. There are a large number of pixel permutations, and currently it is not easy to design a device for all of them. Instead, we investigated twenty different pixel distributions. Pixel distributions were chosen randomly, but were chosen to be sufficiently different from the original device by satisfying two criteria: 1) all pixels had to be moved from their original location, 2) no pixel could be rotated by $90^\circ$, $180^\circ$, or $270^\circ$ degrees, thus ensuring that devices were not similar by any rotational symmetry. Each device was optimized until the device was more than 15% binary according to the following equation, where $B$ is the binarization, $N$ is the total number of permittivity points in the design region, $\epsilon _i$ is the permittivity at the $i$th point, and $\epsilon _{min}$ and $\epsilon _{max}$ are the smallest and largest values of permittivity allowed in the design.

 

Fig. 4. The dependence of device performance on pixel distribution. For each functionality the overlap integral in Eq. (3) is evaluated on all combinations of pixel distributions. The average overlap integral is greater than 96.9% for all functionalities, suggesting that the arbitrary choice of mapping functionalities to specific pixels has little effect on the behavior of the device.

Download Full Size | PDF

4. Discussion

We presented simulations of a multi-functional 20-layer device made of SiO$_2$ and TiO$_2$ for use as a compressed sensing platform for the fundamental properties of an incident light wave: polarization, wavelength, and incident angle. The ability to classify wavelength and angle is contingent upon assumptions of the incident light. Monochromaticity is required for uniquely classifying the wavelength, which is a common assumption made in laser applications. The polarization functionality only measures the amplitudes of x- and y-polarized components, so a linear polarization is required for unambiguous classification. Future work could investigate introducing two additional polarization-classifying pixels to enable full-Stokes polarimetry. The last requirement is a locally linear phase front for classifying the incident angle, which is a commonly employed assumption used with Shack-Hartman and plenoptics sensors for applications such as wavefront sensing and computational imaging. This assumption is valid when the device is laterally small. Microlens arrays are often used for these applications, but the physical size of a microlens tends to be much larger than the devices presented here in order to span a sufficiently large number of sensor pixels.

The device was designed using topology optimization, aided by adjoint-based inverse design. The multi-objective optimization used techniques to approach a nearly binary solution during the continuous permittivity design phase, then concluded with a level-set approach that imposed fabrication constraints. It was found that the device functions at many more input conditions than those that were explicitly optimized for, which is a useful generalizability that can ease the computational burden of these large-scale optimizations. It also speaks to the robustness of the device, since its behavior is broadband in both spectral and angular space thus avoiding highly resonant and angularly sensitive designs.

A desirable tool in inverse-design of volumetric structures is a good preliminary estimate of device thickness and number of layers needed prior to optimization. While some recent efforts have focused on this task [22], it remains difficult to predict. In this work, the hyper-parameters were chosen based on prior experience with inverse-design optimizations and based roughly on computationally faster 2D simulations of a similar structure. Ultimately, the structure discussed here represents only one set of hyperparameters in a broad design landscape that may contain better-performing devices.

The fabrication of this device would be ambitious, but possibly achievable with state-of-the-art current technology. The optimization procedure ensures that the device obeys a minimum feature size fabrication constraint of 50 nm, which is obtainable with optimized electron-beam lithography. A twenty-layer device with precise layer-to-layer alignment could be fabricated using similar processes used in CMOS foundries, although such processes are primarily used for electronics rather than optics and the materials used are not suitable for low-loss optical devices. A promising recent advancement involves designing a photonic device in a CMOS platform, then removing the metal layers with a wet etch to leave only materials made of the dielectric layers and air [41]. This subtractive photonics platform could potentially be used to fabricate the devices presented here, particularly if the air was backfilled with TiO$_2$. In the near future, scaled-up versions of this device could be fabricated at infrared or terahertz frequencies to enable new imaging technologies. Scaling up by approximately 4x to near-IR/mid-IR wavelengths will raise the minimum features size to 200 nm, which is achievable with commercially available two-photon polymerization 3D printing techniques [7,42].