1. Introduction

Advances in 3D-printing techniques like the two photon lithography (2PL) keep extending the range of application for microoptics in areas such as optical imaging and metrology for technical and biological applications [1–3], quantum communication [4,5], photonic circuits [6,7], and optical trapping [8] to name a few. 2PL enables highly flexible and rapid design and production of microoptics.

Until now, few reports like [9], where polarization is used as a degree of freedom in the optical design of 3D-printed microoptics, can be found. As a reason, we suspect the lack of fast and accurate wave optical simulation methods that can handle sufficiently large simulation volumes while considering polarization, which requires vector wave treatment. In most cases, rigorous methods are not useful due to excessive runtime and hardware demands, therefore simplified simulation routines have been proposed. Unfortunately, those typically exhibit downsides like limitation to the paraxial case [10], restrictions on the properties of the simulated optics, e.g. small element thickness, low numerical aperture (NA) [11], or trade fast runtime with even more excessive hardware requirements [12,13]. These factors impair the utility of those simulation routines for the practical applicability in optical design of microoptics, which usually extend several hundreds of microns per dimension.

To overcome these restrictions, this work aims to develop a simulation algorithm that carries out vector wave light propagation through large volume microoptics. Moreover, sufficiently fast runtime is desired in order to complement a rapid prototyping process. As basis for our fast vector wave optical simulation algorithm, we use the relatively slow vector wave propagation method (VWPM) [14]. The VWPM is a vectorial extension to the scalar WPM, which itself has routinely been used as a simulation tool for 3D-printed microoptics [2,8,15,16]. For simulation of realistic 3D-printed microoptics, which extend several hundreds of microns per dimension, the runtime of the VWPM-algorithm is not feasible in most cases. However, an approach for runtime optimization, which has been developed for the scalar WPM [17], can be adapted and applied to the VWPM-algorithm for optimization.

In this paper, we present the Fast Polarized Wave Propagation Method (FPWPM), which is a runtime optimized version of the VWPM-algorithm. Accuracy of the FPWPM is validated by comparison with rigorous methods in simulation examples. Furthermore, we demonstrate its capability of vector wave optical simulation in comparably large volumes by evaluation of realistic microoptical designs, and show the runtime improvement by several orders of magnitude compared to the VWPM-algorithm.

2. Results

2.1 VWPM-algorithm

First, we briefly summarize the basic principle of the VWPM-algorithm, which served as basis for this work, and introduce notations in order to facilitate reading of our later modifications to the algorithm. For a detailed description, the reader may be referred to [14].

The VWPM-algorithm carries out a step wise propagation of an incident vector electric field through a given refractive index distribution of a volume, constituted by a number of $N_x \times N_y \times N_z$ voxels with sizes $d_x \times d_y \times d_z$, each assigned a discrete complex refractive index value $\tilde {n}$. The $z$-axis is chosen as the main axis of propagation. Using the transverse coordinate vector $\mathbf {r}_\perp = (x,y)$ and the discrete position $z$ along the $z$-axis, the spatial dependence of the refractive index can be expressed as $\tilde {n}(\mathbf {r}_\perp,z)$.

The electric field propagation from a plane $z$ to the subsequent plane $z+dz$ is computed in Fourier space by multiplication of the incident electric field vector with the matrix $\mathbf {T}(\mathbf {r}_\perp,\mathbf {k}_\perp,z,z+dz)$. This matrix is obtained by solving the vector Helmholtz equation for an infinitely extended interface perpendicular to the $z$-axis between the media $\tilde {n}(\mathbf {r}_\perp,z)$ and $\tilde {n}(\mathbf {r}_\perp,z+dz)$. The explicit expression for $\mathbf {T}(\mathbf {r}_\perp,\mathbf {k}_\perp,z,z+dz)$ will be discussed in section 2.2. The spatial distribution $\mathbf {E}(\mathbf {r}_\perp,z+dz)$ of the propagated electric field can be obtained by local superposition of all plane wave components, i.e. summation over the transverse spatial frequency components $\mathbf {k}_\perp = (k_x, k_y)$ at each point $\mathbf {r}_\perp$ in the transversal plane:

2.2 Runtime optimization of the VWPM-algorithm

