Flexible lensless endoscope with a conformationally invariant multi-core fiber

The lensless endoscope represents the ultimate limit in miniaturization of imaging tools: an image can be transmitted through a (multi-mode or multi-core) fiber by numerical or physical inversion of the fiber's pre-measured transmission matrix. However, the transmission matrix changes completely with only minute conformational changes of the fiber, which has so far limited lensless endoscopes to fibers that must be kept static. In this letter we report for the first time a lensless endoscope which is exempt from the requirement of static fiber by designing and employing a custom-designed conformationally invariant fiber. We give experimental and theoretical validations and determine the parameter space over which the invariance is maintained.


Introduction
Cellular-level microscopic imaging has long been a vital tool in biomedical research. Recent years have seen numerous efforts to miniaturize imaging instruments to enable cellular-level imaging in behaving animals. A recent example is the miniaturized head mounted microscope for fluorescence imaging [1] (2 g, 9×15×22 mm). This approach has light source, filters, imaging optics, and CMOS camera integrated into a head mounted device. A different approach bases the light delivery and collection on optical fiber which brings the advantage that light source and detectors can be remote rather than integrated in the head mounted device. As an added benefit, optical fiber-based approaches can use pulsed laser sources and perform non-linear imaging-to date demonstrations have been using piezo-electric actuators [2] (few grammes, Ø3 mm × 40 mm) or micro-electromechanical systems (MEMS) mirrors [3] (2.15 g, Ø10 × 40 mm) to perform point-scanning imaging. For an overview of the most "conventional" optical fiber-based approaches, see Ref. [4].
A new approach to fiber-based endoscopes came about in 2011 [5,6,7,8,9,10] which does away with imaging optics between fiber and sample and consequently is often termed "lensless endoscopes". This represents the ultimate limit in miniaturization, since the head-mounted device can be as small as an optical fiber itself. Cellular-level imaging in live mice by lensless endoscopes has recently been demonstrated [11,12], but imaging modality was restricted to one-photon fluorescence imaging and fiber length to few cm, and the animal subject was fixed.
The lensless endoscope is most commonly implemented with image acquisition by point-scanning [5,11,12,7,9], although other strategies exist [8,13,14]. The generic concept of the lensless endoscope is: (i) The transmission matrix (TM) in the basis of localized modes [15] of the fiber is measured in a preliminary step; (ii) The columns of the TM are identified as the input fields that give rise to focused output fields; (iii) The phase masks that convert the laser beam into said input fields are calculated; (iv) Said phase masks are displayed sequentially on a spatial light modulator (SLM), resulting in a two-dimensional scan of the output focus, all the while back-scattered fluorescence signal is recorded, as in a point-scanning microscope. The fiber can be either a multi-mode fiber (MMF) [11,12,8,13,16] or a multi-core fiber (MCF) [9,10,14]. In either case, its transmission matrix H can be thought of as where H 0 is the pre-measured TM, measured in a reference conformation, and X is the "extrinsic contribution" to the TM when the fiber conformation departs from the reference conformation, in general it is represented as an operator [16]. In both the MMF and MCF case, the great challenge remains that following each conformational change of the fiber either the additional extrinsic contribution (X) or the new TM (H) must be experimentally quantified whether directly [17,18] or indirectly [16] in order for aberation-free imaging to continue. This is the main obstacle standing between us and a flexible lensless endoscope which would open the possibility for minimally-invasive imaging in behaving animals.

Twisted multi-core fiber (MCF)
In our earlier work we have shown how lensless endoscopes based on MCF (as opposed to MMF) simplify many of the considerations pertaining to the TM [19]. In particular, in Ref. [20], we showed that the extrinsic contributionX to the TM of a MCF is simply a diagonal matrix with complex elements of unit norm and argument which is linear in the transverse coordinate, i.e. a matrix with only two free parameters. In the present article, we take an entirely new approach to overcoming this challenge: We seek to design a MCF for which the extrinsic contribution is identity (X = I) which would render real-time tracking of the TM unnecessary and finally allow lensless endoscopes to reach their full potential as the endoscope fiber could be allowed to flex freely. As we show below a MCF with a twisted geometry can fulfill this requirement. We designed and fabricated [see the Supplemental Information and Fig. 1(a)] two MCFs, "MCF1" and "MCF2" which have the following parameters in common: Number of cores N = 489; core-to-core distance (pitch) Λ = 16µm; parabolic core index profile with index difference relative to the cladding ∆n = 30·10 −3 ; core diameter D = 3.5 µm (single mode). MCF1 is twisted with a helical period of P 1 = 33 mm, MCF2 with P 2 = 8.2 mm. See Fig. 1(b) for a schematic of the twisted MCFs as well as the coordinates that we will employ. A non-twisted MCF, "MCF0" with the same parameters but P 0 = ∞ was drawn as well to serve as a reference.

