Fast vector wave optical simulation methods for application on 3D-printed microoptics

Abstract. Additive manufacturing of microoptics using two-photon-lithography has been a rapidly advancing field of technology. Striving for ever more sophisticated optical systems prerequisites the access to appropriately fast and accurate wave-optical simulation methods to predict their optical performance. A simulation routine, which has been proven well suited for simulation of a vast range of 3D-printed microoptical systems, is the wave propagation method (WPM). Nevertheless, limitations in applicability remain due to the restriction on scalar electromagnetic fields, which prohibits consideration of polarization and thereby also the calculation of backward reflection at optical interfaces. Capabilities for design and analyses are, therefore, impaired for 3D-printed optical systems using those properties as key features in their design. As a first step to overcome those limitations, we presented new simulation methods based on the WPM in previous publications, extending its applicability toward simulation of vector electric fields, while maintaining short-simulation runtime. We focus on elaborating the practical application and integration of previously presented simulation methods in the design of complex 3D-printed optical systems. With it, we demonstrate the consideration of polarization and backward reflections in simulations far beyond paraxial and thin element approximations.


Introduction
2][3] Taking advantage of the versatile 2PL process, manufactured microoptical systems quickly increased in complexity to maximize their optical performance by demonstration of a highly miniaturized spectrometer 4 or optics for quantum dot fiber coupling. 5Inevitably, routines had to be developed to allow for design and simulation of those microoptical systems in order to allow for increasingly sophisticated optical designs.
7][8] Although ray optics can be sufficiently accurate to predict the performance of 3D-printed microoptics in numerous cases, overall small system size can cause considerable diffraction effects.Furthermore, quantities, such as focal spot size and shape, diffraction efficiency of gratings and polarization are wave effects inherently and thus necessitate consideration of wave optics for their assessment.For these cases, in which the validity of ray optics breaks down, wave optical simulation methods have been established for analyses to complement the design process of 3D-printed microoptics.Short-simulation runtime in the range of seconds up to a few hours at maximum is essential to allow for practical integration of simulation routines into the rapid prototyping approach of 3D-printed microoptics.To this end, the scalar wave propagation method (WPM) has been considered well fit for wave optical simulation of 3D-printed microoptics, offering a good compromise between accuracy, applicability beyond the paraxial-and thin-element approximation, and quick simulation runtime. 9,10][13] Despite its merits, potential application cases of 3D-printed microoptics remain, which cannot be covered using the WPM.More precisely, polarization generally has to be neglected in the WPM due to restrictions on scalar electromagnetic waves.Furthermore, the WPM is regarded as unidirectional, i.e., only light that propagates into the positive direction of the optical axis can be considered, neglecting light propagating in the opposite direction, i.e., backward reflected light, entirely.][16][17] In a first step to overcome those limitations, we have previously developed the fast polarized wave propagation method (FPWPM), extending the WPM toward consideration of polarized vector waves 18 while maintaining a short simulation runtime.Further, extension of the algorithm toward consideration of backward propagating vector waves has led to the derivation of the bidirectional wave propagation methods (BWPM). 19The feasibility of the FPWPM and the BWPM for fast and accurate simulation of typical microoptics manufactured using 2PL has been demonstrated, e.g., by simulation of polarized light focused by a high numerical aperture (NA) microlens and simulation of a catadioptric microoptical system.
In this paper, we present a methodology for a practical application of the BWPM for evaluation and the performance prediction of complex microoptical systems.To introduce the basic concepts of the algorithm, we first summarize the algorithms of the FPWPM and the BWPM.Then using two case studies, we demonstrate the application of the BWPM to support the development cycle of 3D-printed microoptics.In the first case study, we investigate the coupling of illumination light into an imaging fiber bundle (IFB) used for single-channel illumination and imaging.For the reduction of disturbing back reflected light at the IFB interface, 3D-printed microoptics are designed and assessed using the BWPM.Measurement data of the real manufactured optical elements are incorporated into the simulation model of the BWPM to predict their performance as realistically as possible.The simulation results allow for thorough analysis of aberrations and conclusions for further procedure in the development process can subsequently be drawn.In the second case study, we investigate a single-mode fiber (SMF)-based microoptical probe with potential applications, e.g., in Raman and Brillouin spectroscopy, fluorescence, and optical coherence tomography imaging 20 or as a Fabry-Pérot sensor. 16,17,21he probe utilizes a 3D-printed microlens to illuminate a sample and collect back reflected light as a measured signal through the same SMF.We perform batch simulations in order to predict the probe's signal strength and investigate issues such as the required sampling for appropriate resolution of the simulated signal.
2 Summary of the FPWPM and the BWPM 2.1 FPWPM