For numerical evaluation of Eq. (1), a summation over $N_x \times N_y$ spatial frequencies at $N_x \times N_y$ spatial coordinates is necessary. Hence, the numerical effort for a propagation step is proportional to $\mathscr {O}\left ( N_x^2 \times N_y^2\right )$. As a consequence of the quadratic scaling per dimension, the runtime escalates rapidly with increasing number of data points, i.e. with increasing simulation volume. To enable the fast simulation of comparably large volumes using the VWPM-algorithm, we address this scaling issue by adapting an approach which has originally been introduced for the scalar WPM [17]. In short, by eliminating the dependence on the transverse spatial coordinate $\mathbf {r}_\perp$ in $\mathbf {T}^{(2)}(\mathbf {r}_\perp,\mathbf {k}_\perp,z,z+dz)$, the integral in Eq. (1) will be expressed in a way that the inverse fast Fourier transform (IFFT) can be applied. Hence, the computational cost for each propagation step is reduced by several orders of magnitude.

Part of the spatial dependence in $\mathbf {T}^{(2)}(\mathbf {r}_\perp,\mathbf {k}_\perp,z,z+dz)$ can be avoided by slight simplification of the VWPM-algorithm: in case that two adjacent voxels in the transversal plane have the same refractive index, the spatial derivatives $\partial /\partial _x \tilde {n}^2(\mathbf {r}_\perp,z)$ or $\partial /\partial _y \tilde {n}^2(\mathbf {r}_\perp,z)$ in Eq. (2) vanish. For the case that the number of interfaces inside the transverse plane and the portion of energy transmitting through those interfaces are reasonably low, the contribution of the terms towards the computation of $E_z(\mathbf {r}_\perp,\mathbf {k}_\perp,z)$ in Eq. (2) is assumed to be small, and hence is chosen to be neglected. Instead of Eq. (2), we obtain

The remaining dependence on $\mathbf {r}_\perp$ in Eq. (5) and all included quantities ultimately results from $\tilde {n}(\mathbf {r}_\perp,z)$ and $\tilde {n}(\mathbf {r}_\perp,z+dz)$ being functions of $\mathbf {r}_\perp$, i.e. the refractive index distribution being inhomogeneous within the transversal plane. Schmidt et al. [17] proposed that an inhomogeneous refractive index distribution $\tilde {n}(\mathbf {r_\perp },z)$ can be decomposed into a discrete number $M$ of homogeneous subregions. Subsequently, $\tilde {n}(\mathbf {r_\perp },z)$ can be expressed as a sum of these homogeneous regions:

Applying the preceding procedure on $\mathbf {T}^{(2)}(\mathbf {r}_\perp,\mathbf {k}_\perp,z,z+dz)$ in Eq. (5) leads to

For $\tilde {n}_r = \tilde {n}_s$, $\mathbf {T}^{(2)}_{r,s}(\mathbf {k}_\perp )$ equals the identity matrix and Eq. (11) further simplifies to

Using fast Fourier transform (FFT)-algorithms, the computational effort for each propagation step is now proportional to $\mathscr {O}\Big ( M^2 \left [N_x \log (N_x) \times N_y \log (N_y)\right ]\Big )$. As long as $M$, i.e. the total number of different refractive indices inside the simulation volume, remains reasonably low, the computation effort is reduced significantly. This condition is easily met for simulation of 3D-printed microoptics, since the optics usually consist of one photoresist surrounded by a homogeneous medium, e.g. air or water. Hence, $M = 2$ can be assumed to be the usual case.

When necessary, $E_z(\mathbf {r}_\perp,z+dz)$ can be computed explicitly using

2.3 Evaluation of the validity of the FPWPM in simulation examples

Due to the substantial modifications to the VWPM-algorithm, we first evaluate the accuracy of the FPWPM. In the subsequent examples, we focus on challenges for the simulation of 3D-printed microoptics, which specifically demand vectorial simulation and hence exceed the scope of the scalar WPM. Therefore, results are compared to analytical expressions where possible, otherwise well-known or rigorous simulation routines are used for reference. Further, as the main result of this work, we investigate runtime advantage of the FPWPM compared to the VWPM algorithm and other references.

