The numerical prediction of the characteristics of directional multimodal couplers for two-photon endoscopy

The first demonstration of 3D imaging using two-photon excitation through a multimode optical fibre took place in 2015 but the presented experimental setup contained a lot of bulk optics which is inconvenient in use. This paper presents a numerical analysis of power distribution transformations in multimode optical fibre couples for the purpose of developing couplers in nonlinear optical systems. The numerical results are proven by experimental results.


Introduction
There are numerous well known nonlinear diagnostic systems used to examine the tissues of the sick. Most frequently single-mode fibres are used in nonlinear systems to deliver an ultrafast signal to the tested sample without any deformations. However, there is an increasing number of works presenting attempts to use multimode fibres in non-linear imaging systems. Multimode fibres usually have a thickness of a few micrometres and thanks to the fact that many modes are propagated in them, they can replace thick laser beams (several mm) composed of many single-mode fibres. On the other hand the big amount of propagating modes generates random speckled patterns at the distal tip of the multimode fibre. Such a speckle field is formed as a result of the interference among all of the propagating modes. To control the propagation of these modes several methods such as wavefront shaping techniques (Č ižmár and Dholakia 2012; Choi et al. 2012;Sivankutty et al. 2016) or a digital phase conjugation (Morales-Delgado et al. 2015a, b, 2016) are used. For example a nonlinear setup using a multimode fibre was first demonstrated in 2015 (Morales-Delgado et al. 2015a, b). In this work a 20-cm section of a multimode fibre was used for pulse signal transfer. Signal fuzziness was compensated by a special system in which the selective phase conjugation method was used (Morales-Delgado et al. 2015a, b), thanks to this modes with similar group velocities are excited. Selectively excited modes are forced to follow the same route inside a fibre and next they constructively interfere on the fibre output to form a spot of minimum time fuzziness. Using digital phase conjugation it is possible to obtain 1-3 lm lateral resolution (Papadopoulos et al. 2013). Another technic which is used to control the propagation of the modes in the multimode fibre relays on using an orthonormal set of illumination pattern (Č ižmár and Dholakia 2012). It was also shown that it is possible to predict the banding angle based on the analysis of the transmission matrix of the fibre and use it for image correction (Plöschner et al. 2015). Recently a very interesting approach was shown where a special light modulator shapes the input excitation wavefront to focus light on the distal tip of the multimode fibre (Caravaca-Aguirre and Piestun 2017). In aforementioned paper the authors demonstrated that proper selection of the multimode fibre is critical for a robust calibration and for high signal-to background ratio performance. The authors obtained the best results for graded-index fibres. In our approach we would like to simplify the construction of an endoscopy setup using fibre couplers made with multimode graded-index fibres. To achieve this goal it is necessary to develop appropriate couplers and a method which would enable forecasting power distribution on the output of the above mentioned couplers. The work presents the numerical analysis of power distribution transformation in multimode couplers and proposes a polished multimode coupler for research on single-and multi-photon fluorescence.