Summary of the FPWPM algorithm
The FPWPM has been derived from the vector wave propagation method (VWPM) 22 in order to allow for consideration of polarization in wave optical simulation of 3D-printed microoptics.Contrary to the VWPM, great emphasis on an efficient algorithmic implementation during the derivation allows for sufficiently short-simulation runtime to effectively apply the FPWPM as a part of the rapid prototyping manufacturing approach in 3D-printed microoptics.This could be achieved by adapting an approach allowing for an efficient numerical implementation, presented in Ref. 9, for the use in the VWPM.In this paper, the algorithm of the FPWPM will be briefly summarized to facilitate understanding of later sections.For a detailed derivation, the reader may be referred to Ref. 18.
The FPWPM carries out a step-wise, forward-directed propagation of a known incident vector electric field.As illustrated in Fig. 1(a), the simulation environment is defined by a discrete refractive index distribution ñðr ⊥ ; zÞ, assigning a complex refractive index, denoted by the tilde, to a total number N x × N y × N z of discrete volumetric pixels (voxels) with d x × d y × d z extension.For the simulation of 3D-printed microoptical systems, the refractive index data of commercial photoresists are regarded from the measured data from Ref. 23. Per definition, the axis z constitutes the optical axis, along which a step-wise, forward-directed propagation is carried out.The vector r ⊥ ¼ ðx; yÞ denotes the spatial coordinate in the transversal plane, the vector k ⊥ ¼ ðk x ; k y Þ denotes the transverse components of the wave vector k.The spatial distribution of the electric field at a discrete position along the z axis is given by E ð2Þ ðr ⊥ ; zÞ.Here the number two in the superscript denotes a reduced tensor that only contains the x and y components of an inherently three-dimensional tensor, omitting the third dimension, i.e., its z components.
Taking advantage of the Fourier transform F and the inverse Fourier transform F −1 , the FPWPM performs the propagation step of an electric field E ð2Þ ðr ⊥ ; zÞ at the position z by a discrete distance d z to the next position z þ d z along the optical axis by solving the equation: e m p : i n t r a l i n k -; e 0 0 1 ; 1 1 7 ; 5 2 5 r;s ðk ⊥ ÞF fE ð2Þ ðr ⊥ ; zÞgg: (1) The tensor T E Q -T A R G E T ; t e m p : i n t r a l i n k -; s e c 2 . 1 . 1 ; 1 1 4 ; 6 8 7 T ; t e m p : i n t r a l i n k -; e 0 0 3 ; 1 1 4 ; 6 0 7 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 4 ; 1 1 4 ; 5 7 1 The propagation factor P s ðk ⊥ Þ considers the phase progression of the electric field when propagated by the distance d z .The amplitude transmission coefficients t TE and t TM for transverse electric (TE)-and transverse magnetic (TM)-polarization, respectively, can be obtained by solving Fresnel's equations.
In order to allow for efficient computation of Eq. ( 1) by taking advantage of fast Fourier-and inverse fast Fourier transforms, all dependencies on the transverse spatial coordinate r ⊥ had to be eliminated upon the derivation of T ð2Þ r;s ðk ⊥ Þ.Since dependencies on r ⊥ were primarily caused by the inhomogeneity of the refractive index distribution ñðr ⊥ ; zÞ, as indicated in Fig. 1(a), this issue could be solved by decomposing the refractive index distribution inside the transverse plane at the positions z as well as z þ d z into a number of homogeneous, spatially independent subregions.To this end, we derived the weight function: In Eq. ( 1), M denotes the number of distinct refractive index values present inside the simulation volume, while summation over M using the indices of summation r and s ensures that all combinations of homogeneous subregions before and behind the interface are accounted for.A graphical illustration of the decomposition of the refractive index distribution in the simple case of M ¼ 2 can be observed in Fig. 1(b), which results in four possible combinations for the material at the interface to be summed up.

Runtime considerations and impact on practical applicability of the FPWPM
Compared to the VWPM algorithm as its direct predecessor, the efficient formulation of the FPWPM reduces the computational effort for a propagation step from OðN 2 x × N 2 y Þ to OðM 2 ½N x logðN x Þ × N y logðN y ÞÞ.Instead of exhibiting quadratic scaling with the number of data points in the transversal plane, the FPWPM scales close to linearly to the number of data points in the simulation volume.For the simulation of a typical 3D-printed microoptical system, this reduces the expected simulation runtime from up to hundreds of years for the VWPM to just a few minutes to hours for the FPWPM.For this consideration, we assumed one homogeneous material for the printed optics and a homogeneous surrounding, hence M ¼ 2, and a size of a few hundred microns up to a millimeter per dimension.Due to this decisive reduction in simulation runtime, application of the FPWPM to consider polarization in wave optical simulation of 3D-printed microoptical systems are feasible.

