Field of the invention This invention relates to improved methods of making use of multimode optical fibres, and to instruments and systems using such fibres. The invention may be used, for example, in imaging instruments such as endoscopes and in communication systems.

Background to the invention

In a similar fashion to diffusers or other highly scattering media, multimode fibres (MMF) deliver coherent light signals in the form of apparently random speckled patterns. In contrast to other random environments in which light propagates along complex paths, MMF features remarkably faithful cylindrical symmetry.

The theoretical description of light transport processes within ideal MMFs has been developed for over half a century [1-4]. This elaborate theoretical model is, however, frequently considered inadequate to describe real life MMFs that are manufactured by drawing melted silica preforms. Such fibres are commonly seen as unreliable and the inherent randomisation of light propagating through them is typically attributed to undetectable deviations from the ideal fibre structure. It is a commonly held belief that this additional chaos is unpredictable and that its influence grows with the length of the fibre. Despite this, light transport through MMFs remains deterministic. The prospect of deterministic light propagation within MMF has only recently been utilised through methods of digital holography and by adopting the concept of empirical measurement of the transformation matrix (TM) [5-1 1 ]. This technique, developed on studies of light propagation through highly turbid media [12-17], has opened a new window of opportunities for MMFs to become extremely narrow and minimally invasive endoscopes, allowing sub-micrometre resolution imaging in deep regions of sensitive tissues [9,18]. However attractive, this technology suffers from several major limitations, the most critical being the lack of flexible operation: any bending or looping of the fibre results in changes to its TM, thus rendering imaging heavily impaired. All current methods exploiting MMFs for imaging require open optical access to the distal end of the fibre during the time consuming measurement of the TM. Furthermore, this characterisation must be repeated upfront for every intended configuration (deformation) and any axial distance of focal plane behind the fibre before the system can be used for imaging [7,19]. The necessity to determine the TM empirically therefore constitutes a major bottleneck of the technology and it would be immensely advantageous to obtain the TM by another route, ideally on the basis of numerical modelling.

Our experimental studies challenge the commonly held notion that classifies multimode fibres as unpredictable optical systems. Instead, we demonstrate that commercially available M Fs are capable of performing as extremely precise optical components. We show that, with a sufficiently accurate theoretical model, light propagation within straight or even significantly deformed segments of MMFs may be predicted up to distances in excess of hundreds of millimetres. Harnessing this newly discovered predictability in imaging, we demonstrate the unparalleled power of M F-based endoscopes, which offer exceptional performance both in terms of resolution and instrument footprint. These results thus pave the way for numerous exciting applications, including high-quality imaging deep inside motile organisms.

We asked whether such modelling is feasible in the supposedly chaotic environment of multimode fibres that however, unlike other random media, feature a remarkable cylindrical symmetry. Indeed, numerous studies have already indicated that at least some aspects of the ordered behaviour (propagation constants of modes) can 'survive' over very large distances [20-23].

Our study represents a major advancement on these efforts. With the availability of experimentally measured TM stored in the memory of a computer, all light transport processes can be emulated with great accuracy and the optical fibre can be subjected to detailed investigation [24]. While working with relatively large numbers of modes («500) we have initially used the shortest practical fibre segmental 0mm) in order to minimise the complexity of the problem. This important intermediate step allowed us to eliminate influences caused by imperfections (unavoidable misalignment) in our experimental settings and importantly, establish parameters of the fibre (core diameter, numerical aperture) with sufficiently high precision.

Only when corrections for these aspects were implemented could we see a perfect agreement with the theoretical model used. This not only confirmed the existence of propagation-invariant modes within the fibre but also very accurately matched their output phases, vitally important for imaging applications. Progressing to =100mm-1ong fibre segments we have faced the first challenges related to polarisation coupling and significant deviations of refractive index from the ideal step-index profile, which required major enhancement of our numerical modelling and experimental methods.

With such newly developed ability to predict light propagation even through such distances within the fibre, we could for the first time rigorously investigate the influence of fibre deformation (bending) on the resulting TM. The study has shown that even significantly bent fibres are perfectly predictable, thus allowing computation of their TM based purely on the observation of the fibre geometry. Finally, all of these new aspects were brought together to test the performance of imaging in such an enhanced system. We have shown that imaging can be achieved without experimental TM acquisition, in both straight and deformed fibres at an arbitrary distance behind the distal fibre facet.

Summary of the invention

Accordingly, the present invention provides a method of designing a light transmission system using a multimode optical fibre, in which the light transmission properties of the fibre are predicted by numerically modelling the transmission matrix of a given fibre.

From another aspect, the invention provides a method of making a light transmission system, comprising providing a multimode optical fibre, predicting the light transmission properties of the fibre by the foregoing method, and coupling the fibre to a light source and to a light detector in a manner determined by said prediction.

The light transmission system may be incorporated in an endoscope, or a signal transmission system. Preferred features of the invention will be apparent from the description and claims. Brief description of the drawings

Embodiments of the invention will now be described, by way of example, with reference to the drawings, in which: Figure 1 illustrates the analysis of a short segment of fibre, a), Organisation of input and output modes, b), experimentally measured TM. c), theoretically predicted LP PIMs. d), Conversion matrix between the representation of FPs and LP PIMs. e) and f), converted TM before and after the optimisation procedure, respectively.

Figure 2 illustrates polarization coupling effects in MMF. a-c), Data for the 10mm long fibre, a, TM of LP PIMs. b), Poincare sphere depicting polarisation change of LP PIMs. c), Projections of LP PIMs polarisation state on the mode pyramid, d-f), Equivalent of (a-c) for the 100mm long fiber, g-i), Data for the 100mm long fibre using CP PIMs. j), Experimentally measured phase difference between CP PIMs having opposite spin (the influence of SOI), k), Numerically simulated equivalent of (j). I, CP PIMs output amplitude, m), Experimental synthesis of input PIMs and their corresponding output. Figure 3 shows optical phases of PIMs.

a-e), Data for the 10mm long fibre, a), Experimentally measured phases, b), and c), Phases predicted by numerical model, unwrapped and 2-rr-wrapped respectively, d), Difference between phases experimentally measured and theoretically predicted, e), Phase agreement.

f-n), Data for the 100mm long fibre, f), g) and h), Assumed profile of refractive index, phase agreement and difference between experimentally obtained and theoretically predicted phases of PIMs respectively for the case of ideal step-index fibre, i), j) and k), corresponds to a model with included dopant diffusion and I), m) and n), represents correction for diffusion as well as fine index modulation across the fibre core.

o), p) and q), Equivalents of l),m) and n) for 300mm long fibre.

Figure 4. illustrates the influence of fibre deformation, a), Arrangements of deformed fibre used in the experiment with their Roman ID number, b) and c), Experimentally measured and theoretically predicted DO corresponding to deformation (V). d), Empiric estimation of scaling factor ξ. e) and f), Experimentally measured and theoretically predicted influence of deformation on PIMs.

Figure 5 shows application to imaging in the form of imaging of USAF 1951 resolution target performed for three lengths of fibre: 0mm (a-d), 100mm (e-h) and 300mm (i-l). In each case we show imaging with experimentally acquired transformation matrix (a, e, i) and three theoretically predicted transformation matrices corresponding to ideal step-index fibre (b, f, j), profile with correction for dopant diffusion (c, g, k) and profile with full correction (d, h, I).

Figure 6. Imaging with deformed fibre, a), Imaging with empirical TM of straight fibre, b), Imaging with TM after full DO correction, c), TM with only diagonal components of DO applied, d), TM with only phases of diagonal components of DO applied.

Figures S1 to S16 are supplementary figures referred to in the text.

Description of preferred embodiments

This description consists of a general description explaining embodiments of the invention, and an Addendum giving additional detail on calculations and methods.

Where reference is made to Supplementary Information and to Online Methods, these can be found at http://www.nature.com/nphoton/journal/v9/n8/full/nphoton.201 5.112.html. The Supplementary Information includes moving images which cannot be reproduced herein.

General description

Identification of Propagation Invariant Modes

Our experimental geometry, introduced in Online Methods, allows measurement of the TM in the representation of diffraction limited focal points (FP) sometimes also called "delta peaks", which from the experimental perspective represents the most convenient choice. The input FPs (at the proximal end of the fibre) can be generated using a spatial light modulator that simply steers a focused laser beam into a required location at the input fibre facet. The output (distal) facet is imaged on a CCD chip and FPs are acquired from values of individual pixels.

