1. Introduction

Two photon lithography (2PL) has become a well established method for additive fabrication of complex micro optical devices and systems [1–3], as it offers sub micron feature size and good optical surface quality. Simulation of these optical devices can prove to be challenging, since wave optical effects regularly cannot be neglected due to small feature size and therefore raytracing is rendered inaccurate, but total system size simultaneously prohibits rigorous electromagnetic treatment due to runtime and hardware constraints. To bridge this gap, wave optical simulation algorithms like the scalar wave propagation method (WPM) [4] and the vector wave propagation method (VWPM) [5] have been developed and later optimized regarding runtime [6,7]. In short, they trade some degree of accuracy and loss in general applicability compared to rigorous methods, e.g., by restricting solely on forward propagation, for far lower runtime, less demand on hardware and ease of use [8]. The WPM has been well established for simulation of a wide range of 3D-printed micro optics [9–13].

Nevertheless, especially the restriction on forward propagating waves for the WPM remains a harsh constraint in case that backwards propagating, i.e., reflected electromagnetic waves, have to be considered. This would, e.g., be the case for 3D-printed parabolic micro mirrors, which are fabricated using a combination of 2PL and silver coating techniques [14,15]. Moreover, backward reflected light plays a crucial role in the design and evaluation of 3D-printed common path OCT-systems for medical imaging [16] as well as 3D-printed Fabry-Perot sensors used, e.g., for acoustic wave detection [17], refractive index [18] or temperature sensing [19], just to name a few examples. Thus, wave optical simulation considering reflected light would facilitate the analysis of such optical devices and enhance the capabilities to carry out targeted design processes.

Bidirectional beam propagation methods (BPMs), that can consider backward reflection and subsequent backward propagation of light, are available in literature, e.g., [20]. Yet, they suffer from typical limitations which motivated the development of the WPM in the first place, mainly the restriction on small angles, and thus are not well suited for simulation of 3D-printed micro optics. Further, the basic framework for a bidirectional VWPM has been reported [21]. While the accuracy of this method meets our requirements very well, simulation runtime is beyond acceptable for treatment of realistic simulation volumes occurring in 3D-printed micro optics. To the best of our knowledge, no runtime optimized implementation is available at the moment to overcome this issue. Hence, in this work, we use the bidirectional VWPM as well as the runtime optimized fast polarized wave propagation method (FPWPM) [7] as a basic framework to derive a fast bidirectional wave propagation method (BWPM). We present the underlying electromagnetic calculus as well as considerations regarding algorithmic implementation of the bidirectional propagation scheme.

Moreover, we study the convergence of the BWPM and demonstrate its applicability in simulation of practical problems. Thorough investigation is conducted to study the simultaneous use of an imaging fiber bundle for illumination as well as for imaging of a sample. While this configuration is desirable to minimize probe diameter as much as possible, e.g., in keyhole access endoscopy, its practical implementation is not straightforward, as backwards reflected illumination light at both facets of the imaging fiber superimpose and dominate the measured imaging signal. In section 2.2, we present the concept of 3D-printed index matching caps (IMCs) at the facets of the imaging fiber for targeted suppression of backwards reflected light, overcoming limitations of other approaches available in literature. Evaluation of the IMC design is decisively aided by the use of the BWPM. Following up the BWPM simulations, we conduct experiments and demonstrate the feasibility of noise suppression at imaging fibers using the previously designed 3D-printed IMCs.

2. Results

2.1 Bidirectional FPWPM-algorithm

2.1.1 Summary of the FPWPM

In this section, the FPWPM as the algorithmic basis of the BWPM is briefly summarized. For detailed description, readers may be referred to [7].

The FPWPM carries out a step wise propagation of the x- and y-components of an incident vector electric field $\mathbf {E}^{(2)}(\mathbf {r}_\bot,z)$ through a discretized simulation volume along the z-axis. The vector $\mathbf {r}_\bot = (x,y)$ expresses the transverse spatial coordinate, $\mathbf {k}_\bot = (k_x,k_y)$ the transverse spatial frequency coordinate with the absolute value $\vert \mathbf {k}_\bot \vert = k_\bot$. The electric field is propagated from the position $z$ to the position $z+dz$ using the propagation equation