Motivation
Basic concepts have been presented, which extend unidirectional step-wise propagation methods, such as the FPWPM toward bidirectional propagation, allowing for consideration of backward-directed reflected electric fields.Yet, available bidirectional simulation methods suffer from slow computation unfeasible for practical application on 3D-printed microoptics. 24Others are sufficiently fast but based on the beam propagation method (BPM) and therefore restricted to scalar electromagnetic waves.This disallows consideration of polarization when calculating the reflectivity coefficients, 25 which would be mandatory for accurate results.Furthermore, upon exceeding the paraxial approximation, the simulation results of the BPM turn out to be imprecise. 22Therefore, the presented bidirectional algorithms are limited either in their applicability or accuracy.The availability of the FPWPM, which allows for consideration of polarization far beyond the paraxial approximation, therefore naturally paves the way for the development of a fast and accurate BWPM, applicable to 3D-printed microoptics.In the following, we will summarize the algorithm of the BWPM and refer readers to Ref. 19, for detailed description and in-depth analyses of the presented simulation method.

Summary of the BWPM algorithm
Analogously to Eq. ( 1) from the FPWPM, the equation E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 6 ; 1 1 7 ; 5 4 3 serves to compute propagated electric fields E ð2Þþ ðr ⊥ ; z þ d z Þ in the BWPM, the superscript positive sign indicating forward-directed propagation.To compute the backward reflected electric field in each propagation step, we assume the existence of a matrix of reflection R ð2Þ r;s ðk ⊥ Þ, which solves the vector Helmholtz equation at a plane interface for reflected electric fields, complementary to the matrix of transmission T ð2Þ r;s .The backward propagating reflected electric field E ð2Þ− ðr ⊥ ; zÞ, with the superscript negative sign denoting backward-directed propagation, then holds E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 7 ; 1 1 7 ; 4 0 9 For R r;s ðk ⊥ Þ, we derived the explicit expression E Q -T A R G E T ; t e m p : i n t r a l i n k -; s e c 2 . 2 . 2 ; 1 1 7 ; 2 9 5 R The BWPM requires a case sensitive definition of the propagation factor P ðiÞ r;s ðk ⊥ Þ denoted by the superscript ðiÞ.The distinct cases will be discussed in Sec.2.3.1.Although traversing the simulation volume in the positive direction along the optical axis, both equations serve to locally calculate the forward propagating electric field and the reflected electric fields occurring at interfaces between two distinct media for each voxel, as depicted in Fig. 2(a).The reflected electric fields are temporally stored at each point in space for later use, which momentarily requires a restriction on two-dimensional simulations due to demands on system memory (RAM).
Upon reaching the end of the simulation volume, the direction of propagation is reversed, and propagation is subsequently carried out in the reverse direction.Previously stored backwardreflected electric fields now contribute to the incident electric field during backward propagation, as depicted in Fig. 2(b).Since Maxwell's equations are invariant to Lorentz transformation, the propagation equations for backward-directed propagation remain similar in shape compared to Eqs. ( 6) and (7).Attention must be paid to insert backward propagating electric fields into those equations correctly to account for backward directed propagation, obtaining E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 0 ; 1 1 4 ; 4 7 5 and E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 1 ; 1 1 4 ; 4 2 3

Choice of the propagation factor
As hinted at in Sec.2.2.2, a case sensitive choice of the propagation factor P ðiÞ r;s ðk ⊥ Þ must be implemented in order to ensure that the electric field is always assumed at the same position inside of a voxel inside the simulation grid.Thus phase errors upon coherent superposition of the fields are avoided.We arbitrarily choose to always place the electric fields at the furthest possible position toward the positive direction of the z axis inside each voxel.To this end, four distinct cases must be distinguished in the BWPM, which are denoted by the superscript ðiÞ and listed in Table 1.The four cases and the effect of the four propagation factors on the field propagation on the level of a single pixel are depicted in Fig. 2(a).

Iterative procedure
Traversing the simulation volume back and forth once with the BWPM algorithm will subsequently be referred to as one iteration.Upon returning to the beginning of the simulation volume, the algorithm can either be truncated if the residual error is considered to be sufficiently small, or another iteration can be initiated, as depicted in Fig. 2(b).The electric fields calculated in each iteration are coherently superposed to obtain the total electric field.Systematic investigations have shown that as little as one single iteration can already be sufficient to approximate the magnitude of backward-directed reflected fields reasonably well, with residual errors just in the range of a few percent. 19Concerning transmitted electric fields, a minimum of two iterations is required for any effect, while around three iterations yield accurate results in many cases.Upon incidence on an interface between two media, reflected electric fields are locally stored.An iterative approach is subsequently used to propagate those electric fields back and forth, until sufficiently small residual error permits truncation of the algorithm.

Application on Complex 3D-Printed Microoptical Systems
In the following section, we present a methodology for practical application of the BWPM in the simulation of complex 3D-printed microoptical systems.In two application examples, the BWPM will be used to predict the performance of optical systems.Furthermore, we showcase the use of the BWPM in analyses by incorporating measured topographies of manufactured microoptics into the simulation model.