The chosen set of FPs is arranged across an orthogonal grid as shown in Fig. 1 a and ordered as indicated by the red line. The experimentally measured transformation matrix M for a 10 mm long fibre segment is shown in Fig. 1 b. Each row of M represents amplitudes and phases of all output FPs for a single input FP mode sent into the fibre. Due to space constraints, here, as well as in Figs. 1d-f, the basis of modes for the TMs shown has been reduced to 1/3 of the full dimension. The complete transformation matrices are presented in Supplementary Figures S5-S7. Following the TM acquisition we can numerically emulate the optical fibre as an optical system and predict the outcome of any optical field being sent into the fibre and thus validate the correctness of any theoretical prediction.

The simplest form of theoretical description (scalar and paraxial approach reviewed in Supplementary Methods S1) predicts the existence of linearly polarised (LP) modes that do not change their field distribution during propagation through the fibre. A series of such propagation invariant modes (PI Ms), also known as eigenvectors of propagation operators or eigenmodes, are shown in Fig. 1c. Again due to space constraints we only show PlMs of a fibre having many fewer modes. PlMs are defined by a pair of indices m and I. Index I refers to orbital angular momentum of a given mode with magnitude equal to Ih/photon [25].

Whether such modes would remain unchanged after travelling through the real fibre can now be tested using the experimentally measured TM: Each of these theoretically predicted PlMs can be constructed as a superposition of input FPs. Such a vector can be then virtually sent into the fibre which is implemented by its matrix multiplication with the TM. The resulting output vector of FPs should contain the identical PIM, differing only by a phase constant. Carrying out such operation for all modes simultaneously is mathematically equivalent to converting the experimentally measured M into the representation of PlMs. This is achieved by constructing a conversion matrix T shown in Fig. 1 d in which each line represents a single theoretically predicted PIM expressed as a superposition of input FPs. PlMs in T are ordered as indicated by the white line in Fig. 1 c. If the theoretically predicted PlMs are the true eigenvectors of the experimentally measured transformation matrix M, its conversion into the representation of PlMs should result in purely diagonal matrix, indicating that each input mode was perfectly conserved.

The converted transformation matrix Mo = TMT† is shown in Fig. 1e. Apparently Mo is not diagonal, which might lead to the conclusion that the optical fibre does not follow the theoretical model. However, off-diagonal components can also appear due to even a very small misalignment of the fibre. The misalignment space of the degrees of freedom is very large: 3-D position, two tilts and one defocus, each on both sides of the fibre. Moreover, these 12 degrees of freedom are intrinsically intertwined with uncertainty in the radius of the fibre core (a) and the value of numerical aperture (NA). We have developed an optimisation procedure described in Online Methods (source is available in online repository [26]) that simultaneously corrects the TM for alignment imperfections and adjusts the values of a and NA. The optimised result finai, shown in Fig. 1f, carries 93% of optical power on the main diagonal showing an excellent match between the scalar theoretical prediction and the corrected experimental data.

Polarisation coupling of modes

Each mode can be defined in two orthogonal polarisation states and only when both are taken into the consideration can the TM be considered complete. Online Methods explains how we can control the polarisation in our geometry in order to take such complete measurements of TM. After the optimisation and conversion into PIMs, such a complete TM will now have four quadrants, those containing the main diagonal indicate that the polarisation of a given mode was conserved, while the remaining ones indicate mutual coupling between polarisation states. The optimised TM with input PIMs defined by two orthogonal linear polarisation states for the 10mm long fibre is presented in Fig. 2a (see Supplementary Figures S7-S10 for complete data sets). Coupling between polarisation states is clearly present but relatively weak. The change of polarisation state of individual LP PIMs can be efficiently visualised on the Poincare sphere as shown in Fig. 2b. It is seen that the polarisations of all modes remain linear but their orientation is rotated by up to 45° (which corresponds to a shift of 90° along the equator on the Poincare sphere). The same data are visualised by placing the polarisation states of output PIMs into the shape of the I -m pyramid as defined in Fig. 1 c. The colour defining the polarisation state corresponds to the colour of the Poincare sphere in Fig. 2b.

Equivalent measurements were taken for the 100mm long fibre (see Figs. 2d-f). Once again we see output polarisation remaining linear (apart from 'misbehaving' modes of I = ±1 discussed later), but this time the rotation of polarisation is much stronger. Although this effect remains highly ordered, the observed polarisation change proves that the LP PIMs can no longer be considered truly propagation-invariant. To account for this effect, we have enhanced the theoretical description to the fully vectorial model (presented in Supplementary Methods S1 .1 ). This advanced description showed that only circularly polarised (CP) PIMs remain unaffected by propagation through the fibre. Coincidentally, CP PIMs have almost identical distribution of the field amplitude as LP PIMs which considerably simplifies their modelling. The experimental study using circularly polarised modes propagating through a 100mm long segment of fibre is summarised in Figs. 2g-i. As expected, with these CP PIMs we no longer see transitions far from the original polarisation state. The CP PIMs thus indeed represent propagation-invariant modes as predicted by our vectorial theoretical model. Interestingly, there is a broken symmetry between the modes having the same I and m indices but opposite circular polarisation state(spin). Their optical fields are no longer symmetric as the wavefront helicity and spin are orientated the same way in one case but opposite in the other. The vectorial description predicts that the degeneracy in such pairs of modes is removed and in consequence they do not propagate with the same phase velocity [4,27]. This effect is well known in quantum optics as spin-orbit interaction (SOI). Such behaviour can be very precisely identified in our experimental data (Fig. 2j) and also predicted by our theoretical model (Fig. 2k).

Figure 2I shows amplitudes of the output PIMs (diagonal components of TM in Fig. 2g). The most dominant losses are seen in modes of I = ±1 which, as explained in Supplementary Methods S1 .1 , is caused by highly exceptional behaviour of these modes, for which the CP PIMs-based model is not valid. To conclusively confirm the correctness of these findings we tested our results in our experimental geometry, took into account all the detected alignment imperfections and experimentally synthesised the full orthogonal set of CP PIMs (see Online Methods for more details).

The PIMs have been recorded in both circular polarisation states prior to entering the fibre as well as after leaving the fibre. One example showing a superposition of two such modes with the same m but opposite I indices is shown in Fig. 2m. The full sets of PIMs for the 10mm and the 100mm long fibre segments are presented in Supplementary Media SM2- SM5).

Examination of Optical Phases

Evaluation of the ability to predict the optical phases of PIMs is of paramount importance to this study and conclusive in answering whether theoretically predicted TMs could be used for imaging. Phases of CP PIMS measured in the 10mm long fibre segment are shown in Fig. 3a. Since the phases are wrapped within an interval of h-π,ττ], apart from mirror symmetry, (visible due to negligible contribution of SOI at this length) no ordered behaviour is immediately apparent. The simulation for phases of PIMs using our numerical model is shown first unwrapped in Fig. 3b and also wrapped in Fig. 3c. The difference between the simulation and the experiment is shown in Fig. 3d. This very good agreement only occurs if the length of the fibre is known with a very high precision (in the order of units of pm). We have determined the fibre length by seeking the highest value of the quantity referred to as total phase agreement (PA):

PA =^∑f ₌₁ e ^iMj ] ² where N is the number of PI Ms and Δφ j is the phase difference between the experimentally measured and theoretically predicted phases of j-th PIM. This quantity equals 1 for perfect agreement. PA as a function of assumed fibre length is plotted in Fig. 3e. Its peak value of ~ 83% is very satisfactory given that the propagation constants of PIMs have been matched with a relative accuracy better than 10 ^~5 Moreover, it is clearly seen that the deviations are not random (noisy) but rather form a smooth surface (Fig. 3d).

Repeating the same experiment with the 100mm long fibre results in the PA peaking around the expected value of fibre length (Fig. 3g), but the PA peak value has decreased dramatically. The difference between the experimental and the theoretical phases is shown in Fig. 3h, clearly demonstrating that the aberration seen in Fig. 3d has grown proportionally with fibre length. This effect could not be explained by any aberration of our optical system. We therefore conjectured that the observed phase differences could be caused by deviations of the refractive index from the ideal step-index profile shown in Fig. 3f, originating e.g. from the diffusion of dopant material between core and cladding [28,29] (Fig. 3i).