$\mathscr {F}$ and $\mathscr {F}^{-1}$ represent the Fourier-operator and the inverse Fourier-operator. $\mathbf {T}^{(2)}_{r,s}(\mathbf {k}_{\bot })$ is a matrix that solves the vector Helmholtz equation for an incident electric field at a plane interface between the media present at the position $z$ and the position $z+dz$. The factor $\mathscr {P}_{s}(\mathbf {k}_{\bot })$ considers the phase progression of the electric field for a propagation by a discrete distance $dz$. $t_{TE}$ and $t_{TM}$ are the amplitude transmission coefficients for transverse electric (TE) and transverse magnetic (TM) polarized electric fields, which can be obtained from Fresnel’s equations [22]. $M$ is the total number of distinct refractive indices present inside the simulation volume, $r$ and $s$ denote indices for summation. During summation, the index $r$ iterates over all the $M$ distinct discrete refractive indices present inside the simulation volume at the position $z$, the index $s$ likewise for $z+dz$. Therefore, the indices $r$ and $s$ are associated with the discrete refractive indices $\tilde {n}_s$ and $\tilde {n}_r$, $k_r$ and $k_s$ denote the magnitude of the wave vector in those media, respectively. The tildes indicate that the refractive indices can assume complex values. $I_{r,s}(\mathbf {r}_{\bot })$ is a function that takes into account the spatial distribution of the refractive index distribution $\tilde {n}(\mathbf {r}_{\bot })$ in the transverse XY-plane.

Crucially, the dependence on the spatial coordinate $\mathbf {r}_{\bot }$ and the frequency coordinate $\mathbf {k}_{\bot }$ is partitioned towards the separate functions $I_{r,s}(\mathbf {r}_{\bot })$ and $\mathbf {T}^{(2)}_{r,s}(\mathbf {k}_{\bot })$. This has been achieved by suitable derivation of the weight function $I_{r,s}(\mathbf {r}_{\bot })$, which holds

Using this weight function, the matrix $\mathbf {T}^{(2)}_{r,s}(\mathbf {k}_{\bot })$ can be formulated to solely depend on the transverse frequency coordinate $\mathbf {k}_{\bot } = (k_x,k_y)$, as can be seen in Eq. (2). As a consequence, the inverse Fourier-operator can be applied in Eq. (1) instead of a full numerical integration over $k_x$ and $k_y$, allowing for very efficient numerical implementation.

The $z$-component of the electric field is derived from the $x$- and $y$-components in Fourier space using

2.1.2 Derivation of the BWPM

To obtain the BWPM, the FPWPM can be extended to consider reflection at interfaces and propagation of backwards directed electric fields. In analogy to Eq. (1), the forwards propagating electric field is calculated using the propagation equation

The positive sign in $\mathbf {E}^{(2)+}(\mathbf {r}_{\bot },z+dz)$ denotes propagation in the positive direction of the $z$-axis, while a negative sign denotes backwards directed propagation. To obtain the reflected electric field $\mathbf {E}^{(2)-}(\mathbf {r}_{\bot },z)$, it is reasonable to assume that a matrix of reflection $\mathbf {R}^{(2)}_{r,s}(\mathbf {k}_{\bot })$ can be found to substitute the matrix of transmission $\mathbf {T}^{(2)}_{r,s}(\mathbf {k}_{\bot })$ in Eq. (6), satisfying the equation

In order to enable fast algorithmic implementation of the BWPM, it is mandatory that $\mathbf {R}^{(2)}_{r,s}(\mathbf {k}_{\bot })$ complies with the efficient computation method we introduced for the FPWPM in [7], i.e., is only dependent on the transverse frequency $\mathbf {k}_{\bot }$. To the best of our knowledge, no such matrix of reflection is known in literature. As accomplished for the matrix of transmission for the FPWPM, we derive $\mathbf {R}^{(2)}_{r,s}(\mathbf {k}_{\bot })$ by solving the vector Helmholtz equation at a plane boundary between $z$ and $z+dz$, computing the backward reflected electric field from the incident electric field using methods from [22–24]. We obtain

The propagation factor $\mathscr {P}^{(i)}_{r,s}(\mathbf {k}_{\bot })$ considers the phase progression of the electric field during propagation. Its explicit form as well as the implications of the superscript (i) will be discussed later. For every propagation step, the calculated reflected electric field $\mathbf {E}^{(2)-}(\mathbf {r}_{\bot },z)$ has to be temporally stored, since it is later going to contribute towards the electric field propagating into the opposite direction. In the most general implementation as used in this work, this requires storing $N_x \times N_y \times N_z$ electric field vectors. The resulting memory consumption necessitates restriction on simulation of two dimensional or small three dimensional volumes. Since reflected fields can only arise in presence of a physical interface, storing reflected fields only at those points in space would substantially reduce memory requirements. However, this would require careful and extensive adaptation of our algorithmic rountines presented in [8] and therefore is beyond the scope of this work. We consequently restrict the simulations conducted in this work to two dimensions, assuming the third dimension to be infinitely extended.