Case Study 1: Manipulation of Illumination Light
Coupling into an Imaging Fiber Bundle

Motivation and previous investigations
In endoscopy, reducing the size of an endoscopic probe is a key to minimize the invasiveness of an endoscopic intervention.To this end, an IFB poses an ideal candidate for the construction of a flexible, submillimeter-diameter endoscope, capable of imaging and illumination.Particularly, a small probe diameter can be achieved using the same IFB for illumination as well as imaging as a single-channel endoscopic probe.7][28] As depicted in Fig. 3(a), the air-fiber interface at the proximal facet of the IFB, i.e., the facet facing the illumination and imaging objective lens, is considered to be particularly delicate in the construction of a single-channel illumination/imaging endoscope.Back reflected straylight originating at the proximal IFB interface directly superposes signal light guided by the IFB and both are subsequently imaged sharply to the imaging sensor, causing potentially critical reduction in imaging contrast.
To reduce loss of contrast due to noise from back reflection at the proximal IFB facet, we applied the BWPM during the development of an approach to add 3D-printed optical elements at the facet of the IFB. 19As illustrated in Fig. 3(b), those elements provide refractive index matching and are therefore referred to as index matching caps (IMCs).The IMCs are designed to shift the back reflection from the IFB facet toward the facet of the IMC.As a consequence, back reflected light would not be imaged sharply to the sensor anymore and a confocal pinhole could be used to filter the reflected noise, allowing only the signal light from the IFB facet plane to pass.In the simulation, the BWPM predicted considerably less noise caused by back reflection at the IFB facet due to the IMCs, matching our experimental observations. 19evertheless, consideration had been restricted to light at the optical axis in this earlier case study, omitting more complex off-axis aberrations due to imperfectly manufactured IMCs.The manuscript at hand aims to use the BWPM to provide in-depth performance analyses of the manufactured IMCs, by consideration of off-axis illumination light and measured topography data in simulations.Results of those analyses subsequently offer a basis for a decision on the further manufacturing procedure of the IMCs, i.e., whether their performance can be considered satisfactory or whether iterative shape correction steps must be undertaken.

Extending simulation toward off-axis fields
First, the simulation environment for the BWPM is set up by modeling the polished facet of an FIGH-10-350S (Fujikura Ltd., Japan) IFB with a diameter of 350 μm, an average core diameter of 2 μm and a core pitch of 5.2 μm, based on data from Chen et al. 29 The illumination light with a wavelength of λ ¼ 0.635 μm is defined as a converging spherical wavefront, the opening angle matching the NA of the single cores of the IFB.Beyond the on-axis illumination field F 0 , we add two telecentric off-axis fields F 1 and F 2 with lateral offsets of h 1 ¼ 52 μm and h 2 ¼ 104 μm in respect to the center of the IFB.The offsets are chosen to ensure perfect coupling of the incident fields into an off-axis IFB core, respectively.Polarization is chosen to be linear with the polarization axis tilted with respect to the observation plane, resulting in equal energy in TE-and TM-polarization, mimicking the energy balance of unpolarized light.For visualization, the absolute value of the TE-polarized part of the electric field is depicted in Fig. 3(c).The forward propagating part of the electric field is coupled well into the individual fiber cores.The backward propagating part of the electric field shows that a significant portion of the illumination light is reflected at the IFB facet, exactly mimicking light emitted by the fiber cores.In Fig. 3(d), an IMC consisting of the photoresist IP-Q (Nanoscribe GmbH, Germany) with a thickness of 50 μm is added to the proximal IFB facet.The forward propagating part of the electric field remains almost unaffected by the presence of the IMC, i.e., light still can be coupled well into the cores of the IFB.Spherical aberrations due to the refractive index mismatch caused by the presence of the IMC seem to be negligible, only minor axial repositioning of the IFB facet by about 17 μm is required to match the best focus position for ideal coupling.The overwhelming part of the back reflected illumination light originates at the IMC's facet, while no visible reflection at the IMC-to-IFB interface remains.The light reflected at the IMC facet is defocused in respect to the IFB facet, i.e., a confocal pinhole is suitable for a separation of reflected illumination light and signal emitted at the IFB facet.Therefore, the BWPM predicts that a 50 μm thick IMC can enable illumination through an IFB, while maintaining high image contrast for backward transported imaging signal.