To quantify this deviation, we employed perturbation calculus explained in Online Methods. The diffusion is commonly quantified by diffusion length d which became an additional free parameter in our optimisation procedure (see Online Methods). With this new model taking diffusion into account the PA was substantially enhanced (Fig. 3j) but still could not fully explain the observed deviation demonstrated in Fig. 3k. Therefore, we searched for further deviations of refractive index within the fibre core and using perturbation calculus we have arrived at their description using zero-order Zernike polynomials (Z° ₄ (r/a), Z°e (r/a), Z°e (r/a), ...) finding that only the first two of this series are required to entirely explain the observed phase deviations.

Even with such minuscule modifications in the refractive index profile (shown in Fig. 31) the PA has been enhanced beyond the value of 0.95 (Fig. 3m). The difference between the experimental and the theoretical phases of PIMs is shown in Fig. 3n. Further extending the length of the fibre to 300mm we have found that for satisfactory agreement shown in Figs. 3o-q we had to supplement one more free parameter of fine modulation. Transformation matrix for deformed fibres

Our experimental geometry was designed to allow bending of optical fibre with curvatures reaching the long-term damage threshold given by the manufacturer. We have experimentally measured TMs for several different fibre geometries by translating and rotating the output end of the fibre which was monitored with high precision. Using this information we reconstructed the shape of the fibre numerically by minimising its elastic energy. These shapes are represented in Fig. 4a by solid lines and complemented by the actual layout of the fibre (marked by black dots) obtained from camera images. The changes introduced to TM by such deformations can be best characterised by a deformation operator (DO,), which describes the transition of the straight fibre TM (Ms) to the TM of the bent fibre (Mb):

(Mb) = ₅

An example of an experimentally measured DO (in the representation of PIMs) is shown in Fig. 4b (only one polarisation component, for full DO see Supplementary Fig. S12), which corresponds to the largest curvature (V) we introduced to the fibre. It is clearly seen that even for a curvature reaching the damage threshold, DO remains highly diagonal thus indicating that majority of PIMs were very well conserved. The main effect of deformation is therefore only a small relative phase shift of the PIMs. Diagonal components of DO (influence on PIMs) for all fibre deformations are presented in Fig. 4e. Apparently, PIMs with low I and m indices are the most vulnerable to such deformations.

Our theoretical description of bent fibres is detailed in Online Methods. This analytical model can directly handle only fibre segments with uniform curvature. Therefore, to model our experimental conditions as authentically as possible, we have constructed the corresponding theoretical TM from 300 segments, each with constant curvature corresponding to the minimal elastic energy fibre layout (Fig. 4a). The results of the simulation indicated effects of the same kind as observed experimentally, but noticeably stronger. Therefore we searched numerically for an empirical scaling factor with which we had to reduce the curvature of the fibre in our model in order to obtain the optimal agreement with the experiment. For all our fibre deformations (ll)-(V) the scaling factor found was very similar (see Fig. 4d) and its value was estimated to be 0.77±0.02. Further theoretical investigation revealed that this effect could be caused by linear variation of refractive index with respect to material density that changes under the influence of deformation-induced stress.

This leads to exactly the same effect as fibre bending, but with opposite action. Assuming that (n-1 ) is directly proportional to the material density [30] the scaling factor can be expressed as ξ = 1- ~ (1 -2a) where σ is Poisson's ratio (see Online Methods). With a value of σ = 0.17, which is a typical material constant of fused silica given by many manufacturers, we reached a value of ξ = 0.79 which is in perfect agreement with our experiment. Theoretically predicted DO, corresponding to fibre deformation (V) (with the empirical scaling factor applied) is shown in Fig. 4c (see Supplementary Fig. S13 for all polarisation components) and diagonal components for theoretically predicted DOs are depicted in Fig. 4f.

Importantly, the theoretical model predicts that the DO does not change if the plane in which the deformed fibre lies is oriented differently. This was confirmed by extending this study into the third dimension, which is presented in Supplementary Results SR1 . Application to imaging

Being able to predict the correct spatial distribution of the optical fields of PIMs together with their correct phases after propagating through the fibre allows us to construct a transformation matrix of a multimode fibre in any representation of modes, purely based on the grounds of theoretical modelling. For imaging purposes the representation of FPs is the natural choice. In this representation, each column of the TM tells us what superposition of input FPs we need to provide in order to get a single diffraction limited FP at the distal end of the fibre. In consequence we can design a corresponding holographic modulation for our SLM (taking into account all detected misalignment imperfections) in order to synthesise such output FP experimentally. Using various forms of theoretically predicted TMs we have accordingly synthesised a set of FPs across an orthogonal grid of 75x75 different positions, which have been sequentially lit-up by the SLM to efficiently raster-scan the object during image acquisition. Each pixel value of the resultant image was obtained by measuring the total optical power transmitted through the object while exposed to a single FP. More details on this approach, as well as the demonstration of proximal end imaging, can be found in the Supplementary results SR2 and the Supplementary Media SM7. For our demonstration of imaging we have chosen a negative USAF 1951 resolution target as the imaging object, which was placed in a close proximity to the distal end of the fibre. In Supplementary Results SR3 we show that imaging is not restricted to the proximity of the fibre facet but with the use of a free-space propagation operator the image plane can be easily shifted to an arbitrary distance behind the fibre facet. Our evaluation starts with the 10mm long fibre. Fig. 5a shows traditional imaging with experimentally measured TM [7-9,18]. The imaging performance is compared with those obtained while using three different theoretically predicted TMs that correspond to the ideal step-index profile (Fig. 5b), profile with diffusion of dopant material (Fig. 5c) and finally the fully corrected profile (Fig. 5d).

We see that the imaging was successful in all cases, and that a small aberration due to refractive index profile imperfections (Fig. 3d) has minimal influence on imaging quality. Equivalent tests were performed for 100mm-long fibre (Fig. 5e-h) and 300mm-long fibre (Fig. 5i-l). Employing the theoretically predicted TM of ideal step-index fibre at these length scales entirely fails to produce a recognisable image, which clearly demonstrates the importance of the phase corrections studied earlier.

Much better but still heavily impaired imaging can be obtained after taking diffusion effect in refractive index profile into account (Fig. 5g,k). But only after implementing additional correction for the fine modulation of refractive index (Fig. 5h,l) does the imaging quality approaches that for the case of an experimentally measured TM (Fig. 5e,i). The influence of bending is demonstrated in Fig. 6a. where the TM measured on a straight segment of fibre was used for imaging after the fibre was bent as illustrated by deformation (V) in Fig. 4a. Results obtained after applying appropriate theoretically predicted DO (Fig. 4c) are shown in Fig. 6b. Further we have tested DO with all off -diagonal components set to zero (Fig. 6c) and additionally also set absolute values of diagonal components to unity (Fig. 6d). Both of these simplifications exhibit negligible effects on the imaging performance, thus confirming that DOs are highly diagonal and therefore almost perfectly commuting with one another. This will immensely simplify predicting DOs of more complex deformations since the order in which individual bends follow along the fibre won't play an important role.

Discussion

We have demonstrated that MMFs are highly predictable optical systems. We have developed new methods to measure and analyse TMs of real-life MMFs as well as a complex theoretical framework to describe the light-transport processes within, taking into account several important deviations of the refractive index profile and fibre curvature. We have not identified any significant randomising processes in 300mm Song fibre segments, and based on the behaviour observed it is very unlikely that highly ordered light propagation would be obstructed in the length scale of even several metres.

From our results it is clear which modes are the most vulnerable to fibre deformations and in turn which are the most immune. This is of crucial importance not only for endoscopy but also for many other disciplines including telecommunications (mode-division multiplexing).

We have found that the strength of deformation-induced TM changes can be tuned by selection of refractive index and Poisson's ratio. This may bring an exciting new possibility to design fibres featuring suppressed bending effects (minimising value of ξ), however the ultimate goal of designing fibres with ξ = 0, for which the light transport processes are entirely invariant to fibre deformation would require either negative refractive index or negative Poisson ratio, both very difficult to achieve experimentally.