2.3.1 Light transmission through a plane dielectric interface with perpendicular orientation to the propagation axis

A simple test for considering of polarization during electric field propagation accurately is offered by studying the transmittance through a plane dielectric interface which is oriented perpendicularly to the propagation axis. For this case, we simulate the transmittance of incident light as a function of the incident angle $\theta$ at an interface between $n= {1.5}$ and vacuum and vice versa. For theoretical reference, the transmittance is calculated using Fresnel’s equations. As can be seen in Fig. 1, FPWPM and theory agree perfectly. This is to be expected for this case, since the FPWPM is derived to solve the vector Helmholtz equation rigorously, as long as only interfaces perpendicular to the propagation axis are considered.

 

Fig. 1. Comparison between FPWPM and Fresnel’s equations for calculation of light transmittance through a dielectric interface. The orientation of the interface is perpendicular to the propagation axis, the materials used for simulation are glass (g) with $n={1.5}$ and vacuum (v).

Download Full Size | PDF

2.3.2 Light incidence at plane dielectric interfaces with oblique orientation in respect to the propagation axis

In general, the interfaces between media can be shaped arbitrarily for microoptics. Especially for 2PL as fabrication method, one of the main advantages is the high flexibility regarding production of surface shapes. Other than the resolution limit of the 2PL device, barely any limitations are imposed. Hence, spheric, aspheric [18,19], diffractive [8,20], combination of e.g. refractive and diffractive surfaces [2] and even metasurfaces [21,22] can be fabricated and thus are of interest in simulation. For the FPWPM, the simulation volume has to be discretized along the Cartesian coordinate system, which naturally leads to staircase approximations for arbitrarily shaped interfaces. As a result, surfaces with parallel orientation in respect to the propagation axis occur. Like in the VWPM-algorithm, rigorous treatment of electromagnetic boundary conditions for light transmission through those interfaces is not ensured in the FPWPM, which may give rise to errors.

In the following, we investigate the magnitude of those errors. Therefore, we simulate the incidence of a Gaussian shaped laser beam with $\lambda = {0.55}\;\mathrm{\mu}\textrm{m}$ at plane interfaces which are not perpendicular to the propagation axis $z$. The polarization of the beam is linear and energy is distributed evenly in TE- and TM- polarization. For the media at the interface, we choose vacuum in front of the interface and glass with $n = {1.5}$ behind the interface.

First, we briefly explore the performance of the FPWPM-algorithm at its limits of validity by simulating total internal reflection (TIR) at an interface which is almost parallel to the propagation axis $z$. The propagation axis of the incident laser beam is parallel to the $z$-axis. As can be seen in Fig. 2(a), which shows ${\vert E_{\mathrm {TM}}\vert }$, the TIR of the laser beam can in principle be simulated. The propagation direction of the beam after the TIR at the interface is correctly predicted within the discretization limits of the simulation. The electric field amplitude in the reflected beam deviates by ${4.05} \%$ for TE- and ${14.57} \%$ for TM-polarization. Consequently, regarding electric field amplitudes, the FPWPM shows limited applicability in simulation of interfaces which are close to being parallel to the $z$-axis, especially for TM-polarization.

 

Fig. 2. Simulation of light incidence on interfaces with oblique orientation to the $z$-axis, using the FPWPM. (a) Illustration of the geometry for subsequent simulations. $\alpha$ represents the angle between the incident laser beam and the propagation axis $z$, $\theta$ represents the angle of incidence of the laser beam in respect to the surface normal. ${\vert E_{\mathrm {TM}}\vert }$ is plotted for $\alpha = {15}^{\circ}$ and $\theta = {40}^{\circ}$ in respect to the $z$-axis. (b) ${\vert E_{\mathrm {TM}}\vert }$ for $\theta = {75}^{\circ}$, showing that the FPWPM enables qualitative simulation of TIR and variable orientation of interfaces.

Download Full Size | PDF