Simulation of aberrations due to shape deviations in the manufactured IMC
][32] Although the viability of IMCs has already been demonstrated in experiment for on-axis illumination, 19 we expect those shape aberrations to affect and perturb off-axis fields more severely, likely negatively impacting the manufactured IMCs' performance.For thorough analysis of this case, the topography of the manufactured IMC, as depicted in Fig. 4(a), is measured with a confocal microscope (MarSurf CM Explorer, Mahr GmbH, Germany).The measured topography, depicted in Fig. 4(b), indeed exhibits notable deviation from its designed flat shape.Incorporating the measured topography into the simulation model of the BWPM, we obtain the forward propagating field depicted in Fig. 4(c).Although the on-axis field is barely affected by the aberrations and coupling into the IFB remains strong, coupling is diminished very clearly for the off-axis fields.For the backward propagating field, depicted in Fig. 4(d), slight residual backscattering of light from the IFB facet can be observed by closely looking at the IFB-IMC interface, especially for the field F 2 with its comparably weak coupling to the IFB core.This leads to the conclusion that optimum light coupling into the IFB cores is not only essential for resourceful use of illumination light, but also impacts the amount of backscattered light.Optimization of light coupling hence also promises reduction of backscattering at the IMC-IFB interface.However, the overwhelming majority of back reflection still occurs at the IMC facet.Thus the spatial separation between reflected illumination light and signal emitted by the IFB cores remains intact in order to allow for spatial filtering using a confocal pinhole.Therefore, we focus further investigation to the coupling of illumination light into the IFB.
To fathom the aberrations causing the diminished coupling efficiency for the off-axis fields, the intensity of forward propagating light close to the IFB facet inside the IMC is calculated and depicted in Fig. 5(a).Two effects can be observed prominently: due to the curved IMC's additional refractive power, the plane of best focus is curved, as illustrated with a dashed purple line.The off-axis foci hence appear increasingly distant from the IFB facet.The defocus between the best focus and IFB facet is z ¼ −1.08 μm for the field F 1 and z ¼ −4.72 μm for F 2 , resulting in radius of curvature of r ≈ 1020 μm for the best focal plane, reducing coupling efficiency to the cores of the flat IMC facet.Further, the image field is slightly distorted, i.e., the IMC's additional refractive power perturbs the magnification of the illumination system.As a consequence, the illumination light misses the IFB cores in case the image distortion is not accounted for.In practice, this issue can be mitigated by scaling the axial offsets h 1 and h 2 of F 1 and F 2 with a suitably chosen factor D ¼ 1.02.Thereby, the distortion caused by the IMC is considered by adjusting the magnification of the telecentric illumination system.A BWPM simulation with the predistorted input fields in Fig. 5(b) suggests that coupling into the IFB can be improved considerably for off-axis fields, although at a lower performance than for the on-axis field due to the remaining focal plane curvature.
Another quantity with decisive impact on the coupling efficiency is the quality of the focus of the illumination light.For visualization and investigation, lateral cross sections through the simulated intensity in the focal planes without IMC, with ideal IMC, and with the manufactured IMC are depicted in Figs.5(c)-5(e) for the three considered fields, respectively.The intensity is normalized to the peak intensity of the simulated focus without IMC, which can be considered perfectly diffraction limited.For the simulations without and with ideal IMC, the focal plane is assumed to be at the first layer of material, right behind the IFB facet.For the simulation of the manufactured IMC, the layer with the smallest lateral focus extension is considered for each field, respectively.In this case, the lateral displacement of the focus position due to image distortion can be observed clearly in Figs.5(d) and 5(e).Moreover, the quality of the focal spot increasingly deteriorates for the off-axis fields F 1 and F 2 .For quantification, we use Strehl's ratio, i.e., the ratio between the peak intensity of the real and the ideal focus, with a Strehl's ratio of SR ≥ 0.8 commonly considered as diffraction limited. 33Since the plots are normalized to the peak of the ideal focal spot obtained from the simulation without IMC, Strehl's ratio directly equals the peak intensity of the respective focus.Although for F 0 and F 1 , the simulated focus with the manufactured IMC remains above or at least close to the diffraction limit, Strehl's ratio falls short of the threshold SR ≥ 0.8 for F 2 and loss of coupling efficiency due to deterioration of the focal spot has to be expected.

Summary and conclusions for further studies
To summarize the findings from the analyses based on the BWPM simulations from Sec. 3.1.3,the manufactured IMC is well suited when restricting the IFB's image circle to a diameter of around 100 μm.Beyond, the shape deviation of the IMC is expected to cause considerable aberrations and therefore losses due to diminished coupling between IFB and imaging objective.Moreover, the BWPM suggests that this loss in coupling efficiency also leads to an increase of back scattered light at the IFB facet, reducing the effectivity of the IMC.Hence, in case an application requires a larger image circle diameter, the simulations highly encourage shape correction 32 of the manufactured IMC to provide higher fidelity to the IMC's design.The simulations suggest special emphasis on the correction of the focal plane curvature as well as the focal spot quality for off-axis fields.This maximizes coupling between illumination objective and IFB cores and minimizes back scattering of illumination light, thereby maximizing the contrast for imaging through the IFB.