Twisted MCF: Intrinsic properties
The theoretical properties of the twisted MCF are derived in the Supplemental Information, a summary is given in the following. Considering first a straight section of twisted MCF of length L. The center core is unaltered by the twist, the effective index of the fundamental mode in the center core at wavelength λ = 1 µm is n (0) eff = 1.4626. Generally, the core i, radially offset from the center by d (i) , has its physical length modified by the twist as As a consequence, the twisted MCF exhibits a radially-dependent native phase delay ∆φ (i) as an intrinsic property: By the same token, the twisted MCF also exhibits a radially-dependent native group delay ∆τ (i) as an intrinsic property In the strict sense, Eq. 4 is approximate, but in the Supplemental Information we demonstrate that group index and refractive index can be used interchangeably without invalidating our conclusions. Figure 2(a) shows the comparison between the native group delays measured in MCF1 and MCF2 by spectral interference [21] as a function of radial core offset d (i) and those calculated by Eq. 4. We observe excellent agreement. Figure 2(b) shows experimentally measured mode field diameters compared to those found by numerical simulation. Both show that the twist induces no appreciable variation in the mode field diameter. The phase profile of the off-center modes, however, is altered as a consequence of  the cores of the twisted MCF no longer being parallel to the MCF axis. This geometrical consideration predicts that the beam exiting core i has an azimuthal component, the angle of the free space beam with the MCF axis being [See Inset in Fig. 2 In Fig. 2(d) we compare the experimentally measured θ (i) , measured as the optimal coupling angle of the free space beam entering core number i, to those predicted by Eq. 5, and we observe excellent agreement. Since the off-center cores follow a curved path in space, they are expected to experience loss when P is small. Our numerical simulations showed that no curved path loss (< 1 dB/m) due to this effect is expected for neither MCF1 nor MCF2. Curved path loss is expected only for twist period P < 2 mm as detailed in the Supplemental Information. Nevertheless, as seen in Fig. 2(c) MCF2 does present a radiallydependent loss, but a cut-back measurement from 900 mm to 345 mm revealed that their origin was not the fiber itself, so we are able to attribute them to coupling loss due to mode mismatch between the input gaussian beam profiles and the MCF modes with non-flat transverse phase cf Eq. 5.

Twisted MCF: Extrinsic contribution
We consider now a bent section of the twisted MCF with length L and bend radius of curvature R cf Fig. 1. Now, in addition to the native phase and group delays, extrinsic contributions to the phase delay ∆[∆φ (i) ] and group delay ∆[∆τ (i) ] appear. As we demonstrate in the Supplemental Information, when L = k · P , k ∈ N , L (i) is independent of R, and as a consequence ∆[∆φ (i) ] = 0 and ∆[∆τ (i) ] = 0. That is, when the fiber length equals an integer number of twist periods the extrinsic contribution is the identity matrix, and so the MCF is conformationally invariant for all fiber shapes with constant radius of curvature. An intuitive explanation for this is to realize that when the condition is satisfied a given core spends as much time on the exterior of the bent MCF as it does on the interior, so the longitudinal compression of the core on the interior compensated the longitudinal dilation on the exterior. To demonstrate the conformational invariance, we cut MCF2 to length L 2 = 345 mm ≈ 42P 2 and MCF0 to L 0 = 344 mm and performed a series of imaging experiments.

Conformationally invariant endoscope
As a first demonstration we show how a focus at the distal end of the MCF is maintained independently of MCF bends. To do so, we use the experimental setup detailed in the Supplemental Information to establish a focus at the output end of MCF2. In Fig. 3(a) we show an image of the focus established for MCF2 held straight (R = ∞). Subsequently, and without changing the injection, we then displaced the distal end to change the radius of curvature of the MCF to R = 0.06 m. The focus, as seen in Fig. 3(b), retains its position to within one spot size. As a base of comparison, Figs. 3(c), 3(d) show that the PSF of MCF0 is translated by a much larger amount for a miniscule conformational change (to R = 0.6 m). A side-by-side comparison is given by Visualization 1 (MCF2) and Visualization 2 (MCF0) which show how the PSF translated when the MCFs undergo the same conformational change. This showcases the dramatic reduction in extrinsic contribution in MCF2 to the level where MCF2 can be considered virtually conformationally invariant.

Discussion
We note that the images presented here [Figs. 4(a),4(b)] are with one-photon contrast. The reason for this is the large native group delay that results from the twist [ Fig. 2(a)]. This native group delay is static and could be compensated with a static pre-compensation so, while not demonstrated here, there are no conceptual difficulties in rendering the setup compatible with ultra-short laser excitation, and to access nonlinear imaging modalities.
We note also that the imaging results presented here were done with forwardcollection of signals, rather than with epi-collection through the MCF. The reason for this is the low fill factor of the cores, which results in poor epicollection efficiency. However, it can be envisioned to increase the epi-collection efficiency of the twisted MCF by adding an external low-index cladding, like in e.g. Ref. [10].
In the above we have shown that the lensless endoscope based on twisted MCF with an integer number of twist periods is conformationally invariant as long as the bend radius of curvature remains constant. In practice there are two main ways that conditions can depart from this ideal condition: (i) Mismatch between length and twist period i.e. L = k · P, k ∈ N ; and (ii) non-constant bend radius of curvature i.e. R varies along the MCF. Concerning (i) we present in the Supplemental Information a map of the extrinsic contribution to phase delay ∆[∆φ (i) ] in the relevant parameters, permitting to identify the parametric sub-space where departure from conformational invariance does not exceed a certain threshold. Concerning (ii) it is not possible to offer insights from analytical considerations, instead an integral must be numerically evaluated, and in the Supplemental Information we calculate ∆[∆φ (i) ] for realistic MCF shapes. We concede that conformational invariance naturally cannot be maintained over the entire parameter space, nevertheless we maintain that the twisted MCF drastically reduces the extrinsic contribution over the majority of fiber conformations likely to be found in practical settings.
Finally we remark that the shorter the twist period of a MCF, the more conformationally invariant it will be. In the limiting case of infinitesimal twist period, the MCF is conformationally invariant over the entire parameter space. In practice, this limit is unattainable because losses due to the curved path of the cores will set in at a certain twist period. In the Supplemental Information we detail numerical calculations which show that, for the MCF parameters used here, twist periods below 2 mm lead to losses above 1 dB/m.
We have demonstrated a lensless endoscope based on a novel, twisted MCF at an operating point where imaging performance is unaffected by the conformation of the fiber. Our findings relax the constraints on fiber shape in lensless endoscopes and pave the way towards imaging tools that reap simultaneously the two principal benefits of using nothing but an optical fiber as imaging probe, namely, its small diameter and its flexibility.

Fabrication of twisted MCF
The MCF was made using the stack and draw method. More precisely, a stack of 487 capillaries drawn from a commercial pure silica tube (Heraeus F300) were first assembled inside a silica tube of ∼50 mm outside diameter. Then 487 rods drawn from a commercial graded-index preform (∆n = 30 · 10 −3 , Prysmian) were inserted into each of these capillaries. This sleeving approach was used in order to guarantee no coupling between adjacent cores of the final fiber by increasing sufficiently the distance between these cores. This stack was drawn into canes of about 5 mm. During this draw, a vacuum was applied between the capillaries and inside them in order to get solid canes free of bubbles that could appear at the different silica interfaces present in the stack. Finally one of these canes was drawn into a fiber of ∼450 µm diameter at low drawing speed (2 m/min), the twist being obtained by rotating the cane at the top of the fiber drawing tower (60 rpm for MCF1 and 250 rpm at MCF2).

Experimental methods
The experimental setup (Fig. 5, Supplemental Information) was designed to measure the MCF properties using two modalities: ultrashort (fs) pulses or continuous wave (CW) in order to perform group delay and imaging measurements, respectively. The first modality was used to perform the phase-stepping spectral interferometry to measure group delays of laser pulses transmitted through different MCF cores (with respect to a reference core, usually the central core), as described in [21]. The second modality is used to measure the PSF and imaging performance during MCF conformational changes.
For phase-stepping spectral interferometry (Fig. 5, Supplemental Information) the laser beam (Amplitude Systèmes t-Pulse, central wavelength 1030 nm, pulse length 170 fs, repetition rate 50 MHz) is expanded with a telescope (lenses L1 and L2) to overfill the clear aperture of a spatial light modulator (SLM, Hamamatsu LCoS-SLM X10468-07). The SLM is used to segment and shape the wavefront into beamlets prior to their injection into the cores of the MCF proximal facet. For group delay measurements only two cores are injected into (the central, reference, core and the core i) resulting in fringes at the MCF output far field. For PSF measurements all cores are injected into in order to produce a focus at the MCF output. The necessary demagnification of the beamlets to fit the core diameters is achieved via a lens L3 and a microscope objective MO1 (Olympus Plan N, 20x NA 0.40). Light, transmitted through the cores, is collected at the fiber distal end using a second microscope objective MO2 (Nikon Plan, 10x NA 0.25) and passes through a linear polarizer LP (Thorlabs LPNIR100) to maximize the contrast of the interference fringes [22]. The image of the distal MCF facet or its near field is monitored with CCD1 (FLIR FL3-U3-32S2M-CS) in order to follow any changes of the transmitted power or the generated interference pattern while the fiber is bent. The MCF output far field is coupled into a multimode fiber (core diameter of 62.5 µm, not shown in Fig. 5, Supplemental Information), linked to an optical spectrum analyzer OSA (Yokogawa AQ-6315A), or imaged entirely onto a camera CCD2 (Thorlabs DCU223M) to record the point spread function (PSF) stability in during MCF bending. Magnification of the objective MO2 and lenses L4, L5 is chosen so that only restricted area of the far field pattern is selected with the MMF probe (fulfilling therefore k x D > 2, where k x is related to the interference fringe spatial frequency and D is diameter of the MMF core collecting the light). The distal end of the MCF is kept fixed on a portable unit that can freely rotate whilst the MCF is bent. The fiber conformation is monitored using a portable camera (not shown). Finally galvanometric scan mirrors conjugated to the SLM active area (telescope L1 and L2) are used to scan the focused distal PSF across the sample for imaging. The galvanometric mirrors impose controllable phase tilts on the input wavefront that are translated to the output distal MCF facet owing to the diagonal MCF transmission matrix.  The experimental setup (Fig. 5) was designed to measure the MCF properties using two modalities: ultrashort (fs) pulses or continuous wave (CW) in order to perform group delay and imaging measurements, respectively. The first modality was used to perform the phase-stepping spectral interferometry to measure group delays of laser pulses transmitted through different MCF cores (with respect to a reference core, usually the central core), as described in [21]. The second modality is used to measure the PSF and imaging performance during MCF conformational changes.

Funding Information
For phase-stepping spectral interferometry (Fig. 5) the laser beam (Amplitude Systèmes t-Pulse, central wavelength 1030 nm, pulse length 170 fs, repetition rate 50 MHz) is expanded with a telescope (lenses L1 and L2) to overfill the clear aperture of a spatial light modulator (SLM, Hamamatsu LCoS-SLM X10468-07). The SLM is used to segment and shape the wavefront into beamlets prior to their injection into the cores of the MCF proximal facet. For group delay measurements only two cores are injected into (the central, reference, core and the core i) resulting in fringes at the MCF output far field. For PSF measurements all cores are injected into in order to produce a focus at the MCF output. The necessary demagnification of the beamlets to fit the core diameters is achieved via a lens L3 and a microscope objective MO1 (Olympus Plan N, 20x NA 0.40). Light, transmitted through the cores, is collected at the fiber distal end using a second microscope objective MO2 (Nikon Plan, 10x NA 0.25) and passes through a linear polarizer LP (Thorlabs LPNIR100) to maximize the contrast of the interference fringes [22]. The image of the distal MCF facet or its near field is monitored with CCD1 (FLIR FL3-U3-32S2M-CS) in order to follow any changes of the transmitted power or the generated interference pattern while the fiber is bent. The MCF output far field is coupled into a multimode fiber (core diameter of 62.5 µm, not shown in Fig. 5), linked to an optical spectrum analyzer OSA (Yokogawa AQ-6315A), or imaged entirely onto a camera CCD2 (Thorlabs DCU223M) to record the point spread function (PSF) stability in during MCF bending. Magnification of the objective MO2 and lenses L4, L5 is chosen so that only restricted area of the far field pattern is selected with the MMF probe (fulfilling therefore k x D > 2, where k x is related to the interference fringe spatial frequency and D is diameter of the MMF core collecting the light). The distal end of the MCF is kept fixed on a portable unit that can freely rotate whilst the MCF is bent. The fiber conformation is monitored using a portable camera (not shown). Finally galvanometric scan mirrors conjugated to the SLM active area (telescope L1 and L2) are used to scan the focused distal PSF across the sample for imaging. The galvanometric mirrors impose controllable phase tilts on the input wavefront that are translated to the output distal MCF facet owing to the diagonal MCF transmission matrix.
6 Bending-induced inter-core group delays for twisted and non-twisted fiber Figure 6 provides additional experimental data for the native and extrinsic contribution to group delay in the cores, presented as a function of the transverse core coordinates. In the case of the non-twisted fiber MCF0 held straight [ Fig. 6(a)] we observe a random distribution of group delays with standard deviation [XXX -VT please specify the standard deviation] that is related to intrinsic imperfections [20,21]. For a relatively small angle (∆α = 25 • ) [ Fig. 6(c)] we observe a linear dependence of the bending-induced group delays on the core position with respect to the bending axis, as detailed in [20]. If MCF0 were employed with ultrashort pulses, increasing ∆α will result in a PSF degradation associated with a non perfect temporal overlap of the pulses arriving at the focal plane of a lensless endoscope. The twisted MCF2, on the other hand, exhibits important native group delays when held straight [ Fig. 6(d)] due to the increased physical lengths of the peripheral cores (Eq. 4). Extrinsic contribution do to conformational change is, however, negligible in this fiber [ Fig. 6(e)] and did not exceed ±20 fs for a MCF geometry bending angle of ∆α = 150 • , which is close to the resolution limit of the used phase-stepping spectral inferferometry technique. Indeed Figs. 3(a),3(b) (main text) showed that the extrinsic contribution to phase delay is less than 2π, this strongly suggests that extrinsic contribution to group delay remains less than an optical cycle or 3.3 fs.

Fabrication of twisted MCF
The MCF was made using the stack and draw method. More precisely, a stack of 487 capillaries drawn from a commercial pure silica tube (Heraeus F300) were first assembled inside a silica tube of ∼50 mm outside diameter. Then 487 rods drawn from a commercial graded-index preform (∆n = 30 · 10 −3 , Prysmian) were inserted into each of these capillaries. This sleeving approach was used in order to guarantee no coupling between adjacent cores of the final fiber by increasing sufficiently the distance between these cores. This stack was drawn into canes of about 5 mm. During this draw, a vacuum was applied between the capillaries and inside them in order to get solid canes free of bubbles that could appear at the different silica interfaces present in the stack. Finally one of these canes was drawn into a fiber of ∼450 µm diameter at low drawing speed (2 m/min), the twist being obtained by rotating the cane at the top of the fiber drawing tower (60 rpm for MCF1 and 250 rpm at MCF2).

Analytical model of twisted MCF
The analytical model is equivalent to the first order of perturbation for which the field distribution is not affected by the twist. Since the refractive index contrast of the fiber is very low, we can consider that the effective index is equal to the group index : with the parameters given in part 2.1A, n ef f = 1.4626 and n g = 1.4909 for the central core at 1 µm wavelength.

Preamble
Using the coordinate system and definitions of Fig. 1 we write the following parametric equations for the coordinates of core i which is at the distance d (i) from the center of the fiber and which makes a ξ (i) angle with the normal vector directed outside the curve described by the fiber: from which we can find the physical length of core i by the line integral

Case: Straight MCF
In the special case of an unbent MCF R = ∞ and the expression reduces to which is how Eq. 2 comes about. The optical path length of core i is We note that we can take n eff because the twist of the MCF is introduced during the fiber drawing process, so no additional stress accompanies the twist. This is different from the case where twist is introduced by mechanical twisting of the MCF in which case the photoelastic effect must be taken into account. The group delay experienced by a pulse of light travelling in core i is then 13) or, relative to the center core with i = 0: = n Which is how Eq. 4 comes about.
And similarly, the phase delay Which is how Eq. 3 comes about.