Upon reaching the end of the simulation volume, the propagation direction is reversed in order to carry out propagation in the negative direction along the $z$-axis. Since Maxwell’s equations are invariant to Lorentz transformation [25], propagation equations for backwards propagating waves remain similar in shape. In that sense, backwards propagating electric fields $\mathbf {E}^{(2)-}(\mathbf {r}_{\bot },z)$ in respect to the global coordinate system can be interpreted as a forwards propagating electric fields in a local coordinate system with inverted main propagation direction, i.e., $z$-axis. Thus, propagation equations for backwards propagating electric fields can be obtained from Eqs. (6) and (7) by inverting the sign of the superscripts and adjusting the positions addressed along the $z$-axis:

Starting the backwards propagation, the backwards directed electric field is initialized as zero. The incident electric field at each propagation step is then constituted by the sum of the backward propagating electric field and the reflected electric field at the current $z$-position, which had been calculated and stored during the last forward propagation. Upon reaching the beginning of the simulation volume again, the algorithm can either be terminated or another back and forth propagation can be initiated. One cycle of back and forth propagation will be further referred to as one iteration. In Fig. 1(a), a simulation volume consisting of three distinct media is used in order to visualize the appearance of backward an forward propagating electric fields as well as their summation during back and forth propagation in the BWPM. In order to enable accurate summation of the backward- and forward propagating electric fields to obtain the total electric field, the positioning of the electric fields at each simulation pixel must be consistent. Therefore, special attention has to be dedicated towards the propagation factor $\mathscr {P}_{r,s}^{(i)}(\mathbf {k}_{\bot })$ and four cases requiring different propagation factors have to be distinguished, which are denoted by the superscript $i = \{1,2,3,4\}$. Those four cases are listed in Table 1 and are illustrated in Fig. 1(b) for three consecutive pixels inside the propagation volume, each assigned a distinct refractive index. The propagation factors are chosen to ensure that the electric fields always are positioned at the furthest to the right of each pixel.

 

Fig. 1. (a) Visualization of the addition of forward- and backward propagating electric fields for simulation of a stack of three materials during the first two iterations of the BWPM. (b) Illustration of the possible cases $i = \{1,2,3,4\}$ of the phase propagation factor $\mathscr {P}_{r,s}^{(i)}(\mathbf {k}_{\bot })$ and the respective covered distance. Three pixels in a row with distinct refractive indices are depicted. The three vertical dashed lines indicate the placement of the electric fields inside each pixel.

Download Full Size | PDF

Table 1. Explicit form and case description of the four distinct propagation factors applied in the BWPM.^a

View Table

Concerning simulation runtime, a single propagation step of the BWPM requires approximately twice the computational effort compared to one step of the FPWPM, since the two propagation Eqs. (6) and (7) have to be evaluated instead of just one single Eq. (1). Moreover, a full iteration of the BWPM necessitates traversing the simulation volume twice. Overall, the computation time of one iteration of the BWPM adds up to take about four times the runtime of the FPWPM. Regardless, the scaling behaviour of both methods remains identical, therefore the BWPM can be applied to approximately the same size of simulation volumes as the FPWPM.

2.1.3 Studying convergence of the BWPM for simulation of a dielectric slab

Over the number of iterations, the BWPM converges towards its final solution. For practical use, an appropriate cutoff threshold for the number of iterations has to be determined as a compromise between simulation runtime and acceptable residual error. To this end, we study the convergence of the BWPM for simulation of oblique light incidence onto a dielectric slab as depicted in Fig. 2. Closed analytical solutions [26,27] are available as reference to calculate the transmittivity and reflectance of the slab, which are used to determine the residual error of the BWPM-simulation, depending on the number of iterations used.

The simulated setup is illustrated in Fig. 2(a). The slab is 2.5 µm thick, has a refractive index of $n = {1.5}$ and is surrounded by vacuum. It is illuminated obliquely from the left with infinitely extended plane waves ($\lambda = {1}\;\mathrm{\mu}\textrm{m}$) propagating in the YZ-plane with main propagation axis along $z$ and incident angles up to $ {80}^{\circ}$. Along the $y$-axis, perfectly transparent boundary conditions are imposed. Using the BWPM with only one single iteration as depicted in 2(b), reflectance curves already approach the theoretical reference. Since no secondary reflections which would contribute to the forwards travelling energy have been taken into account with only one iteration, the transmittivity curves are comparably inaccurate. Yet, in case that only the reflectivity is of interest, as little as one single iteration suffices to obtain a useful approximation of the solution. The evolution of residual absolute error size will be later discussed in detail. Three iterations as shown in 2(c) already agree very precisely with the theory for the majority of the angular range. Therefore, three iterations can be considered an appropriate cutoff for most practical cases, even when high accuracy is desired. Using ten iterations as depicted in 2(d), no visible deviation remains. However, considering the already small error size remaining after the third iteration, the absolute impact of error reduction is marginal when using more iterations and has to be carefully weighted against the increasing simulation runtime.

 