Subsequently, we systematically analyze the accuracy of the FPWPM-algorithm for simulation of light transmission through interfaces with oblique orientation in respect to the propagation axis. Therefore, we carry out simulations with incident beam propagation angles of $\alpha = {0}^{\circ}$ and $\alpha = {15}^{\circ}$ in respect to the $z$-axis for an incident angle range of ${0}^{\circ} \leq \theta \leq {60}^{\circ}$ respectively. Both angles are illustrated in Fig. 2(b), while ${\vert E_{\mathrm {TM}}\vert }$ is plotted for visualization. Additionally, taking into account possible effects of simulation discretization, we run each simulation with both $dy = dz = \lambda /10$ and $dy = dz = \lambda /30$.

The resulting relative errors in simulated transmitted energy compared to analytical reference obtained using Fresnel’s equation are plotted in Fig. 3(a). For $\alpha ={0}^{\circ}$, the relative error at $\theta = {0}^{\circ}$ vanishes, since the FPWPM is rigorous in this particular case. With increasing angle of incidence $\theta$, the relative error tends to grow, but remains in the order of magnitude of a few percent, with its maximum of approximately ${5} \%$ at $\theta = {60}^{\circ}$ for TM-polarization and $dy = dz = \lambda /30$. For $\alpha = {15}^{\circ}$, the relative error tends to be larger than in the previous simulations. For TE-polarization, the relative transmitted energy tends to be slightly too high for low values of $\theta$, progressing to be too low while $\theta$ increases. Contrary behaviour can be observed for TM-polarization. While for TE-polarization, the absolute value of the error remains in the range of approximately ${5} \%$, the error for TM-polarization initially is in the order of a few percent, then deteriorates up to approximately ${15} \%$ for very large angles of $\theta$. Thus, the applicability of the FPWPM to quantitatively simulate energy transmission for TM-polarized light is limited for very large angles. Notably, finer discretization does not lead to lower relative error. Rather, finer discretization leads to higher transmitted energy in simulation. Possibly, this is due to the larger absolute number of propagation steps for finer discretization along the $z$-axis. Numerical gain due to the application of Fourier transforms to discontinuous electric fields, as reported by other authors [14,17], can hence grow more impactful. As a consequence, unnecessarily fine discretization should be avoided in simulation. To summarize, the errors are within a few percent for most of the parameter space under investigation. Given the benefits in computation speed and applicability of the FPWPM in practice, which will be demonstrated in the subsequent sections, the mostly minor errors are well acceptable, but should be noted nonetheless.

 

Fig. 3. (a) Relative error in transmitted energy for FPWPM-simulations compared to analytical reference from Fresnel’s equations, depending on the angle of incidence $\theta$, tilt angle $\alpha$ and the simulation discretization. (b) Deviation of the angle of refraction computed with FPWPM and analytical reference using Snell’s law. The dotted orange lines indicate the angular resultion limit of the FPWPM-simulation imposed by the discretization in Fourier space.

Download Full Size | PDF

Figure 3(b) shows the deviation between the angle of refraction simulated using the FPWPM and analytical reference obtained using Snell’s law. For the FPWPM, we computed the angle of refraction from the Fourier spectrum by selecting the peak energy frequency component and calculating the corresponding propagation angle of light. To outline the order of magnitude of the deviation, the orange dotted lines indicate the angular discretization in Fourier space of the simulation for the chosen parameters. To maintain clear visualization, the orange lines only depict the minimum value of angular discretization $\Delta \beta$, which strictly speaking only holds for $k_y = 0$ and its adjacent spatial frequencies and can be computed using $\Delta \beta = \arcsin \left [2\pi /\left (N_y\hspace {0.1mm}dy\hspace {0.1mm}n\hspace {0.1mm}k_0\right )\right ]$. Larger spatial frequencies correspond to increasingly coarse angular discretization. Clearly, for all conducted simulations, the deviation in computed angle of refraction is in the range of simulation discretization. Hence, the angle of refraction can be simulated precisely up to the Fourier space discretization limit of the simulation. Fixing all other simulation parameters, TE- and TM- polarization show identical angle of refraction as would be expected, with one exception for the simulation with $\alpha = {15}^{\circ}$, $\theta = {33}^{\circ}$ and $dy = dz = \lambda /30$. This exception might be caused by our method of determining the angle of refraction in terms of the frequency component with the largest energy content, which in principle is sensitive to amplitude errors. Furthermore, the evolution of the deviation in angle of refraction over the angle $\theta$ appears to be identical for both investigated simulation discretizations, apart from a constant offset. For all simulations with $dy = dz = \lambda /10$, $N_y = 2727$, which is odd, while $dy = dz = \lambda /30$ leads to $N_y = 8182$ which is even. The different symmetry properties lead to an offset of half a sampling point in Fourier space, which agrees well to the offset of calculated angles of refraction in Fig. 3(b). A more elaborate, yet robust method to avoid both those sources of inaccuracy for determining the angle of refraction might be found in calculating the angle of refraction in terms of the center of gravity of the Fourier spectrum.