Case Study 2: Signal Prediction of a Single-Mode Fiber-Based
Single-Channel Sensing Probe

Introduction
Subject of our second case study is a 3D-printed microoptical system based on an SMF.This probe ought to serve as a point illumination system, which simultaneously collects backward reflected or scattered light as a measured signal through the same SMF.Illumination and signal collection through a single SMF allows for ultra-compact realization of a sensing probe with potential applications, e.g., in fluorescence imaging and OCT, 20 as well as Raman-or Brillouin spectroscopy.

BWPM simulation
In this section, we showcase the use of the BWPM to investigate and predict the performance of a probe as outlined in Sec.3.2.1.To this end, the probe's design is depicted in Fig. 6(a).It consists of a 780 HP (Thorlabs GmbH, United States) SMF, with a piece of coreless fused silica fiber (CFSF) spliced on top to allow for mode expansion.The SMF is modeled with a ratio of n clad ¼ 0.996 • n core to obtain NA ¼ 0.13 according to the specifications. 34At the facet of the CFSF, a 3D-printed aspheric lens with NA ¼ 0.4 is attached, designed to tightly focus light emitted by the SMF core into the vacuum behind the lens.Therefore, the probe acts as a point illumination system.Simultaneously, back-reflected light from the focal plane of the probe is coupled backward into the core of the SMF, acting as a collected signal to be evaluated at the proximal facet of the SMF.We restrict considerations to monochromatic light, i.e., illumination as well as the collected signal share the wavelength λ ¼ 0.78 μm.Simulation discretization is chosen d y ¼ d z ¼ λ∕ð2 • n max Þ to satisfy the sampling theorem for the comparably large simulated region of 135.5 μm × 1256 μm.To generate a signal, we place a dielectric target with a refractive index n ¼ 1.5 in the focal plane of the probe.In a batch simulation with a total of 101 single simulations and approximately 90 min cumulative simulation runtime on an AMD Ryzen 7 5700G CPU, we gradually displace the target relative to the probe's focal plane over a range of Δz ¼ AE6 μm in steps of 0.12 μm.For visualization, Fig. 6(b) depicts the simulated forward propagating illumination light for a target displacement of Δz ¼ 6 μm, and Fig. 6(c) shows the associated reflected electric field, i.e., the collected signal.In this particular case, only minimal coupling to the SMF core can be observed.For each value of Δz, we integrate the simulated intensity of the back coupled light into the SMF core and calculate its ratio to the energy of input illumination light.The resulting integrated energy as a function of the target displacement is depicted in Fig. 6(d).
As expected, the strongest signal is predicted in cases in which the target is placed close to the focal plane of the probe.In this scenario, the 3D-printed asphere focuses back reflected light from the target interface directly onto the SMF core, ensuring high back coupling efficiency.The displacement of the target, independent of its sign, leads to a defocus, declining back coupling efficiency and therefore a decrease in signal strength.The resulting envelope of the collected energy curve in Fig. 6(d) is Gaussian-like.This Gaussian-like envelope curve is notably superposed by a higher frequency carrier function.Harsh edges in the curve indicate undersampling of the carrier function.To thoroughly investigate the carrier function, another batch simulation is carried out, with a smaller target displacement step size of 0.067 μm to achieve finer resolution.In order to resolve such a fine displacement in the simulation model, we likewise set d z ¼ 0.067 μm.To avoid excessive simulation runtime of the batch simulation, we restrict on a region of interest (ROI) of Δz ¼ AE1 μm.The result of the second batch simulation for the ROI is plotted in Fig. 6(e).Now, the carrier function appears to be well resolved.The period width of the carrier function is around 0.39 μm ¼ λ∕2.Therefore, the oscillation is likely caused by the interference between reflected light originating at the photoresist-vacuum interface and the vacuum-target interface, respectively.Since the target displacement Δz leads to a change in optical path length by 2Δz, an interference pattern with a period with of λ∕2 is generated, requiring twice the sampling rate for appropriate resolution compared to the simulation wavelength.This finding explains why fulfilling the sampling theorem, which is commonly considered to ensure appropriate resolution in wave optical simulation, has not been sufficient to resolve the carrier function in our first batch simulation depicted in Fig. 6(d).For the future case studies that investigate interference patterns of electric fields with the BWPM, we suggest fulfilling the sampling theorem at least twice as an initial choice for simulation discretization, i.e., using a sampling step size d ≤ λ∕ð4 • n max Þ.Even finer sampling might be advisable in cases in which more than two electric fields contribute to the interference pattern to be resolved.
With sufficient resolution, the amplitudes of the carrier function and the envelope curve can now be analyzed properly.The envelope curve has a maximum value of around 0.04 at z ¼ 0 μm, while the carrier function oscillates with an amplitude of around 0.015.The knowledge about the mechanism behind the formation of the carrier function enables targeted tuning of its amplitude in the optical design.Larger amplitude can be achieved by optimizing for ideal coupling of back reflected light from the 3D-printed lens-vacuum interface into the SMF core, lower amplitude likewise by minimizing said quantity.