Fig. 2. Simulation of transmittivity and reflectance of a dielectric slab for oblique illumination. (a) Illustration of the geometry of the simulated slab. (b)–(d) Simulation of the slab using the BWPM with 1, 3 and 10 iterations, respectively. Convergence notably increases with number of iterations. (e) and (f) Rate of convergence over the number of iterations for the angles marked with arrows in (b) for TE- and TM-polarization, respectively.

Download Full Size | PDF

To more thoroughly investigate the rate of convergence as a function of the incident angle and number of iterations, we choose the incident angles $ {5}^{\circ }, {20}^{\circ}, {32.5}^{\circ}, {51}|^{\circ }$ and $ {64}^{\circ }$, which accompany distinctive features on the simulated transmittivity and reflectance curves, marked with small arrows in Fig. 2(b). Results are depicted in Fig. 2(e) + (f) for TE- and TM-polarization, respectively. As criteria to determine the residual error, we use $\varepsilon = 1 - \vert \text {R}+\text {T}\vert$, which represents the residual deviation of transmittivity and reflectance from energy conservation, omitting the sign. While in theory, deviations of transmittivity and reflectance with opposite sign could cancel out each other and therefore wrongfully satisfy energy conservation, studying the error in transmittivity and reflectance separately proved this effect to be insignificant. Consequently, the deviation from energy conservation is found to be a meaningful quantity to estimate the rate of convergence for this simulation method.

With increasing number of iterations, all curves decline close to linearly on the logarithmic scale, indicating that the relative error reduction per iteration is close to constant for each angle individually. At around $10^{-13}$ residual error, the curves level out, likely due to rounding errors or machine precision. It is further evident that the higher the absolute initial error, the slower the relative convergence rate, i.e., the flatter the convergence curves. The initial error magnitude as a function of the incident angle cannot generally be predicted prior to simulations, since it heavily depends on the shape of the simulated as well as the theoretical curves of the first iteration as depicted in Fig. 2. Those again are strongly impacted in a nontrivial way by parameters like the thickness of the slab and its refractive index.

As previously proposed, three iterations as a cutoff point results in residual error well below $ {1}{\%}$ except for the single angle $ {64}^{\circ }$ with particularly large initial error. Therefore, it can be confirmed as a reasonable number of iterations.

2.1.4 Simulation of a complex catadioptric micro optical system

Feasibility of fabrication of reflective 3D-printed micro optics such as silver coated parabolic mirrors has been previously reported [14,15]. For those micro optics, simulation capabilities are extended by the BWPM beyond rigorous methods or raytracing as used in the works previously mentioned. This particularly facilitates simulation of reflective optics inside of more complex integrated micro optical systems.