2.3.3 Investigation of a high frequency blazed grating

Next, the performance of the FPWPM for simulation of diffractive mircooptics is analyzed. We therefore investigate diffraction at a high frequency grating with period $p = {2}\;\mathrm{\mu}\textrm{m}$ and compare the results to rigorous coupled wave analysis (RCWA). The refractive index of the grating material is $n={1.5}$. Each grating period is composed of three equidistant steps with linearly increasing thickness, approximating a sawtooth-shape. With $\lambda = {0.55}\;\mathrm{\mu}\textrm{m}$, the lateral feature size of the grating is close to the size of the simulation wavelength. The peak thickness of the grating is 1.1 µm. To ensure sufficient resolution in Fourier space when computing diffraction spectra in subsequent analyses, a total number of 20 grating periods is simulated with both FPWPM and RCWA.

For comparison, $\Re \{E_{y}\}$ obtained with the FPWPM and RCWA are depicted in Figs. 4(a) and (b), respectively. In contrast to the absolute value of the electric field, which has been used for investigation in prior examples, observation of the physical electric field $\Re \{E_{y}\}$ reveals additional information about the electric field oscillation. Both simulations clearly agree very well. Notable difference is limited to the area inside the grating material, most likely caused by omission of back reflected waves in the FPWPM, contrary to the fully rigorous RCWA simulation.

 

Fig. 4. Simulation of a high frequency grating with FPWPM and RCWA, with simulation wavelength $\lambda = {0.55}\;\mathrm{\mu}\textrm{m}$ and grating period 2 µm. (a) $\Re \{ E_{\mathrm {TE}}\}$, FPWPM. The dashed lines frame the grating material. (b) $\Re \{ E_{\mathrm {TE}}\}$, RCWA. (c) Diffraction efficiencies of FPWPM and RCWA, TE-polarization. (d) Diffraction efficiencies of FPWPM and RCWA, TM-polarization.

Download Full Size | PDF

For quantitative analysis, the diffraction efficiencies of TE- and TM-polarized fields at $z = {5}\;\mathrm{\mu}\textrm{m}$ are computed and plotted in Figs. 4(c) and (d) respectively. The locations of the diffraction orders in the frequency spectrum agree perfectly between FPWPM and RCWA, while slight deviations in amplitude can be observed, especially for TE-polarization.

2.3.4 Fast focal intensity simulation of a large, high NA microlens

Next, we demonstrate the feasibility of the FPWPM to carry out fast vector wave simulation of 3D-printed microoptics in large simulation volumes. Therefore, we design a large microlens which tightly focuses normally incident light and study the focal intensity using the FPWPM.

A microlens with diameter $d = {305}\;\mathrm{\mu}\textrm{m}$, thickness $t = {130}\;\mathrm{\mu}\textrm{m}$ and $\text {NA} = {0.72}$ is designed and optimized using the commercial software Optic Studio using ray optics optimization. For illumination, we use normally incident, circularly polarized plane waves with $\lambda = {0.635}\;\mathrm{\mu}\textrm{m}$. In order to obtain prominent polarization effects in the focal plane of the microlens, orbital angular momentum (OAM) is added to the plane wavefront via transmission through a vortex phase plate prior to focusing. An illustrated cross section through the simulated volume, constituted by the vortex phase plate and the microlens in a surrounding of vacuum, is depicted in Fig. 5(a). Depending on the handedness of the circular polarization of the incident light, the z-component of the electric field is expected to interfere either constructively or destructively on the optical axis in the focal region of the microlens.

 