Summary of the findings and implications of this case study
In the preceding case study, we used the BWPM to investigate an SMF-based microoptical system that simultaneously serves for illumination and collection of backward-reflected signal light.A first batch simulation revealed the shape of the expected signal curve, however, it exhibited insufficient resolution.Further investigation allowed us to narrow down the cause of the insufficient resolution to the interference of different reflected electric fields, leading to interference patterns with a period width smaller than the simulation wavelength.Based on this finding, we were able to propose sampling suggestions for BWPM simulations in order to avoid such resolution issues in future studies.
Further, the simulations reveal possible strategies for optimization of the probe's signal strength, depending on its purpose in practical application.For example, the probe could be utilized as a Fabry-Pérot-sensor by interpreting the amplitude of the carrier function as a signal.In this case, the amplitude of the carrier function should be optimized as outlined in Sec.3.2.2.A second option would be the interpretation of the backward reflected light from the target as a signal, as would be the case in an OCT probe.In that case, the carrier function would act as a disturbance of the signal and its amplitude should be minimized.

Outlook: Applicability of the FPWPM on 3D-Printed-Inhomogeneous
Refractive Index Media Recently, by taking advantage of carefully tailored modulation of the laser dose during the 2PL printing process, local variation in the refractive index of 3D-printed structures has been demonstrated. 35,36This technique, referred to as (3 + 1)D laser writing, is promising for the fabrication of waveguides, volume holograms, and GRIN lenses.Furthermore, a similar effect of exposure-modulated variation of the refractive index occurs as a side effect in the emerging field of grayscale 2PL, where the variation of laser dose is primarily used for dynamic adjustment of the printing voxel size 37,38 in order to improve printing speed and surface quality.In all those cases, wave optical simulation is desirable to predict the influence of the inhomogeneous refractive index on the performance of manufactured optical systems.
In respect to the underlying physics, the FPWPM is well suited to treat those types of simulation problems.To begin with, spatial variation of the refractive index can quickly and easily be considered in the creation of the simulation environment, i.e., the definition of the refractive index distribution ñðr ⊥ ; zÞ, as outlined in Sec.2.1.1.Furthermore, the calculus that has been subsequently presented in Eqs. ( 1)-( 5) fully supports inhomogeneous refractive index media.Nevertheless, a large number M of distinct refractive index values can pose a challenge regarding simulation runtime, since the calculation effort per propagation step is proportional to OðM 2 ½N x logðN x Þ × N y logðN y ÞÞ.If no further measures are taken, the simulation runtime can such become unacceptable for practical application when considering a GRIN lens, which has to be described by a large number of distinct local refractive indices M. As a compromise to alleviate this issue, appropriate rounding of the refractive index could reduce the required number of distinct refractive index values with negligible effect on the simulation accuracy.More importantly, a large number of M can also be very well and comparably easily addressed by massive parallelization of the simulation algorithm, since calculations are mutually independent for each distinct refractive index value.An acceptable downside would be that requirements would likely exceed the means of a workstation computer and hence must be provided by a cluster.Overall, for future research, the application of the FPWPM on microoptics manufactured by (3 + 1)D laser writing poses a promising field.

Summary and Conclusion
In this work, we showcased the application of the BWPM to simulate microoptics in a system context, which could not be assessed with customary simulation routines like the WPM before.To this end, we carried out two case studies on 3D-printed microoptical systems with regard to back-reflections and possibilities for their mitigation and targeted manipulation.
In the first case study, we investigated the performance of optical elements 3D-printed at the tip of an IFB, providing refractive image matching, in order to mitigate back reflection of illumination light during coupling into the IFB.Simulations with the BWPM validate this concept for the ideal IMC design.Regarding the measured topography of a manufactured IMC in the simulation, the BWPM predicts notable aberrations for off-axis fields.Simulations further imply that those aberrations not only reduce the coupling efficiency of illumination light into the IFB but also increase the amount of back-scattered light at the IFB-IMC facet.Consequently, in this first case study, the BWPM allowed us to gather detailed information on the presumable performance of a reflective microoptical system.Using that information allows for a solid decision on the further procedure in the development process, i.e., whether the system performance is considered sufficient, or optimization is necessary to achieve acceptable performance.
In the second case study, we investigated an ultra-compact 3D-printed microoptical probe, which combines illumination and signal collection through just one SMF.In the course of this study, we observed the superposition of different reflected electric fields in simulations, causing interference patterns with subwavelength period width.As a consequence, fulfilling the commonly considered sampling theorem proved insufficient in order to properly resolve these interference patterns.We were therefore able to propose sampling suggestions for bidirectional simulations in order to help avoid such issues in the future.Solving the sampling issue enabled analyses of the interference pattern of the reflected electric fields and strategies could be proposed for targeted optimization of the probe.Thus the maximum signal-to-noise ratio can be striven for the use as a Fabry-Pérot-sensor, as well as an OCT-or fluorescence imaging probe.