For demonstration, we designed a micro optical system consisting of an optical fiber delivering TE-polarized illumination light with $\text {NA} = {0.4}$ and $\lambda = {0.78}\;\mathrm{\mu}\textrm{m}$, followed by a microlens for collimation, which is constituted of the photoresist IP-S (Nanoscribe GmbH, Germany). For lens optimization, the software Zemax Optic Studio was used. The lens surface type is spherical with conic constants. The optimized radii are $r_1 = {0.098}\;\textrm{mm}$ and $r_2 = {-0.201}\;\textrm{mm}$, the conic constants are $c_1 = {-3,494}{}$ and $c_2 = {-4,248}{}$ and the thickness of the lens is $t = {0.076}\;\textrm{mm}$. After passing through the lens, the collimated light is then focused by a parabolic reflector, which was calculated using the function $z = y^2/(4f')$ with the focal length $f' = {50}\;\mathrm{\mu}\textrm{m}$ [28]. The parabolic reflector consists of a silver coating with $ {400}\;\textrm{nm}$ thickness, the complex refractive index $\tilde {n}_{Ag} = 0.213+6.697i$ and is modeled to be placed on a substrate consisting of IP-S with thickness $d = {5}\;\mathrm{\mu}\textrm{m}$. For simulation discretization, we chose $dy = dz =\lambda /(10\cdot \vert {\tilde {n}_{max}}\vert )\approx {0.052}\;\mathrm{\mu}\textrm{m}$, while the x-axis is assumed to be infinitely extended. As a result, $N_x \times N_y \times N_z = 1 \times 5403 \times 8711$ data points constitute the simulation volume with approximately $ {280}\;\mathrm{\mu}\textrm{m}$ lateral and $ {451.5}\;\mathrm{\mu}\textrm{m}$ axial extension. Due to the high reflectivity of the silver mirror, only one iteration is deemed necessary for the BWPM simulation. Relying on routines thoroughly described in [8], model creation can be quickly and easily carried out by automated extraction of optical design parameters from the Zemax Optic Studio file and creating a two dimensional matrix containing the refractive index data, using the programming language python. The refractive index positional data of the parabolic reflector as well as the reflector’s substrate were subsequently integrated into this matrix. Simulation runtime is ${170}\;\textrm{s}$ using an AMD Ryzen 7 5700G CPU. Figure 3(a) shows the simulated forward propagating electric field emitted by the fiber and collimated by the microlens. Upon incidence on the silver mirror, the electric field is completely reflected, hence no further forward propagating field remains. Figure 3(b) shows the backwards propagating electric field focused by the silver mirror, while Fig. 3(c) shows the total electric field. Figure 3(d) allows for close inspection of the focal region of the parabolic reflector. Interference of the forward and backward propagating fields gives rise to a fringe pattern. The focus is located $ {49.85}\;\mathrm{\mu}\textrm{m}$ from the mirror vertex, close to the design focal length $f' = {50}\;\mathrm{\mu}\textrm{m}$ of the parabolic reflector.

 

Fig. 3. Simulation of a catadioptric micro optical system, consisting of a light source, collimating lens made of IP-S ($n = {1.505}{}$) and a silver-coated parabolic mirror. (a) Absolute value of the forward propagating electric field. (b) Absolute value of the backward propagating electric field. (c) Absolute value of the total electric field. (d) Magnified image of the focal region of the parabolic mirror as simulated in (c).

Download Full Size | PDF

2.2 Simulation and experimental evaluation of 3D-printed index matching caps to reduce back reflection noise in imaging fibers

2.2.1 Basic concept of index matching caps

In the following, as briefly introduced in section 1, we use the BWPM to study light reflection at the facets of an imaging fiber and to design of 3D-printed optical components, so called index matching caps (IMCs), for targeted suppression of back reflected light. The underlying problem is illustrated in Fig. 4(a) and (b). Coupling light into the imaging fiber leads to backwards reflection at the proximal fiber facet, i.e., the fiber facet facing the coupling objective, due to the change of refractive index at the interface between air and the core material. This backwards reflected light carries no imaging information and can be thus considered noise. Meanwhile, the successfully coupled light propagates inside the fiber up to the distal fiber facet, i.e., the fiber facet facing the sample to be imaged. Upon exiting the fiber at the distal facet, again a part of the light is reflected backwards due to the index contrast between fiber core and the environment. Like the backwards reflected light at the proximal fiber facet, it acts as a source of noise. The portion of light leaving the fiber subsequently is used to illuminate the sample. Scattered and backwards reflected light from the sample is collected and re-coupled into the distal facet of the imaging fiber and propagates back to the proximal fiber facet. In order to generate an image which can be recorded by a sensor, the proximal fiber facet serves as an intermediary image plane, from which light is collected by the same objective lens originally used for illumination. The single fiber cores of the imaging fiber act as pixels of the intermediary image. The collected light can finally be imaged onto a sensor. Unfortunately, the back reflection takes place directly at the intermediary image plane at the proximal fiber facet and further, the propagation directions of back reflected light and light emitted by the fiber core are identical, thus both are imaged sharply onto the sensor. Especially in biological applications, the refractive index contrast between sample and aqueous environment is typically around one order of magnitude lower [29,30] than the index contrast between air and the fiber core of around $\Delta n \approx {0.5}{}$. Therefore, the reflected noise reduces the imaging contrast up to the point of completely dominating the signal, prohibiting the use of imaging fiber bundles in this setup.

 

Fig. 4. (a) Illustration of light coupling between a single core of a bare imaging fiber and an objective lens for epi-illumination. To illuminate the whole sample plane as well as to generate a full image, pointwise scanning over the whole fiber facet is necessary. The arrows indicate the propagation direction of the wavefront. (b) For a bare fiber, light is reflected directly at the fiber facet, i.e., the intermediary image plane, and thus superimposes the signal transported from the distal end of the fiber. (c) Light coupling into an imaging fiber with an 3D-printed IMC at the fiber facet. (d) The reflection takes place at the facet of the IMC and is thus spatially separated from the intermediary image plane.

Download Full Size | PDF

To address this issue, several methods have been proposed to dispose of the disturbing reflected noise. For example, some authors demonstrated an angled polish of the proximal fiber facet and oblique mounting of under yet another angle in respect to the optical axis of the illumination light, combined with glycerine for index matching [31]. The required manual adjustments render this approach sensitive and prone to errors and prohibit quick and easy mounting of imaging fibers inside an experimental setup. Other authors use an index matching oil to cover the gap between the optics and both fiber facets [32,33]. While this approach is decisively more straightforward, it requires application of a potentially hazardous chemical at both facets of the imaging fiber, which could prohibit any biomedical applications. Critically, oil immersion objectives matching the numerical aperture (NA) of the fiber cores in the range of $\text {NA} \approx {0.3}{}- {0.4}{}$ are usually not commercially available and had to be custom made for the mentioned work. Further approaches like single core illumination and wide field detection [34] require no further manual preparation and are well suited in principle, but are obviously restricted to cases where those illumination conditions can be met.

As an alternative approach, we propose 3D-printed index matching caps (IMCs) at both fiber facets, while using a commercially available dry objective without cover glass correction for light coupling and imaging of the proximal fiber facet. As depicted in Fig. 4(c) and (d), when coupling light into a single core of the imaging fiber, the majority of back reflection takes place at the air facet of the IMC, while index matching to the fiber core leads to negligible reflection between IMC and fiber core. Illumination of the whole sample area as well as image generation can subsequently be carried out by telecentric pointwise scanning of the laser beam over the proximal fiber facet, a concept well known from laser scanning microscopy, which also has been demonstrated for an imaging fiber endoscope in [32]. For the fabrication of the IMC, biocompatible photoresists are available [35,36] which allow for particularly close index matching in the range of $\Delta n \approx {0.01}$. As a consequence of the close index matching, the intermediary image plane and the location of significant noise generation are spatially separated and spatial filtering using a confocal pinhole in the imaging beam path is enabled, allowing only the collected signal from the intermediary image plane to pass. Even in case no confocal pinhole is used, the noise would still be mitigated since the back reflected light from the IMC facet is out of focus and therefore not imaged sharply onto the sensor.

The thickness of the 3D-printed IMC has to be chosen carefully. On the one hand, it has to be large enough to ensure that the reflected light is shifted out of focus sufficiently far. On the other hand, a low thickness of the IMC is desired to minimize degradation of focus quality and therefore coupling efficiency due to refractive index mismatch caused by the IMC when using the air corrected objective lens. We use the BWPM to predict the feasibility of an IMC regarding those two contrary requirements in simulation. Furthermore, we use the BWPM to determine the efficacy of the IMC in suppressing back reflected light at the proximal fiber facet.

2.2.2 Wave optical simulation

We first use the BWPM to perform wave optical simulation of coupling illumination light into an unmodified imaging fiber. In this work, we use FIGH-10-350S (Fujikura Ltd., Japan) fiber bundles with an average core diameter $d_{\text {Core}} = {2}\;\mathrm{\mu}\textrm{m}$ and $\text {NA}_{\text {Core}} = {0.39}{}$ according to [37]. Monochromatic light with a wavelength $\lambda = {0.635}\;\mathrm{\mu}\textrm{m}$ and TE-polarization is used in simulation. Discretization is chosen $dy = dz = \lambda /(5 \cdot n_{max})\approx {0.085}\;\mathrm{\mu}\textrm{m}$. The forwards as well as backwards propagating electric fields are depicted in Fig. 5(a) and (b), respectively. For visual clarity, only one single fiber core is simulated at the center of the optical fiber. Considering the adjacent cores has been tested in simulation and does not impact the simulation results significantly. As mentioned in section 2.2, the image is generated by scanning the laser focus over the whole proximal fiber facet. Assuming telecentric illumination, results obtained for the central fiber core of the imaging fiber can be assumed to be valid for all cores over the whole fiber facet. Fig. 5(a) shows that the forward propagating field is successfully coupled into the fiber core. According to the V-parameter, the cores of the imaging fiber are expected to be multi-moded for $\lambda = {0.635}\;\mathrm{\mu}\textrm{m}$, i.e., $V \approx {3.86}{} > {2.405}{}$ [38]. Excitation of more than one fiber mode explains the oscillating electric field pattern observed inside the fiber core. As can be seen in Fig. 5(b), reflection takes place at the interface between air and the fiber core. Since this reflected electric field directly originates in the focal region of the objective lens, it closely resembles the incident electric field, but with inverted propagation direction. More critically, it thereby also closely resembles incoming signal transported from the distal facet of the imaging fiber, emitted from the core at the proximal fiber facet and then collected by the objective lens. The back reflected electric field and the incoming signal from the distal fiber facet superimpose, with the back reflected electric field acting as a source of noise.

 

Fig. 5. BWPM-simulation of light coupling into a single fiber core of an imaging fiber. (a) Bare fiber, forwards propagating electric field. (b) Bare fiber, backwards reflected electric field. (c) Fiber with IMC, forwards propagating electric field. (d) Fiber with IMC, backwards reflected electric field. The white dashed lines outline the imaging fiber and indicate the shape of the IMC in (c) and (d). Further, white dashed lines indicate the shape of the central fiber core in (b) and (d). In (a) and (c), the outlines around the core are omitted to allow for clear observation of the electric field propagating inside the fiber core. (e) Simulated point spread function of the air corrected objective without IMC, with ideal IMC and with the measured surface of the fabricated IMC. Simulated Strehl ratio remains $ {0.97}{}$ for the ideal and $ {0.91}{}$ for the manufactured IMC (see section 2.2.3). (f) Intensity of simulated back reflected light along the $z$-axis for the bare fiber and a fiber with an ideal IMC.

Download Full Size | PDF

In comparison, simulation is carried out with a 3D-printed IMC added on top of the imaging fiber facet in order to prevent reflections directly at the fiber facet. The IMC is $d = {50}\;\mathrm{\mu}\textrm{m}$ thick and made of the photoresist IP-Q (Nanoscribe GmbH, Germany) with $n \approx {1.511}{}$, closely matching the refractive index of the fiber core with $n \approx {1.5}{}$. The simulated forward propagating field in Fig. 5(c) suggests that the light is still coupled well into the fiber core. Meanwhile, analysis of the reflected electric field as depicted in Fig. 5(d) shows that the focus of back reflected light is shifted by $\Delta z\approx {83}\;\mathrm{\mu}\textrm{m}$ away from the intermediary image plane at the fiber facet. According to the intensity of the backwards propagating electric field along the $z$-axis as plotted in Fig. 5(f), using 3D-printed IMCs reduces back reflection at the fiber facet from approximately $ {4}{\% }$ to well below 0.1 % of the incoming light. Cross sections through the intensity at the focal plane of the objective lens at the fiber facet, indicated with white arrows in a) and c), are plotted in e). They confirm that the quality of focus remains nearly ideal with a Strehl ratio of 0.97 % despite spherical aberrations introduced by the refractive index mismatch due to the IMC. Therefore, the BWPM simulation suggests that the IMC is capable of strongly suppressing the back reflection at the intermediary image plane, while coupling of illumination light into the fiber core remains largely unaffected by the presence of the IMC.