Fig. 5. Simulation of electric field propagation through a microlens for normal plane wave incidence at $\lambda = {0.635}\;\mathrm{\mu}\textrm{m}$. (a) Illustrated YZ-cross section of the refractive index distribution used for simulation. Upper right corner: Profile of the vortex phase plate profile, XY-plane. (b) $\vert E_{\mathrm {TE}} \vert$ for 2D-simulation with TE-polarized light ($dy=dz = {0.05}\;\mathrm{\mu}\textrm{m}$). (c)–(h) Simulated intensity in the focal plane of the microlens. The handedness of the circular polarization used in the simulations is indicated in the upper right corner. Bars: 0.5 µm. (c) + (d) FPWPM, $dx=dy=dz = {0.127}\;\mathrm{\mu}\textrm{m}$. (e) + (f) FPWPM, $dx=dy=dz = {0.05}\;\mathrm{\mu}\textrm{m}$. (g) + (h) RW-Integral, $dx=dy= {0.05}\;\mathrm{\mu}\textrm{m}$. (i) Comparison of cross sections through foci in (f) and (h) and the focus obtained with the scalar WPM.

Download Full Size | PDF

For quick estimation and illustration of the focal position as the region of interest for subsequent focal intensity studies, a 2D-simulation is run using TE-polarized light, as depicted in 5(b). Next, full 3D-simulations with the FPWPM algorithm are carried out and normalized focal intensity distributions are extracted for left-handed-circular (LHC) and right-handed-circular (RHC) polarization, respectively. For first simulations, a large sampling step size of $dx=dy=dz = \lambda /5 = {0.127}\;\mathrm{\mu}\textrm{m}$ is chosen to minimize runtime, resulting in $N_x\times N_y\times N_z ={2400} \times {2400} \times {2330}$ discrete data points to constitute the whole 3D simulation volume of ${305}\;\mathrm{\mu}\textrm{m} \times {305}\;\mathrm{\mu}\textrm{m} \times {296}\;\mathrm{\mu}\textrm{m}$. The simulation resolution in Figs. 5(c) and (d) is rather coarse, nevertheless effects of illumination polarization on intensity at the center of the foci are clearly visible. Finer sampling step size of $dx=dy=dz = \lambda /5 = {0.05}\;\mathrm{\mu}\textrm{m}$ leads to ${6100} \times {6100} \times {5920}$ data points for the whole simulation volume. The resulting focal intensity distributions, which are plotted in Figs. 5(e) and (f), qualitatively agree well with the coarse simulation in Figs. 5(c) and (d).

Assuming an ideal aplanatic lens, the well-known Richards-Wolf (RW)-integral [23,24] can be used as reference to compute the focal intensity. Results are depicted in Figs. 5(g) and (h). The RW-integral simulations qualitatively show very good agreement compared to the results obtained with the FPWPM in Figs. 5(e) and (f). 1D-cross sections through the focal spots in (f) and (h) are plotted in Fig. 5(i) and confirm the good agreement. The FPWPM predicts a marginally wider focal spot and marginally more energy in the higher order sidelobes than the RW-integral. This might be caused by the more realistic, full propagation of the electric field through the microlens in case of the FPWPM, taking into account lens aberrations and thus slightly reducing the quality of the focal spot. In contrast, the RW-integral does not offer full propagation of the electric field through the simulation volume. Instead, since we refrained from introducing aberrations manually, the lens is assumed to create a perfectly spherical wavefront and therefore an ideal focal spot. Further, the focal intensity cross section obtained with the scalar WPM strongly deviates due to neglection of polarization. As a consequence, vectorial treatment is mandatory for accurate simulation of the focal intensity distribution due to the high NA as well as the polarized illumination. Therefore, the FPWPM complements the scalar WPM as an established simulation routine for 3D-printed microoptics, in case polarization treatment is required for accurate results.

2.4 Evaluation of simulation runtime of the FPWPM