The calculations
In the Fig. 1 the investigated coupler is presented where Port 1 is used to provide excitation light to a sample which is placed on port 4. The fluorescence signal can be measured at port 3 of the coupler. In this simplified numerical approach phase is neglected and only intensity distribution is considered.
The numerical examination of power distribution in the X coupler is based on the method proposed by Frieden (1980). As it has not been widely applied to date (to our knowledge), the Authors will present a short outline.
The intensity distribution within the multimode fibre may be referred to as the joint distribution of a set of random variables. Whatever the light source is, one has to take into account its spatial and angular emission probability density. In the following considerations the energy coefficients assigned to the rays at particular locations and directions will be proportional to the distribution densities of respective random variables. We assume the Gaussian distribution of the light source in the input (port 1).
The transformations of light (ray) distributions in various components of an optical fibre link may be analysed according to Frieden (1980) by means of a functional analysis of multi-dimensional random variables and their functions (Heimrath 1985). The initial ray coordinates are: x; y; zÀspatial; a; b; cÀangular ð1Þ The output ray coordinates are: The input and output joint distribution densities will be referred to, respectively, as: The formal description of ray image transformation is: (the 6-D pseudo-vector functional = of 6-D pseudo-vector arguments) Due to the energy conservation law, the basic relation between input and output densities is: (in the case of unique roots of = only). The general problem here is the prediction of the joint distribution density of six random variables on the output of the optical medium, given the distribution density of input rays and (at least partial) knowledge of = transformation. The usual complication is the lack of unique-roots solutions of =. This becomes particularly obvious when the considerations are limited to simplified distributions in terms of spatial coordinates only, neglecting the angular dependence. The incident rays at various points of the input plane may intersect at one output point, thus the roots of = at this point are not unique. If the number of the roots is r [ 1, then the right side of Eq. (5) must be a sum of r input densities contributing to the output intensity at this point. If the coordinates of the i-th contributing input ray are x i , …, c i , then finally: where J x i ; . . .; c i =x 0 ; . . .; c 0 ð Þ is the Jacobian of transformation = between the respective pairs of input/output points: By the law of large numbers the above densities of random variables closely approximate the corresponding densities of ray strikes.
The major setback for the analytical solution of Eq. (6) is that the analytical form of = should be known, so that all partial derivatives encountered in Eq. (7) might be calculated. In many cases of practical interest the explicit functional form of = cannot be specified. The partial derivatives should be approximated in some other way. The proposed solution is to numerically trace the ray (Frieden 1980): and its differential counterpart: Such a ray-trace produces two sets of output ray coordinates: Thus, one obtains the numerical values of respective differentials. For sufficiently small values of dx i , …, dc i the partial derivatives of Eq. (7) may now be approximated. Each such pair of numerically traced rays yields one point on the output spot intensity curve which is the merit of the whole procedure.
One strong requirement appears when choosing the ray racing method: it must be accurate to a very high degree. This is an indispensable condition-the whole procedure must be sensitive to minute variations between the ray and its differential counterpart. From among a number of practical procedures for numerical ray tracing, the Authors chose two methods: one proposed by Sharma et al. (1982) and the one proposed by Heimrath (1985). Both methods yield similar results, however, Sharma's method based on the Runge-Kutta integration of differential equations offered faster convergence and the results presented below were obtained by this procedure.
The distribution of the refractive index for a graded-index fibre was approximated according to the following expression (Ankiewicz and Pask 1977) The initial illuminance conditions are highly dependent on the measurement setup configuration. Since the Authors used the He-Ne laser, the initial light distribution was assumed to have the Gaussian density function and to be highly directional. The density of ray distribution in the main fibre was calculated: both the refraction upon incidence on input plane and the radiation of non-guided rays were taken into account following the results presented by Ankiewicz and Pask (1977). The rays and their differential counterparts were scanned over the input plane. They were distributed uniformly, forming an 30 9 30 array (900 rays). A value of intensity is assigned to each ray according to Gaussian distribution. Thus in this simplified model we can think about rays as pixels which are transported through the coupler and create intensity pattern at the port 4 of the coupler. In the Fig. 2a the initial distribution of light lunched at port 1 and in the Fig. 2b obtained after propagating through the coupler are presented.

The measurements
The Authors examined the X-couplers manufactured at the Wroclaw University of Science and Technology, Poland. The couplers were made using the side polishing method to facilitate the evaluation of the distribution of power propagation. The main components of the couplers were the multimode GI optical fibres (50/125 core/clad diameters). During experiments and numerical analysis we considered 0.5 m long fibres. The fibre branches were fitted into the supporting glass cubes with a rift of radius curvature R and then polished to depth h (see Fig. 1). Both similar components of the coupler were then mounted together with an immersion liquid of refraction index n i , and thickness d i . The distance between the centres of cores depends on both polishing depths and immersion film thickness and is referred to as d. The relative longitudinal shift of the centres of the polished zones in both branches is referred to as Dz. This process allows for relatively easy and accurate control of the dividing factor. The optical setup for measuring the dividing factors is presented in Fig. 3. The light source was the He-Ne laser 630 nm. The light beam was formed to fit according to the numerical aperture of the fibre. The beam was focused on the centre of the input fibre. In order to obtain Gaussian distribution at the port 1 of the coupler we use a long (500 m) GI multimode fibre.
The output beam was detected and measured by means of photodiodes in a homodyne setup.

Results
The Authors calculated the coupling coefficients or normalized power P for a number of possible configurations. In all cases the radius R of fibre curvature was fixed at 125 mm. One of the most interesting results is the prediction of the normalized power and its dependence on such technical parameters as immersion liquid index n i and polishing depth h.  The experiment was to measure and calculate the influence of longitudinal shift Dz on the normalized power P 4 /(P 2 ? P 4 ) coupled to the coupling fibre. The immersion index was fixed with the value n i = 1.468. The results (the numerical simulation curves vs. the measured points) are presented in Fig. 4, for two different polishing depths (h 1 = 25 lm and h 2 = 17 lm). The plots show quite good accordance of measured and calculated results.

Conclusions
The presented evidence leads to the conclusion that the proposed method of numerical simulation enables a relatively easy and inexpensive (once the computer program is executable) way of predicting the optical characteristics of multimode, gradient-index of optical components. The analysis was limited to X couplers only, but the method may be tested against other devices as well. These results are promising toward the implantation of such methods in predicting light distribution on the ports of the multimode GI coupler. In the next manuscript we would like to consider also phase of the lunched signals in our numerical model and predict the speckle pattern on the port 4 of the investigated coupler. Using a properly designed coupler can significantly simplify a two photon endoscopy setup.