2.2.3 Experimental investigation

Since using IMCs proved viable in simulation, experimental investigation is carried out using a measurement setup as shown in Fig. 6. An illumination laser source with $\lambda = {635}\;\textrm{nm }$ is collimated by lens 1 (L1) and focused onto the proximal facet of the imaging fiber using a microscope with magnification $m_p = {20}{\times }$, consisting of a tube lens and an objective lens. The objective lens $\text {NA}_{Obj} = {0.4}{}$ ($ {20}{}{\times }$ Mitutoyo Plan Apo NIR) is chosen to closely match the NA of the individual fiber cores. The distal fiber facet is imaged to a camera chip by another microscope with $m_{d} = {20}{\times }$, which enables measurement of transmission through the fiber for light coupled into the proximal fiber facet. At the proximal side of the imaging fiber, back reflected light from both fiber facets propagates backwards through the microscope setup, is collimated by lens 2 (L2) and incides onto a 50:50 beamsplitter. Here, a portion of the incident light is redirected towards lens 3 (L3) and focused on a pinhole with diameter $d_{PH} = {15}\;\mathrm{\mu}\textrm{m}$. Since $f'_{L2} = f'_{L3} = {50}\;\textrm{mm }$, L2 and L3 act as a 1:1 relay system. Light transmitted through the pinhole is reimaged on a highly sensitive photodiode (Thorlabs APD410A2) using lens 4 (L4).