In this section, we investigate the runtime of the FPWPM in comparison to different simulation algorithms used in computing the prior simulation examples from sections 2.3.3 and 2.3.4. Table 1 offers an overview of the results. For our simulations, an AMD Ryzen 7 5700G CPU and 32 GB RAM were used. The FPWPM-algorithm was implemented as a plugin written in C++ for the open source software ITOM [25]. The FFTW-library [26] was used for numerical calculations of fast Fourier transforms. Further reading about the integration of our algorithm into ITOM and about the modeling of the refractive index distribution can be found in [15], where a previous algorithmic implementation of the scalar WPM-algorithm is presented.

Table 1. Comparison of runtime for simulation methods and simulation examples investigated in this work.^a

View Table

In the left column of Table 1, the runtimes for simulation of the grating from section 2.3.3 with RCWA, FPWPM, VWPM and the scalar WPM are listed. The runtime of the FPWPM is lower than the runtime of the RCWA by two orders of magnitude. Furthermore, the preparation of the simulation model for the RCWA is time consuming and not straightforward, while for the FPWPM, a simple matrix can be used to characterize the spatial distribution of the refractive index. Therefore, as long as slight loss in amplitude accuracy can be tolerated, the faster and more user friendly of both simulation routines can be found in the FPWPM. For comparison, runtime of the VWPM can be estimated by multiplying the runtime of the FPWPM with the scaling factor $N_xN_y/\left [\log (N_x)\log (N_y)M^2\right ]$ (with $M = 2$), assuming identical level of parallelization. The runtime of the FPWPM is significantly shorter than the estimated runtime of the VWPM.

Concerning the microlens in section 2.3.4, the two dimensional simulation for quick estimation of the focal position takes approximately ${15}\;\textrm{s}$ for $dy = dz = \lambda /15 \approx {0.042}\;\mathrm{\mu}\textrm{m}$. Further, the runtimes for the full 3D simulations are compared in the center and the right column in Table 1. For coarse sampling with $dx = dy = dz = {0.127}\;\mathrm{\mu}\textrm{m}$, the runtime of the FPWPM is ${53}\;\textrm{min}$, which poses a significant reduction compared to the estimated runtime of ${2.4}$ years for the VWPM. Finer sampling step size of 0.05 µm leads to about 12.4 h runtime for the FPWPM, while the VWPM is estimated to take approximately 173 years. Hence, FPWPM offers a simulation runtime which is acceptable in real-world application, while the VWPM for simulation of microoptics as large as in this example would not be feasible. The RW-integral, on the one hand, offers a comparably fast way to compute the focal plane of the microlens within a few minutes runtime, since simulation can be restricted to a very small area. On the other hand, as mentioned in section 2.3.4, the validity of the method is restricted to computing tho focal area instead of performing a full electric field propagation. Moreover, even if it would be physically sensible, the extrapolated runtime of the RW-integral would exceed hundreds of years for simulating the whole volume. To summarize, from the methods investigated in this work, the FPWPM is the only offering full propagation of a polarized electric field through the simulation volume within reasonable runtime.

In both the simulation of the grating and the microlens, the runtime of the FPWPM exceeds the runtime of the scalar WPM about fourfold. This can be seen as a tradeoff between considering the polarization for higher accuracy as well as extended validity for the FPWPM, and faster runtime for scalar WPM. Since the runtime remains in the same order of magnitude for both simulation routines, approximately the same size of simulation volume can be treated with the FPWPM algorithm as with the scalar WPM in practical application.

3. Conclusion

In this work, the FPWPM as a fast vectorial simulation algorithm for evaluation of complex microoptics was developed based on the VWPM-algorithm. In simulation examples, the runtime of the FPWPM was reduced by several orders of magnitude compared to the VWPM or the RCWA. The accuracy of the FPWPM for simulating typical 3D-printed microoptics was investigated and compared to analytical references, rigorous methods or other well-known simulation routines. We obtained very good agreement in the cases which are relevant for design of 3D-printed microoptics, such as the simulation of high frequency diffraction gratings and refractive microlenses. The feasibility of the FPWPM for simulation of large volumes with several hundreds of microns per dimension and visible wavelengths was demonstrated. Runtime proved to be viable for simulation of 3D-printed microoptics in context of the rapid prototyping cycle. This improvement in simulation capabilities enables consideration of polarization as a degree of freedom in optical design of microoptics. The performance of future optical designs can thus be enhanced, e.g., by considering high NA applications and using tailored polarization of illumination light.