Bent MCF
In the general case R = ∞, L (i) we must take into account the photoelastic effect of silica. Indeed, the refractive index of the core i is modified as: where p e 0.22 for fused silica is a coefficient depending on Poisson's ratio and components of the photoelastic tensor [23].
If the radius of the curvature is much smaller than the fiber diameter, we can use the first order theory of perturbation and assume that the electric field E is not affected by the bending and remains confined in the core. In the scalar approximation [24]: Since n eff n core in the scalar approximation, the effective index is modified as: The group index is modified in the same manner since [24]: hence As a result, the optical path length of core i of the bent, twisted fiber is We know that d (i) /R 1, so we can get some initial insights from a linear expansion, From which it can be seen that when L = k · P , k ∈ N , OPL (i) does not change with R. This is the basis of the conformationally invariant lensless endoscope as demonstrated in the main text. In practice, however, it is not possible to attain this condition with arbitrary precision. In order to remain within the regime of conformational invariance, the following criterion on the error phase δφ (i) should respect the criterion δφ (i) < π/2, ∀i. Eq. 24 can in principle be evaluated to arbitrary precision as a function of all involved parameters (P , ξ (i) , L (i) , R, d (i) ). Here, we evaluate δφ (i) for an excursion δL in a range around the optimal value of L and bend radius of curvature R, mapped in Fig. 7 for the extreme core with d (i) = 12·16 = 192 µm which experiences the largest error phase of all the cores.