We speculate that the latter could possibly be approached by providing a highly anisotropic sleeve around the fibre. The startling agreement of theoretical predictions with experimental reality allowed us to perform imaging using theoretically predicted TMs in straight and curved fibre segments sufficiently long for many practical applications in life sciences, medicine and other disciplines. These new possibilities bring a step-change for this technology by eliminating the most critical drawbacks related to the requirement for empirically measured TMs as well as lack of flexibility. We believe there are currently no fundamental problems standing in a way of further fast development of this exciting technology, and with sufficiently sophisticated engineering we might soon see a multitude of exciting applications such as observations of single or multiple neurons in unrestrained and awake animal models.

Addendum

In this section, we give further detail on methods and experimental setups we have used in developing this invention, followed by a section on supplementary methods.

Methods

Calculating the changes of propagation constants due to deviations in refractive index profile Light transport processes in multimode fibres are theoretically well understood and the necessary theoretical considerations are briefly reviewed in Supplementary Methods S1. In order to find corrections to the propagation constants arising from slight deviations of the refractive index profile we used perturbation calculus with respect to the basic scalar theory. These corrections are then added to propagation constants calculated by the weak guidance approximation of the full vector mode! to obtain corrected propagation constants.

We denote the refractive index in the idea! step fibre by n(r) and the slightly modified index by n 0 (r). Starting from the scalar wave equation ΛΨ-(η/ο) ² St ² Ψ = 0 and separating the time and z-coordinate along the fibre as Ψ(Γ,φ,ζ,ί) = ψ(Γ,φ)βχρ[ί(βζ-ωί)]), we arrive at the He!mholtz equation for the transversal part of the wave ψ(Γ,φ) in the step refractive index n(r):

[A± +ko ² n ² -β ² 1ψ(Γ,φ) = 0, (3)

where A± is the transversal part of the Laplacian. The solutions can be indexed by the angular index 1 (angular momentum, or topological charge) and radial index m as mentioned before. A similar equation holds for the perturbed mode functions ψ 0 (Γ,φ) ,

but with n replaced by n' and β replaced by the perturbed propagation constants β'. Next we write the square of the modified refractive index as n' ² (r) = n ² (r)+g(r), where g(r) is a small perturbation, express the mode functions ip'im as superpositions of the unperturbed mode functions ψ im, and perform the standard perturbation calculation. In this way we arrive at the following first-order correction to the propagation constant of the mode

ψΊπι:

k0 ²

ΔβΐΓΠ =— <ψ im |g|lJJ lm >

For a given modified index n 0 (r), these corrections can be readily calculated by numerical integration. The polarisation index σ of the mode does not influence the correction. Calculating the changes to transformation matrix due to fibre bending

Eigenmodes of a bent fibre are different from those of a straight fibre. Taking into account the fact that in practice the radius of curvature is larger than the core radius by several orders of magnitude, one could be tempted to think that the two sets of modes will differ only slightly. This is not true, however, due to the very small index difference between the core and the cladding. In fact, the bending introduces effective index changes that are comparable to this index difference, and therefore perturbation theory would yield inaccurate results. Fortunately, one can still describe modes of the bent fibre approximately using a calculation that reminds of perturbation theory, as we show in the following.

Consider a short element of the fibre centred at the origin of Cartesian coordinates with the fibre axis oriented along the z axis. Suppose now that the fibre is bent in the xz plane with the radius of curvature p » a with the centre of curvature at the point (ρ,Ο,Ο). The local longitudinal wavenumber (along z axis) will no longer be constant across the fibre cross section but will depend on x as k _z (χ)=β' ² (1 +x/p), where β' is the value of the longitudinal wavenumber on the axis— the propagation constant. This yields the following equation for the mode in the xy plane:

Next we express ψ' as a superposition of the modes of the straight fibre, ψ' =∑i c i ψ i (for simplicity we have replaced the pair of indices l,m by a single index i), and substitute into Eq. (5). Using also Eq. (3), we get

Multiplying this equation by, integrating over the xy plane, using the orthogonality of arranging the terms slightly and interchanging the indices i and j, we get

1

Here we have defined the matrix element of the x coordinate Ai _j ≡ <≠\xM) = ^{!o rdr !o d} *^ ^rc →J ₍₈₎ τ, φ I

12

where ¹ _2π ] ^1/2 is the normalisation constant of the i th mode. Equation (7) is in fact an rdr _Q d<p

eigenvalue problem for 1/β' ² The eigenvectors (sequences of coefficients c , ) then determine the corresponding modes of the bent fibre. This equation can be further simplified by taking advantage of the fact that the propagation constants β, are limited by the condition n cladding ko < pi < Π core ko and therefore are all very close to n core k o. Making the substitutions β i = n core ko +Δβί and β' = n core ko +Δβ' , inserting into Eq. (7), neglecting the term containing the product A Δβ , and returning back from Δβ, ,Δβ' to βί ,β',

we finally get the equation iC _i - ^f ∑ _j A _ij C _j = pc _i (9)

This equation shows that the propagation constants are eigenvalues of the matrix B with entries

B _tj = - The advantage of Eq. (9) compared to Eq. (7) is that its

eigenvalues are directly the new propagation constants. Therefore the evolution operator of the state along the fibre segment of length L is simply e iBL. Evolution along a fibre that is bent non-uniformly can be expressed as a product of such operators corresponding to sufficiently short fibre segments in which the curvature can be considered constant.

In the previous we have assumed that the refractive index profile of the deformed fibre remains identical to the original straight fibre. The fibre deformation, however, causes local density changes and consequently also refractive index changes, which influences light propagation. When the fibre is bent, its outer side becomes longer and the inner side becomes shorter. Such longitudinal changes of the length of infinitesimal fibre elements cause corresponding lateral changes of their width of opposite sign and smaller in magnitude by the factor of the Poisson ratio σ. In other words, the diagonal elements of the deformation tensor are related as εχχ = £yy = -σε∑ζ. The relative changes of the element volume and density are then Ττέ = ε _ζζ (1 - 2σ) and -Ττέ = -ε _ζζ (1 - 2σ), respectively, provided that the deformation is small. If we assume that the refractive index has the property that n-1 is proportional to the density, which seems to be reasonable for fused silica the fibre is made from, we obtain a modified refractive index n' = η-(η-1)ε zz (1 -2σ). Inserting this modified index n' instead of n into Eq. (5) together with ε _zz = -x/p, we find after neglecting one insignificant term:

[Δ + kW - β' ² (l + [l - (1 - 2σ) =^])] ψ' =0 We see that the resulting equation differs from Eq. (5) by the factor ξ = 1-(1-2σ)(η-1)/η standing in front of the curvature 1/p. This means that the fibre behaves equally as if the index did not change but the bending were weaker by the factor of ξ. This also means that the effect of index change alone (due to fibre deformation) is of a similar character as the effect of the bending itself, just of opposite sign.

In general, the transverse deformation leads to a change of the shape of the cross section of the fibre core, which could also affect the modes. However, it is not hard to show that for the particular deformation in question, i.e., ε xx — ε yy = -σε zz = σχ/ρ, the cross section remains circular (still with radius a) up to the first order in a/p; the deviation from the circular shape is of the second order in a/p and therefore completely negligible.

The experimental system

All of our studies have been accomplished using a step-index MMF with radius 3=25.0±0.5μΐ7ΐ and numerical aperture (NA) 0.22 ± 0.02 (these values were provided by the fibre supplier Thorlabs). Our system is based on a standard digital holographic geometry in a Fourier regime where an SLM is imaged onto the back aperture of a microscope objective MO1 (20x, NA=0.25) by a 4-f system formed by lenses L3 and L4 (see Supplementary Figure S14 for a simplified scheme). The signal from the SLM (1064nm) is sent to three different zones within this telescope simultaneously. One part is separated by mirror M1 and used as a reference signal during the acquisition of TM. The remaining two are sent into a system of half-wave plates HWP1 -3 and polarising beam displacers PBD1 ,2 to create and overlap two fields with orthogonal linear polarisation states.

Even though this is sufficient for the generation of optical fields with an arbitrary spatial distribution of amplitude, phase and polarisation, a quarter-wave plate QWP1 is inserted into the geometry to gain higher efficiency and reduce noise, since most of the experimental work requires circular polarisation. Part of the signal sent to 01 is separated by the non-polarising beam splitter NPBS1 and imaged onto CCD1 to eliminate optical aberrations and measure SLM irradiance (detailed in the following section) and to monitor the field being sent into the MMF (the field recorded at CCD1 is a scaled copy of the field at the input facet of MMF). The optical signal leaving MMF is collected by microscope objective M02 (same parameters as M01) and focussed by tube lens L6. Circularly polarised modes are converted to linearly polarised ones by HWP2 and individual polarisation components are separated on polarising beamsplitter PBS and imaged on CCD2 and CCD3.

The reference signal is delivered to both CCD2 and CCD3 by single mode fibre SMF with collimator lenses L7 and L8 at each end. The reference optical pathway is merged with the imaging pathway before the polarisation components are separated by NPBS2. The geometry enables us to generate any optical field allowed to propagate in the optical fibre and to observe how it is transformed by its transport through the fibre. Crucially, the geometry allows us to synthesise and observe all aspects of the input and output optical fields respectively, i.e. their intensity, phase and polarisation. CCD4 was only employed in tests of proximal end imaging explained below in the section SR2.

Wavefront correction

All the experimental studies in this paper require beam control with extremely high fidelity, virtually free of optical aberrations. To eliminate aberrations in our very complex optical system we use a method developed in [14]. This algorithm efficiently manipulates phases of small segments of the SLM to gain the highest signal in the selected pixel of the CCD1 , signifying optimal focussing. The quality of the resulting wavefront correction can be influenced by CCD noise and SLM flicker noise, both of which can be statistically averaged out by longer data collection. However, with increased duration of the procedure we risk a strong influence of systematic drift. To overcome this trade off we measure optical aberrations with large number of different pixels (10 000) simultaneously and combine the results. As pixels are organised across an orthogonal grid, each of these measurements carries a specific linear phase-modulation (equivalent of a prism) which can be predicted and subtracted from the data before averaging. The enhancement is demonstrated in Supplementary Figure S15. Measurement of transformation matrix

In our previous studies [6, 7] we have defined input modes for the transformation matrix as beamlets originating from series of SLM segments. One considerable problem of this approach is a relatively weak signal that is often significantly affected by various sources of scattered light. For the purposes of this study we use a representation of optimally focussed spots at the input facets of the MMF (as explained above). This is beneficial not only for gaining much stronger signals (the whole power reflected from SLM is utilised for analysis of every mode) but also for simpler comparison with theoretical prediction. During the TM measurement in this regime the SLM splits the optical power between two pathways - one is sent to the SMF which acts as a reference signal, and the other results in an FP input mode at a selected location of the input facet of MMF. Due to the series of polarisation optical elements shown in Fig. S14, each FP can be generated by two different SLM modulations of the signal, leading to generation of a pair of mutually orthogonal polarisation states. Output FP modes are monitored on selected (grouped) pixels of CCD2 and CCD3, where the output fibre facet is imaged.

Their phases are established from interference with the reference signal delivered by a SMF. Employing the SLM to uniformly alter phase difference ΔΘ between the tested input mode and the reference field leads to a harmonic signal for each CCD pixel:

where Fm and Fref are the fields of the tested mode and the reference field, respectively. Typically we use more that one period of 2π in such measurement, e.g. Δθ =2π/4·[0,1 ,...,7].

Analysing the detected signal we can immediately isolate the distribution of phase of the tested input mode φ™ across the output modes. From the measured harmonic signal we can also establish the magnitude of AC and DC components that based on equation (10) allow for recovery of the output mode amplitudes. This way the measurement is not affected by the Gaussian envelope and small spatial intensity nonuniformity of the reference signal (resulting from interference effects at various optical components of the reference pathway), which can significantly affect measurement of the TM. Similarly to [6] we use a feedback loop in order to eliminate drift between MMF optical pathway and the reference optical pathway. The initial measurement of the TM contains 1089 input modes and 5625 output modes. While the input modes are measured sequentially the output modes are measured simultaneously, and therefore they are favourably over sampled in order to reduce noise. Processing raw data

The data from direct measurement of the TM ( ) are heavily affected by various sources of misalignment. Moreover, input and output modes are obtained with different sampling. Before such a TM can be used in comparison with theoretical prediction the data must be processed to address these issues. The amount of misalignment can roughly be estimated from total transmission profiles of the TM which in the ideal case (no misalignment) would be perfectly symmetric and focused on both sides of the optical system. These profiles are obtained by averaging absolute squares of either rows or columns of the TM and they signify position misalignment of the fibre end with respect to the optical pathway. Using a Fourier transform we can convert the measured TM into the representation of plane-waves which analogously provide transmission profiles, this time signifying tilts and defocus. The transmission profiles are shown in Supplementary Figure. S16 (a-d).

All misalignment operators can be represented by phase-only linear or quadratic modulations of the modes or their Fourier spectra and combined into one matrix for input Am and one matrix for output Aout of the fibre. The correction is then implemented by matrix multiplication of the TM from either side:$ = A _0Ut MA . The resulting pre-optimisation giving the best image symmetry and sharpness is presented in Supplementary Figure S16 (e-h).

Following this procedure the output modes are resampled using the scaling property of the Fourier transform in order to match the sampling of the input modes. Finally, redundant data (manifesting itself as dark boundaries in Supplementary Figure S16) are removed leading to a matrix with 729 input and 729 output modes presented in Figure S5. This operation is also achieved by matrix multiplication: Optimisation procedure

Our optimisation procedure gradually evolved over the course of this project to take into account a series of newly recognised factors. In the following we describe only the final one, optimised to study CP PIMs propagating through 100mm long fibre segments. With the well established theoretical description of light transport in MMFs [1 ] (reviewed in Supplementary Methods S1) and its modifications introduced above we can compute conversion matrix T that allows conversion of the measured TM into the representation of PIMs: M= TMT† . Due to remaining misalignment imperfections and insufficient precision in estimation of fibre parameters this does not lead to satisfactory results and further optimisation is necessary. While compensation for misalignment can be introduced by matrix multiplications of M, correction of fibre parameters requires calculation of new conversion matrix T at every iteration, which is time-consuming and M does not have fixed size which requires more elaborate optimisation criteria to converge. Therefore these two parts of the optimisation are decoupled. In each iteration (indexed by j) of the misalignment optimisation procedure we implement misalignment corrections to the experimentally measured TM before its conversion into PIMs: y = TC _0Ut MC T† . The quality metric that the algorithm optimises across the 12-dimensional space of free parameters is defined on the resulting matrix My as the fraction of the total optical power (absolute squares) carried by the TM's diagonal components, multiplied by PA (phase agreement introduced above).

During the optimisation of fibre parameters (that directly influence the conversion matrix T) the quality metric defined on the resulting My is a combination of four factors. This time we could not use the total power ratio carried by diagonal components as this strongly favoured small-sized matrices corresponding to low NA and smaller core diameters, which led to divergence of the optimisation algorithm. Neither the total power of diagonal components proved itself to be a suitable criterion as it favoured large-sized matrices corresponding to high NA and larger core diameters. The optimal power conversion was therefore achieved by reverse conversion of the My back into the space of diffraction limited points: = T^M _j T and comparison with the original experimentally measured one using a least squares . This value steeply falls while increasing both the NA and the diameter of fibre core but after passing the expected value it remains stationary. This was complemented by the introduction of a second criterion ensuring maximal uniformity of diagonal components. The standard deviation of diagonal components of My remains stationary until reaching the expected values and than rises steeply. Similarly, as in the previous case, we use PA as the third optimisation factor (PA ~ ¹ as we search for minimum in this case). Finally, we have found that much faster convergence is ensured by minimising the value of a fourth criterion which represents a total power of elements directly neighbouring to the main diagonal. The combination of these four criteria features a sharp minimum in the multidimensional space of all fibre parameters: NA, core diameter, diffusion lengths, and further refractive index corrections as described above. Both parts of the optimisation procedure are cycled several times until a stationary point is reached. The result of the algorithm is the set of optimal fibre parameters, optimised M and T and fine misalignment matrices Cm and Com . The commented Matlab code for this part together with detailed instructions and samples of pre-optimised experimental data samples (M) is available for download from [26].

The whole optimisation procedure is executed separately for polarisation components of M

We only seek optimal fibre parameters in the first step of analysing Mil together with misalignment factors. A record of this procedure is presented in Supplementary media SM1. This is followed by analysis of M+l where we only seek misalignment parameters (light follows different optical pathway as shown in Fig. S14) and for conversion into PIMs we use fibre parameters obtained in the previous step. For conversion of components Mil and M+l we know all the misalignment parameters of the systems from analysis of and M+l , therefore no further optimisation is needed.

Computing complex holographic modulations

Our LabView-based controlling system is designed to synthesise complex holographic modulations to be applied across the SLM. It uses the experimentally measured TM (M) to combine a series of input modes (in the representation of FPs) in order to achieve any output optical field leaving the optical fibre, with arbitrary distribution of amplitude, phase and polarisation. This is achieved using a modified version [31] of the Gerchberg-Saxton algorithm [32].

Synthesis of holograms for imaging

Since imaging is achieved by scanning a FP behind the output fibre facet, performing imaging with an experimentally measured TM is accomplished by sequential indexing of rows of M leading to synthesis of a single FP output mode at a time. For imaging with theoretically predicted transformation matrices M _th we can utilise pairs of matrices A, B and C obtained during optimisation of the data: M _th = B ^f A* .Here M _th represents any of the modifications introduced in the main text, including those for deformed fibres and axially shifted imaging planes. Generation of PIMs

Generation of individual PIMs uses the same procedure as described in the previous case. The only difference is in the transformation matrix that we load into the controlling interface. Instead of the fully reversed conversion we only convert the input modes as Af _h ^M = _th TC ^† 5 ^{† †} .The matrix f _h ^/ is a hybrid transformation between the input, represented using FPs and the output represented by PIMs. The complete set of PIMs recorded is presented in Supplementary Media SM2-SM5.

Supplementary methods

Theoretical description of light propagation in multimode fibres

The most common description of light propagation in MMF is based on the assumption that the polarisation state of the optical fields propagating through them remains unchanged. Such scalar and paraxial modelling based on the theory of weakly guiding modes leads to the identification of a set of orthogonal linearly polarised (LP) propagation-invariant modes (PIMs) with radial dependence described simply by splicing a Bessel function of the first kind and a modified Bessel function of the second kind defined across the fibre core and the cladding, respectively. The requirement of continuity of the function together with its first derivative leads to a set of discrete solutions (being indexed by letter m in this paper) for every order of the Bessel functions /. The intensities of the fields are azimuthal!y independent and their phase varies with the azimuthal coordinate φ as e' ¹ * . This simplified approach serves well to explain several important light transport processes, however linearly polarised light in an MMF with parameters such as those used in our experiments will only retain its polarisation state over distances of a few tens of millimetres. A rigorous vectorial modelling also predicts the existence of a set of PIMs; their polarisation is however no longer uniform, with orientation periodically changing across the mode cross-section. Very fortunately, for the majority of the modes the theory predicts degeneracy in the values of propagation constants (axial components of the k vector). A linear combination of such degenerate modes therefore also leads to propagation-invariant optical fields that have an almost identical distribution of the field to the simplified scalar LP PIMs and uniform polarisation across the mode cross-section, but the polarisation state is circular. Due to these qualities, we chose this representation for our experimental study since it is much easier to generate circularly polarised (CP) modes than modes with spatially non-uniform distribution of polarisation. The only exception to this are azimuthally and radially polarised modes sometimes termed as hedgehogs and bagels. While their superposition leads to modes with / = ±1 , here such degeneracy does not occur. In our model they are approximated by circularly polarised modes with averaged propagation constants. Intensity of such modes is still azimuthally independent but their energy periodically oscillates between opposite polarisation states, therefore they are no longer propagation invariant. As shown in Fig. 6, this insufficiency has a negligible influence on imaging performance.

The representation of discrete circularly polarised modes can be conveniently expressed using a quantum formalism describing each mode with the trio of integer quantum numbers /, m and a. In this formalism, / and a denote the amount of orbital and spin angular momentum, respectively, in the units of fi per photon [25]. This representation is also very convenient to observe the effect of spin-orbit interaction. This effect, known from quantum physics, causes two modes with the same but nonzero / and opposite spin states σ = -1 and σ = 1 to propagate with slightly different phase velocities.

The similarity of the vectorial CP modes with the scalar LP modes makes this representation very suitable for implementing corrections for deviations from the ideal step-index profile that can strongly influence the experimental observations. Again, the influence of these aberrations does not considerably affect the distribution of optical fields of the modes. The most pronounced effect is a small variation of propagation constants (relative deviations are in the order of 10 ^"6 ) which can however significantly affect the output phases of the modes over distances of only a few tens of millimetres. As we will see, these small deviations of propagation constants can be modelled sufficiently precisely and efficiently using scalar LP modes by employing perturbation calculus. Our modelling technique is therefore a hybrid solution employing exact vectorial CP modes with propagation constants being corrected for imperfections in the step-index profile.

Vector modes

CP PIMs constitute the simplest set of modal fields that take into account the vector character of light propagation. In a step-index profile circular fibre the transverse electric field components of CP PIMs can be expressed as:

ι,τη,σ = (* + ia )F _lim { exp{ilo , ? = _c + δβ _ι± (S1)

where R = r/a is a dimensionless radial parameter with r being the radial coordinate and a the radius of the fibre core, β are the propagation constants of the corresponding modes and we have denoted the polarisation corrections to scalar propagation constants sc by δβι,+ for Ι,σ with equal signs and by δβι,- for Ι,σ with opposite signs. £ and† are unit Cartesian coordinate vectors and Fim(R) is the radial profile of the mode defined by Bessel functions of the first (Ji ) and second kind (Ki ) :

Parameter u is proportional to the transverse wavevector of a considered mode, u = a. ]kln _ore - β ² , where ko = ω/c is the vacuum wavenumber and n core is the refractive index of the fibre core.

Equation (S1 ) describes the CP modes. As we have mentioned, they are the eigenmodes of the fibre with the exception of the modes for which Ι+σ = 0, i.e., the ones whose total angular momentum is zero. These exceptional modes are not even approximate solutions of the vector wave equation and as such they are not the eigenmodes of the fibre; however, we can construct the eigenmodes by their equal superpositions:

ψ hedgehog, m = (ψ l ,m,-1 +ψ -l ,m, i )/2

ψ bagel.m≡ (ψ 1 ,m,-1 ~ ψ -1 ,m, 1 )/2. (S3)

Therefore, strictly speaking, in the case of / +σ = 0, the propagation constant β in Eq. (S1 ) cannot be defined. For the same reason we cannot define the correction δβι ,- as a correction to the propagation constant of the mode ψ i , _m ,-/ . We will therefore reserve this symbol for a correction defined with respect to the hedgehog mode, and the correction for the bage! mode will be denoted by δβ' i - .

The fact that the modes ψ and ψ themselves are not the eigenmodes of the fibre shows itself in a very specific way. if the mode ψ i , _m , -? is injected into the fibre, then after a certain distance (say Lm ) it is completely transformed into the mode ψ -i , _m , i . Then after propagating through an additional distance of Lm , the original mode ψ i . _m ,-? is recovered, and the state continues to oscillate between these two modes along the fibre. This way, the case of Ι+σ = 0 is the only one where orbital and spin angular momenta are not conserved separately along the fibre, however their sum is still conserved due to the rotational symmetry of the fibre. The bagel and hedgehog modes are in fact members of a complete, although slightly more complicated, representation of the modes of the fibre which we can call "generalised bagel and hedgehog" (GBH PIMs). This representation of modes in a step index profile fibre is constructed as a linear superposition of CP PIMs: +HEf™ il≥ 0)≡ + l/>-i,m,-l), /? = Psc + δβ

I -

EHf∑l ⁿ _m {l≥ 1) = - O m,-! + Ψ-Ι,ηΛβ = fi _sc + δβ I +HE? _m (l≥0)≡ 0Km,l + Ψ-1,τη,-ΐ),β = Psc + δβ ι,-