First, a bare imaging fiber is mounted in the measurement setup, and light is coupled into a single core at the proximal end of the fiber. Movement of the proximal fiber facet along the $z$-axis as depicted in Fig. 6 is carried out using a Standa manual XYZ translation stage (Standa Ltd., Lithuania). With scale graduation of $ {10}\;\mathrm{\mu}\textrm{m}$, we deemed a displacement of $ {5}\;\mathrm{\mu}\textrm{m}$ between each measurement feasible. Figure 7(a) shows the measured diode signal as a function of displacement of the proximal fiber end along the $z$-axis, as well as the intensity emitted from the same core, which is simultaneously detected by the camera at the distal fiber facet. As expected according to prior simulation, both reflected light and intensity at the proximal fiber end peak off at the same $z$-position, i.e. when the fiber facet is in focus of the objective lens. Switching to configuration b) as depicted in Fig. 6, the camera is replaced by a laser which is coupled into the distal facet of the previously illuminated core of the imaging fiber. This laser acts as an artificial signal, imitating a signal collected with the distal facet of the fiber. As can be observed in Fig. 7(a), detected artificial signal using the photodiode peaks at the same $z$-position as both other signals, i.e. when the proximal fiber facet is in the focal plane of the objective lens. Therefore, back reflected light from the fiber facets and a real signal would be indistinguishable during imaging. Next, 3D-printed IMCs as described in section 2.2.2 are attached to both facets of the imaging fiber. The IMCs are fabricated using the Photonic Professional GT2 (Nanoscribe GmbH, Germany) 2PL printer. Microscopic images of an IMC at the tip of the imaging fiber are shown in Fig. 8(a) and (b). Due to typical shrinkage during the development process, the facet of the IMC deviates from the initial flat design and is slightly curved, as can be observed in Fig. 8(a). Based on confocal microscope measurements of the manufactured surface, a BWPM simulation of the fabricated system is carried out, obtaining a Strehl ratio of $ {0.91}{}$ for the focal spot of the objective lens, depicted in Fig. 5(f) in comparison to the ideal focal spot. For higher performance, an iterative correction of surface shape as reported in [39] can be carried out to increase the optical performance of the IMC. For this work, the minor decrease in quality of the focal spot is considered insignificant for the proof of concept of the IMCs, hence we refrain from shape correction at this point.

 