Non-constant bend radius of curvature
All the results obtained with a constant radius can be easily extended by replacing the constant radius R with the radius of curvature of the curve described by the fiber. No general results can be extracted from the analytical formulae and, as an example, we will study the case of a fiber maintained fixed at one extremity and whose other end is subjected to a normal force. It is well-known that, in this case, the curve of the fiber is described by a cubic function y = ax 3 (a depends on the force, the Young modulus and the second moment of area of the cross section of the fiber). Some different bending configuration are depicted in Fig 8 for different value of a and a constant fiber length of 0.8 m. The variation of phase ∆φ (i) (defined in Eq. 3) between the core i and the central core was numerically computed and plotted if d (i) = 160 µm for a nontwisted fiber (Fig 9a) and for a twisted fiber with a helical period of 8 mm (Fig  9b). For the non-twisted MCF, we have already demonstrated that [20]: (1 − p e ) d (i) | cos(ξ (i) )| λ ∆α (25) where ∆α is the angular increase along the curve formed by the fiber. This gives a phase difference of approximately 360π between the core i and the central core if a = 1 m −2 and ξ (i) = 0 while the maximum phase difference in the case of the twisted MCF is only 0.47π. This clearly demonstrates the bending resilience of the twisted MCF.

FEM model of twisted MCF
To determine the properties of the modes with the finite elements method (FEM), we used method describe by Nicolet et al. [25]. This method uses the transformation optics formalism to obtain an equivalent translation invariant problem for which the permittivity and the permeability are anisotropic. The effective indices calculated with this method must be modified according to the mode azimutal number to obtain the real effective indices [26]. Note than this method is not applicable to a bent twisted fiber since the problem can not be reduced to a translation invariant problem in this case. However, this allows to determine mode losses due to the twist which are neglected in the analytical model. Note that the effective index and the group index calculated by this method are defined as: In Fig. 10a, values of effective index and group index obtained with Eq. 9 (blue and red lines) are compared with those obtained with the FEM method (blue and red crosses). Losses of the fundamental mode computed with the FEM method are plotted in orange (right axis). We can verify that the approximations used to obtained Eq. 9 are valid with the parameters used in this article even when losses are very high (> 1dB/m). Note that losses are are lower than 1 dB/km when the twist period equals 8 mm. The evolution of the phase of the electric field with the distance is plotted in Fig. 10b. The linear evolution is due to the beam angle with the MCF axis : according to Eq. 5 when 2πd (i) P 1, the slope is and losses (orange) of a twisted fiber whose parameters are described in section 2 with a helical period P = 2 mm. Losses were computed with the FEM method (Comsol Multiphysics). Blue and red lines were obtained with Eq. 9 while blue and red crosses were obtained with the FEM method.(b) Phase of the electric field computed with the FEM method with d (i) /Λ = 6.