EHffinQ≥ 1)≡ (Vi,m,-i + Ψ-Ι,τη,ι),β = sc + δβ'

The first index (/+1 or M) of these modes denotes the magnitude of the total angular momentum, and δβ i ,± are again corrections of the propagation constants with respect to the scalar theory. However, in the case of / +σ = 0 the corrections split into two values, different for the even and odd modes (hedgehogs and bagels, respectively), which is emphasised by

I - ' I - ' the prime in ' ^ in the last equation. For /+σ≠ 0, the corrections βι,- and ' ^ are equal.

Specifically, the corrections can be expressed in a simple form: for all / > 0, δβι,+ = h - ; for

/≠ 1 , δβι,- = ^l ' ^~ ' = Ii + h , and for / = 1 , δβι - = 2{h + ) and ^l ' ^~ ' = 0. Here h and l ₂ are simple integrals [1]:

l(2A) ^{3 2} f ^∞ I r ^∞

^{/2 = "} ^^ J ₀ ^RF Lm R)dR/l RFf _m {R)dR, where V = ako NA is a waveguide parameter with NA being the numerical aperture of the fibre, Δ= (n _ore - n _lad ) / (2n _ore ) is a profile height parameter and f(R) is a dimensionless profile shape function. For the case of step-index fibre f (R) is defined as:

f _(R) _ [°f° ^rR < ¹

K ⁾ ~ \lforR≥ 1

This makes the calculation of the above integrals really simple as (R) = 5(R-1). As a consequence, only the normalisation integrals have to be calculated. Table 1 summarises the corrections for the CP modes as well as their corresponding GBH modes in terms of // and h integrals. In case of / = 0 that is shown separately in the table, = 0 and therefore the correction is just . The correction h is opposite for the modes ψ i, _m , i and ψ i, _m ,-i ,

Table 1. Corrections δβί to the scalar propagation constants εο .

CP modes Corresponding δβι / values

GBH

ψ 0,m,±1 jt peven rj codd

n ^c l,m > ^{n D} l,m h /=0

ψ±Ι,πι,±1 u neven jj cOdd I1 - I2 l >1

(hedgehog) E ¾ ^e ", ΤΜ _0ιΤη 2(h - l ₂ ) I = 1

(bagel) rp iiodd T*p

c ⁿ 0,m > ^{1 D} 0,m 0 I = 1

xP which means that the propagation constants of these modes are slightly different. This is a demonstration of the spin-orbit interaction, which causes rotation of the polarisation vector direction if linearly polarised light with a given I propagates through the fibre, or rotation of the intensity patterns when the superposition ψ ι,η,,σ +ψ -im.a propagates.

For completeness, we have also compared the weak guidance corrections δβί (Figs.

S1a,b) calculated using the integrals li and I2 and used them to evaluate higher order vectorial corrections (Figs. S1 e,f) from:

βνβο -βεο = δβί + (higher order corrections), (S5)

shown in(Figs.S1c,d), where βνβο is a fully vectorial propagation constant. The higher order corrections are on the order of =10 ^~6 .

The modal fields above constitute a complete approximate set of solutions of the vector- wave equation and have similar ΕΗ,ΗΕ-nomenclature to that used with the full vectorial approach. Even though the GBH PIMs form a complete basis (including the problematic I +σ = 0 case), their experimental drawback is that their non-uniform polarisation changes orientation periodically across the mode cross-section. As mentioned before, this led us to choose the CP PIMs representation for our experimental study since it is much easier to generate circularly polarised modes than modes with spatially non-uniform distribution of polarisation. This choice of basis has a negligible influence on imaging performance.

Due to the very high accuracy of weak guidance approximation in our case, the propagation constants β for the above GBH modes, and also for CP modes, obtained directly from a full vectorial approach, yield virtually identical results to the scalar wave equation with weak guidance polarisation corrections applied. The latter approach, however, has significant advantages over the full vectorial approach. For example, the fibre exhibits small variations of refractive index. These small variations can be readily accounted for by utilising the perturbation theory applied to scalar modes, as we will do in the following. The same cannot be done easily in the full vectorial description. Due to this clear advantage, we choose the scalar theory with weak guidance and perturbation corrections as a main tool in our calculations.

Deformation operator in 3-D

The theoretical model predicts that DO does not change if the plane in which the deformed fibre lies is oriented differently, provided that the output imaging system remains oriented accordingly to the deformation imposed. If this condition is not fulfilled, the DO is accompanied by an orientation operator that rotates all output modes by the differing angle y. To confirm this behaviour we have revolved the deformed fibre (V) together with the imaging apparatus around the axis of the initially straight fibre as shown in Fig. 4a by various angles β. For this case the model predicts that y =β. Diagonal components of DO for these cases are presented in Fig. S2b. To isolate the orientation operator from these data we can efficiently utilise the fact that for γ =0, DO features symmetry of ύ(+1, m) = ΰ(-1, m).Therefore we have searched for the common angle with which output modes have to be rotated in order to maximise this symmetry. The isolated values of γ against the actual orientation β are shown in Fig. S2c. The resulting DOs with the orientation operator eliminated are presented in Fig. S2d, clearly confirming the orientation invariance of deformation operators.

Proximal end imaging

During image acquisition of images the total transmitted signal was detected using a bucket detector (integrating signal across a CCD 3) for each of the output FP generated across the grid of 75x75 positions (Fig. S3a). As some may argue that this is not a real endoscopic method (the detector is not at the proximal side of the instrument), we have also integrated the light returning from the fibre at the proximal end using CCD4.

The corresponding image (Fig. S3b) clearly exhibits unwanted interference effects present in our imaging pathway. These originate from the light being retro-reflected on many surfaces of our geometry (including fibre facets). This does not pose any obstacles for methods relying on fluorescence emission where the retroreflected excitation wavelength can be filtered out and only the emission signal is detected. To minimise this effect we have taken a reference image of a uniform reflective surface (Fig. S3c). The difference of the two, shown in Fig. S3d, gives an image analogous to distal end imaging (Fig. S3a) but due to the more elaborate procedure it is more affected by noise. For this reason we have chosen to use the distal end imaging in our demonstrations as it makes observation of subtle effects on imaging quality more evident.

Imaging in arbitrary distance from the fibre end

Here we show that with highly precise knowledge of the fibre TM (particularly phases of output modes) we can complement the TM with a free-space propagating operator and thus modify the location of imaging plane numerically. As shown in the sequence Fig. S4, the further the image plane is set from the fibre, the larger is the accessible field of view but accordingly the resolution obtained is lower. The amount of information imaged remains the same in all the cases shown.

Figure S1. a-b, Vector corrections to β calculated from li and integrals, c-d, Difference between full vector ee and scalar p _sc . e-f, Higher order vectorial corrections deduced by subtracting (a-b) from (c-d).

Figure S2. Influence of deformation orientation, a, Arrangements of deformed fibre used in the experiment with their roman ID number, b, Influence of curvature orientation on PIMs. c, Argument of isolated rotation operator, d, Influence of DO on PIMs after removing rotation operator.

Figure S3. Demonstration of proximal end imaging, a, Direct distal end imaging, b, Proximal end imaging with USAF target, c, reference image of reflective surface, d, Difference between (c) and (b).

Figure S4. Imaging with numerically updated TM to image at various distances behind the fibre end. The white scale bar corresponds to the size of fibre core, 50 pm.

Figure S5. Experimental transformation matrix data after processing introduced in Online Methods This example represents TM for 10mm long segment of fibre in the representation of linearly polarised FPs. a, Polarization component S→ S. b, Polarization component S→ P. c, Polarization component P→ S. d, Polarization component P→ P. Figure S6. Conversion matrix, λ = 1064nm, NA=0.22 and core diameter d = 50pm.

Figure S7. Converted transformation matrix for 10mm long segment of fibre. TM is expressed in the representation of LP PIMs

Figure S8. Converted transformation matrix for 100mm long segment of fibre. TM is expressed in the representation of LP PIMs

Figure S9. Converted transformation matrix for 100mm long segment of fibre. TM is expressed in the representation of CP PIMs

Figure S10. Converted transformation matrix for 300mm long segment of fibre. TM is expressed in the representation of CP PIMs

Figure S11. Optical phases of PIMs, data for the 300mm long fibre, a, b and c, Assumed profile of refractive index, phase agreement and difference between experimentally obtained and theoretically predicted phases of PIMs respectively for the case of ideal step-index fibre. d, e and f, corresponds to a model with included dopant diffusion, g, h and i, represents correction for diffusion as well as fine index modulation (three parameters) across the fibre core. Figure S12. Experimentally measured deformation operator corresponding to fibre deformation (V) from Figure 4.

Figure S13. Theoretically predicted deformation operator corresponding to fibre deformation (V) from Figure 4.

Figure S1 . The experimental geometry.

Figure S15. Enhancement of wavefront correction, a, Single pixel correction, b, Enhanced wavefront correction obtained averaging 10 000 single pixel corrections.

Figure S16. Coarse misalignment compensation, a, Transmission profile of the input facet area, scanned during calibration, b, Transmission profile of input angular spectrum, c, Transmission profile of output facet, d, Transmission profile of output angular spectrum. Both ( a ) and ( b ) were up-sampled to match the resolution of output transmission profiles. The large shift of the transmission spectrum towards the top left corner in ( d ) was introduced into the system deliberately by adjusting the angle of incidence of the reference signal onto the CCDs. This was done in order to eliminate strong interference effects caused by multiple reflections of the reference signal within our system, e-h, The equivalent of ( a- d ) after applying coarse correction for misalignment.

Media SMI . Record of progression of optimisation procedure introduced in Online Methods

Media SM2. Experimentally generated full set of PIMs as introduced in Online Methods. PIMs are generated in circular polarisation and propagate through 10mm long segment of fibre. Colour corresponds to individual polarisation components. The left column shows images of the modes prior they enter the optical fibre (recorded by CCD1). The right column shows PIMs at the output, recorded by CCD1 and CCD2.

Media SM3. Experimentally generated interfering pairs of PIMs with opposite I indices.

PIMs are generated in circular polarisation and propagate through 10mm long segment of fibre.

Media SM4. Experimentally generated full set of PIMs as introduced in Online Methods. PIMs are generated in circular polarisation and propagate through 100mm long segment of fibre.

Media SM5. Experimentally generated interfering pairs of PIMs with opposite I indices.

PIMs are generated in circular polarisation and propagate through 100mm long segment of fibre.

Media SM6. Invariance of deformation orientation.

Media SM7. Imaging on the distal and the proximal side of the multimode fibre, a,

Record of CCD3 (at the distal end of the fibre) during imaging. We deliberately image the FPs onto the CCD chip out of focus. This is to spread the total power over a larger area thus allowing larger power to be detected and reducing overall noise, b, Progress of the image acquisition from data obtained by CCD3. c, Record of CCD4 (at the proximal side of the fibre) during imaging. Here we collect the light returning from the fibre after reflecting off the object. Due to its propagation through the fibre the signal appears scrambled, d, Progress of the image acquisition from data obtained by CCD4. Due to space constraints both (a) and (c) represent only a fraction of the whole CCD chips' areas.

References

[I] Snyder, A. W. & Love, J. Optical Waveguide Theory (Springer Science & Business Media, 1983).

[2]Gloge, D. Weakly guiding fibres. Appl. Optics 10, 2252-2258 (1971 ).

[3]Snitzer, E. Cylindrical dielectric waveguide modes. J. Opt.Soc. Am. 51 , 491-498 (1961 ).

[4] Liberman, V. S. & Zel'dovich, B. Y. Spin-orbit polarization effects in isotropic multimode fibres. Pure Appl. Opt. 2, 367-382 (1993).

[5] Di Leonardo, R. & Bianchi, S. Hologram transmission through multi-mode optical fibers. Opt. Express 19, 247-254 (201 1 ).

[6] Cizmar, T. & Dholakia, K. Shaping the light transmission through a multimode optical fibre: complex transformation analysis and applications in biophotonics. Opt. Express 19,18871 -18884 (201 1 ).

[7] Cizmar, T. & Dholakia, K. Exploiting multimode waveguides for pure fibre-based imaging. Nat. Commun. 3, 1027 (2012).

[8] Bianchi, S. & Di Leonardo, R. A multi-mode fiber probe for holographic micromanipulation and microscopy. Lab Chip 12, 635-639 (2012).

[9] Choi, Y. et al. Scanner-free and wide-field endoscopic imaging by using a single multimode optical fiber. Phys. Rev. Lett. 109, 203901 (2012).

[10] Papadopoulos, I. N., Farahi, S., Moser, C. & Psaltis, D. Focusing and scanning light through a multimode optical fiber using digital phase conjugation. Opt. Express 20, 10583- 10590 (2012).

[I I] Nasiri, R., Mahalati, Gu, R. Y. & Kahn, J. M. Resolution limits for imaging through multi- mode fiber. Opt. Express 21 , 1656-1668 (2013).

[12] Vellekoop, I. M. & Mosk, A. P. Focusing coherent light through opaque strongly scattering media. Opt. Lett. 32, 2309-231 1 (2007).

[13]Popoff, S. M. et al. Measuring the transmission matrix in optics: An approach to the study and control of light propagation in disordered media. Phys. Rev. Lett. 104, 100601 (2010).

[14] Cizmar, T., Mazilu, M. & Dholakia, K. In situ wavefront correction and its application to micromanipulation. Nature Photon. 4, 388-394 (2010).

[15] Popoff, S., Lerosey, G., Fink, M., Boccara, A. C. & Gigan, S. Image transmission through an opaque material. Nat. Commun. 1 , 81 (2010).

[16] Vellekoop, I. M., Lagendijk, A. & Mosk, A. P. Exploiting disorder for perfect focusing. Nature Photon. 4, 320-322 (2010).

[17] Ji, N., Milkie, D. E. & Betzig, E. Adaptive optics via pupil segmentation for high-resolution imaging in biological tissues. Nat. Methods 7, 141-147 (2009).

[18] Papadopoulos, I. N., Farahi, S., Moser, C. & Psaltis, D. High-resolution, lensless endoscope based on digital scanning through a multimode optical fiber. Biomed. Opt. Express 4, 260 (2013).

[19] Farahi, S., Ziegler, D., Papadopoulos, I. N., Psaltis, D. & Moser, C. Dynamic bending compensation while focusing through a multimode fiber. Opt. Express 21 , 22504-22514 (2013).

[20] Gambling, W. A., Payne, D. N. & Matsumura, H. Mode conversion coefficients in optical fibers. Appl. Optics 14, 1538-1542 (1975).

[21] Friesem, A. A., Levy, U. & Silberberg, Y. Parallel transmission of images through single optical fibers, in Proc. IEEE, 208-221 (1983).

[22] Kreysing, M. et al. Dynamic operation of optical fibres beyond the single-mode regime facilitates the orientation of biological cells. Nat. Commun. 5, 5481 (2014).

[23] von Hoyningen-Huene, J., Ryf, R. & Winzer, P. LCoS-based mode shaper for few-mode fiber. Opt. Express 21 , 18097-18110 (2013).

[24] Carpenter, J., Eggleton, B. J. & Schroder, J. 110x110 optical mode transfer matrix inversion. Opt. Express 22, 96-101 (2014).

[25] Leach, J., Padgett, M., Barnett, S., Franke-Arnold, S. & Courtial, J. Measuring the orbital angular momentum of a single photon. Phys. Rev. Lett. 88, 257901 (2002).

[26] http://complexphotonics.dundee.ac.uk/ .

[27] Bliokh, K. Y., Niv, A., Kleiner, V. & Hasman, E. Geometrodynamics of spinning light. Nature Photon. 2, 748-753 (2008).

[28] Lyytikainen, K. et al. Dopant diffusion during optical fibre drawing. Opt. Express 12, 972-977 (2004).

[29] Gibson, B. C. et al. Controlled modification and direct characterization of multimode- fiber refractive-index profiles. Appl. Optics 42, 627-633 (2003).

[30] Skinner, B. J. & Appleman, D. E. Melanophlogite, a Cubic Polymorph of Silica. Am. Mineral. 48, 854-867 (1963).

[31] Cizmar, T. & Dholakia, K. Tunable bessel light modes: engineering the axial propagation. Opt. Express 17, 15558-15570 (2009).

[32] Gerchberg, R. W. & Saxton, W. O. A practical algorithm for the determination of the phase from image and diffraction plane pictures. Optik 35, 237-246 (1972).