Fig. 6. Schematic depiction of the experimental setup for analysis of light transmission and reflection at an imaging fiber in this work. L1-4: Lens 1-4. BS: Beamsplitter. TL: Tube lens. OL: Objective lens. IF: Imaging fiber. Cam: Camera. Configuration (a) is used for measuring transmitted intensity at the distal fiber end using a camera and back reflected light at the proximal end of the imaging fiber with a photodiode. Configuration (b) is used to transmit an artificial signal through the fiber by coupling a laser source into the distal facet and measuring signal strength with the photodiode. All data is obtained as a function of $z$-position of the proximal fiber end, using a manual translation stage for displacement along the $z$-axis.

Download Full Size | PDF

 

Fig. 7. Measurements obtained using the setup depicted in Fig. 6. (a) Bare fiber. (b) Fiber with a 50 µm thick IMCs at each facet. The use of IMCs leads to a notable shift of the back reflected noise in regard to the plane of optimal coupling and signal detection.

Download Full Size | PDF

 

Fig. 8. Microscopic images of a 3D-printed IMC on an imaging fiber. (a) Side view, the facet of the IMC is slightly curved due to shrinkage of the material, typical for the 2PL fabrication process. (b) Top view, the imaging fiber facet and the individual fiber cores can be imaged sharply through the IMC, despite imperfections resulting from fabrication.

Download Full Size | PDF

In paraxial approximation, $\Delta z' = m^2 \Delta z$ can be used to relate the simulated focus offset $\Delta z = {83}\;\mathrm{\mu}\textrm{m}$ in object space caused by the IMC, as obtained in section 2.2.2, to a defocus of $\Delta z' \approx {33}\;\textrm{mm }$ in image space [28]. With focus spot length of approximately $2z_R = \frac {2\lambda }{\pi (\text {NA'})^2} \approx {1.01}\;\textrm{mm }$ in image space past L3, $\Delta z'$ is large enough to allow elimination of back reflected light from the facet of the IMC using a pinhole, while the artificial signal from the distal fiber facet is allowed to reach the sensor as desired.

For verification, the fiber with IMCs at both facets is mounted into the experimental setup and undertaken the same measurement procedure as the bare fiber. The measurement results are depicted in Fig. 7(b). The transmitted intensity and detected artificial signal peak at the same $z$-position, indicating that optimal coupling efficiency as well as optimal signal detection both take place when the fiber core is at the best focus position of the objective lens. Most notably though, the peak position of back reflected light is clearly shifted by about $ {50}\;\mathrm{\mu}\textrm{m}$ in regard to the ideal coupling position between objective lens and fiber core. Since the measured back reflected light is expected to peak as soon as the facet of the IMC is in the focal plane of the objective lens, the measured shift of $ {50}\;\mathrm{\mu}\textrm{m}$ agrees very well with the design thickness $d = {50}\;\mathrm{\mu}\textrm{m}$ of the IMC. As a consequence of this shift, the signal can be detected without noise contribution from back reflections originating at the fiber facets. This proves that the 3D-printed IMCs are viable for suppression of back reflected light in epi-illuminated imaging fibers and therefore enable simultaneous epi-illumination and imaging using the same fiber bundle.

3. Conclusion

In this work, we presented the BWPM, which is a fast simulation method capable of propagating vector wave electric fields in both directions along a main propagation axis inside of comparably large simulation volumes. The BWPM extends the capabilities for simulation of 3D-printed micro optics towards complex reflective optical systems, e.g., catadioptric systems and resonators, and also allows for investigation of reflected light, e.g., in interferometers and OCT systems.

Further, we successfully applied the BWPM to investigate the performance of adding 3D-printed IMCs at the facet of an imaging fiber in order to suppress backwards reflected light. Obtaining promising results in simulation, we subsequently prepared an imaging fiber with IMCs and evaluated the performance experimentally. The measurements proved the viability of using 3D-printed IMCs for noise suppression when using epi-illumination on imaging fibers. Based on these results, future work can combine IMCs and 3D-printed imaging optics at the tip of an imaging fiber to enable illumination and imaging with the same fiber, enabling a very low probe diameter well suited keyhole access endoscopy.