Disclosures
The authors declare no conflicts of interest.

2 n 1 < n 2 xFig. 1
Fig.1(a) Illustration of simulation model constitution in the FPWPM algorithm.Refractive index data of the simulated microoptics are defined on a discrete grid through which a step-wise propagation of an electric field is carried out along the z axis.(b) Possible material combinations at an interface for a number of M ¼ 2 distinct refractive indices.Summation over the indices r and s is required to consider all possible material combinations in front of and behind the interface between the planes z and z þ d z , respectively, for the definition of the weight function I r ;s ðr ⊥ Þ.

Fig. 2
Fig.2(a) Consideration of transmitted and reflected electric fields in the BWPM at the single-pixel level.The arrows indicate the effects on the electric field for application of the distinct four cases of the propagation factor, indicated by the superscript i.(b) Consideration of transmitted and reflected electric fields in the BWPM at the top level of the simulation.Upon incidence on an interface between two media, reflected electric fields are locally stored.An iterative approach is subsequently used to propagate those electric fields back and forth, until sufficiently small residual error permits truncation of the algorithm.

¼ 1 3
ðk ⊥ Þ 1 Propagation of the transmitted electric field during forward directed progression through the simulation volume e ik z;s dz 2 Propagation of the reflected electric field during forward directed progression through the simulation volume e 0 Propagation of the transmitted electric field during backward directed progression through the simulation volume e ik z;r dz 4 Propagation of the reflected electric field during backward directed progression through the simulation volume e 2ik z;r dz

Fig. 3
Fig. 3 (a) Setup for single-channel illumination/imaging through an IFB.For areal imaging, the focused laser beam is scanned pointwise over the IFB facet via a telecentric laser scanning confocal microscope.Back reflected illumination light, which is considered as noise, and signal collected from the IFB spatially overlap at the IFB facet, decreasing imaging contrast.(b) Adding a 3D-printed IMC at the fiber facet results in spatial separation of signal and noise, allowing for filtering of the noise via the confocal pinhole of the microscope.(c) BWPM simulation of the TE-polarized electric fields during coupling into a polished IFB.The white arrow indicates the propagation direction of the electric field, the dashed and solid white lines indicate the outlines of the IFB and its single cores, respectively.Although jE TE j þ reveals good coupling between IFB and illumination objective lens, jE TE j − indicates significant back reflection at the IFB facet.(d) Adding a 50 μm thick 3D-printed IMC at the IFB, no visible back reflection at the IFB facet remains.

Fig. 4
Fig. 4 (a) Light microscopic image of the manufactured IMC at the facet of an IFB.(b) Topography of the manufactured IMC, obtained by confocal microscope measurement.(c) and (d) BWPMsimulation, incorporating the measured topography of the manufactured IMC into the simulation model.The white arrows indicate the propagation direction of the electric fields.The white dashed and solid lines indicate the outlines of the IFB, its cores, and the IMC, respectively.

Fig. 5
Fig. 5 BWPM simulation and analyses of aberrations caused by shape error of the manufactured IMC.(a) Intensity of the forward propagating light.The purple dashed line indicates the curved plane of focus for off-axis fields.The white dashed lines indicate the outlines of the IFB and IMC, and the white solid line outlines the cores of the IMC.The white arrow indicates to direction of light propagation.(b) Intensity of the forward propagating light with precompensation of image distortion.Compared to (a), overlap between off-axis foci and IFB cores and therefore light coupling is increased notably.(c)-(e) Comparison of simulated foci for all three field positions without IMC, with ideal IMC and with measured manufactured IMC.The plots depict lateral cross sections along the y axis.For the measured IMC, the cross sections are taken at the z position with the best focus, respectively.

Fig. 6
Fig. 6 Design and simulation of a SMF-based single-channel illumination/imaging probe.(a) Structure of the probe.Close to the focal area of the aspheric lens, a displaceable target structure is placed, surrounded by vacuum.(b) Simulated forward propagating electric field for a target displacement of 6 μm.(c) Simulated reflected electric field.The white arrow indicates the propagation direction of the depicted electric field.(d) Integrated energy collected by the core as a function of the target displacement relative to the probe's focal plane.The envelope of the collected energy curve can be well observed, but the simulation resolution is too coarse to properly resolve a higher frequency carrier function.(e) Restricting on an ROI to avoid excessive simulation runtime, a higher simulation resolution allows to resolve the carrier function well.

Table 1
List and description of the four distinct propagation factors used in the BWPM, taken from Ref. 19.