Wavefront shaping in multimode fibers by transmission matrix engineering

One of the greatest challenges in utilizing multimode optical fibers is mode-mixing and inter-modal interference, which scramble the information delivered by the fiber. A common approach for canceling these effects is to tailor the optical field at the input of the fiber to obtain a desired field at its output. In this work, we present a new approach which relies on modulating the transmission matrix of the fiber rather than the incident light. We apply computer-controlled mechanical perturbations to the fiber to obtain a desired intensity pattern at its output. Using an all-fiber apparatus, we demonstrate focusing light at the distal end of the fiber and conversion between fiber modes. Since in this approach the number of degrees of control can be larger than the number of fiber modes, it allows simultaneous control over multiple inputs and multiple wavelengths.


Introduction
In recent years, multimode optical fibers (MMFs) are at the focus of numerous studies aiming at enhancing the capacity of optical communications and endoscopic imaging systems [1,2]. Ideally, one would like to utilize the transverse modes of the fiber to deliver information via multiple channels, simultaneously. However, inter-modal interference and coupling between the guided modes of the fiber result in scrambling between channels. One of the most promising approaches for unscrambling the transmitted information is by shaping the optical wavefront at the proximal end of the fiber in order to get a desired output at the distal end. Demonstrations include compensation of modal dispersion [3][4][5], focusing at the distal end [6][7][8][9][10], and delivering images [11][12][13] or an orthogonal set of modes [14,15] through the fiber.
Typically in wavefront shaping, the incident wavefront is controlled using spatial light modulators (SLMs), digital micromirror devices (DMDs) or nonlinear crystals. In all cases, the shaped wavefront sets the superposition of guided modes that is coupled into the fiber. For a fixed transmission matrix (TM) of the fiber, this superposition determines the field at the output of the fiber, as depicted in Fig. 1(a)). Hence, in a fiber that supports N guided modes, wavefront shaping provides at most N complex degrees of control. However, many applications require the number of degrees of control to be larger than the number of modes. For example, one of the key ingredients for spatial division multiplexing is mode converters, which require simultaneous control over the output field of multiple incident wavefronts. To this end, complex multimode transformations were previously demonstrated by applying phase modulations at multiple planes [16][17][18][19]. However, this requires free-space propagation between the modulators, thus limiting the stability of the system and increasing its footprint.
In this work we propose and demonstrate a new method for controlling light at the output of MMF, which does not rely on shaping the incident light and that can be implemented in an all-fiber configuration. Inspired by the ongoing efforts to generate on-chip mode converters by manipulating modal interference in multimode interferometers [20][21][22], we directly control the light propagation inside the fiber to manipulate its TM, allowing us to generate a desired field at its output ( Fig. 1(b)). Since the TM is determined by O N 2 complex parameters, TM-shaping provides access to much more degrees of control than shaping the incident wavefront.
To control the fiber's TM, we apply computer controlled bends at multiple positions along the fiber. Since the stress induced by the bends changes the boundary conditions of the system, it modifies the TM such that different bends yield different speckle patterns at the distal end (Fig. 1(c)). We can therefore obtain a desired field at the output of the fiber by imposing a set of controlled bends, without modifying the incident wavefront. Since in this approach the input field is fixed, it does not require an SLM or any other free-space component. Such an all-fiber configuration is especially attractive for MMF-based applications that require high throughput and an efficient control over the field at the output of the fiber. As a proof-of-concept demonstration of TMshaping, we demonstrate focusing at the distal end of the fiber, and conversion between the fiber modes. Figure 1: Shaping the transmission matrix of multimode optical fibers. (a) The conventional method for wavefront shaping in complex media, performed e.g. by using an SLM and free space optics to tailor the incoming wavefront at the proximal end of the multimode fiber. (b) Proposed method for light modulation, in which the transmission matrix of the medium is altered, e.g. by performing perturbations on the fiber itself. (c) Illustration of the sensitivity of the output pattern on the fiber geometry. Three different configurations of the fiber (depicted by red, green and blue curves), correspond to three different speckle patterns at the output of the fiber. Since the input field coupled into the fiber is fixed, the different output patterns correspond to different transmission matrices of the fiber.

Experimental Techniques Principle
Our method relies on applying controlled weak local bends along the fiber. To this end, we use an array of computer-controlled piezoelectric actuators to locally apply pressure on the fiber at multiple positions [23,24].
The TM of the fiber depends of the curvatures of the bends, which are determined by the travel of each actuator. To obtain a target pattern at the distal end, we compare the intensity pattern recorded at the output of the fiber with a desired target pattern. Using an iterative algorithm, we search for the optimal configuration of the actuators, i.e. the optimal travel of each actuator, that maximizes the overlap of the output and target patterns.

Experimental Setup
The experimental setup is depicted in Fig. 2. A HeNe laser (wavelength of λ = 632.8 nm) is coupled to an optical fiber, overfilling its core. We placed 37 piezoelectric actuators along the fiber. By applying a set of computer-controlled voltages to each actuator, we controlled the vertical displacement of the actuators. Each actuator bends the fiber by a three-point contact, creating a bell-shaped local deformation of the fiber, with a curvature that depends on the vertical travel of the actuator (see Figs. 2(b,c)). For the maximal curvature we applied (R ≈ 10 mm), we measured an attenuation of few percent per actuator due to bending loss. The spacing between nearby actuators was set to be at least 3 cm, which is larger than d 2 λ for d the core's diameter, such that the interference pattern inside the fiber between two adjacent actuators is uncorrelated. At the distal end, a CMOS camera records the intensity distribution of both the horizontally and vertically polarized light.
We used two types of multimode fibers: a fiber supporting few modes for demonstrating mode conversion, and a fiber supporting numerous modes for demonstrating focusing. For the focusing experiment, we used a 2 meter-long graded-index (GRIN) multimode optical fiber with numerical aperture (NA) of 0.275 and core diameter of d M M F = 62.5 µm (InfiCor OM1, Corning). The fiber supports approximately 900 transverse modes per polarization at λ = 632.8 nm (V ≈ 85), yet we used weak focusing at the fiber's input facet to excite only ≈ 280 modes. For the experiments with the few mode fiber (FMF), we used a 5 meter-long step-index (SI) fiber, with an NA of 0.1 and core diameter of d F M F = 10 µm (FG010LDA, Thorlabs). In principle, at our wavelength the fiber supports 6 modes per polarization, (V ≈ 5).

Optimization Process
The curvature of the bends, set by the travel of each actuator, modifies how light propagates through the fiber and thus determines the speckle pattern that is received at the distal end. We can therefore define an optimization problem of finding the voltages that should be applied to the actuators, to receive a given target pattern at the output of the fiber. The distance between the target and each measured pattern is quantified by a cost function, which the algorithm iteratively attempts to minimize. The fiber is pressed by two pins that are attached to each actuator, and one pin which is placed below it, creating a three-point contact. A computer-controlled voltage that is applied on each actuator sets its travel and defines the curvature of the local deformation it poses on the fiber. L, lens; M, mirror; PBS, polarizing beamsplitter; CMOS, camera.
For M actuators, the solution space is an Mdimensional sub-space, defined by the voltages range and the algorithm's step intervals, and can be searched using an optimization algorithm. While the optical system is linear in the optical field, the response of the actuators, i.e. the modulation they pose on the complex light field, is not linear in the voltages.
Moreover, since a change in the curvature of an actuator at one point along the fiber affects the interference pattern at all of the following actuators positions, the actuators cannot be regarded as independent degrees of control. Similar nonlinear dependence between degrees of control is obtained, for example, for wave control in chaotic microwave cavities [25]. Out of the wide range of iterative optimization algorithms that can efficiently find a solution to such nonlinear optimization problems, we chose to use Particle Swarm Optimization (PSO) [26], as on average it achieved the best results out of the algorithms we tested (See the Supplementary Material for more details regarding the use of PSO).

Focusing at the Distal End of the Fiber
To illustrate the concept of shaping the intensity patterns at the output of the fiber by controlling its TM, we first demonstrate focusing the light to a sharp spot at the distal end of the fiber. We excite a subset of the fiber modes by weakly focusing the input light on the proximal end of the fiber. Due to inter-modal interference and mode mixing, at the output of the fiber the modes interfere in a random manner, exhibiting a fully developed speckle pattern ( Fig. 3(a)). Based on the number of speckle grains in the output pattern, we estimate that we excite the first 280 guided fiber modes.
To focus the light to some region of interest (ROI) in the recorded image, we run the optimization algorithm to enhance the total intensity at the target area. We define the enhancement factor η by the total intensity in the ROI after the optimization, divided by the ensemble average of the total intensity in the ROI before the optimization. The ensemble average is computed by averaging the output intensity over random configurations of the actuators, and applying an additional azimuthal integration to improve the averaging.
We start by choosing an arbitrary spot in the output speckle pattern of one of the polarizations. We define a small ROI surrounding the chosen position, in an area that is roughly the area of a single speckle grain, and run the optimization scheme to maximize the total intensity of that area. Fig. 3 depicts the output speckle pattern of the horizontal polarization before ( Fig. 3(a)) and after ( Fig. 3(b) the optimization, using all 37 actuators. The enhanced speckle grain is clearly visible and has a much higher intensity than its surroundings, corresponding to an enhancement factor of η = 25.
We repeat the focusing experiment described above with a varying number of actuators M . When a subset of actuators is used, the remaining are left idle throughout the optimization. Fig. 3(d) summarizes the results of this set of experiments, showing the obtained enhancement factor η grows linearly with the number of active actuators M . It is instructive to compare this linear scaling with the well-known results for focusing light through random media using SLMs or DMDs. Vellekoop and Mosk have shown that when the number of degrees of control (i.e. independent SLM or DMD pixels) is small compared to the effective number of transverse modes of the sample, the enhancement scales linearly with the number of degrees of control. The slope of the linear scaling α depends on the speckle statistics and on the modulation mode [27][28][29]. For Rayleigh speckle statis-tics, as in our system (see Supplementary Material), the slopes predicted by theory are α = 1 for perfect amplitude and phase modulation, α = π 4 ≈ 0.78 for phaseonly modulation [29]. Experimentally measured slopes, however, are typically smaller, mainly due to technical limitations such as finite persistence time of the system, unequal contribution of the degrees of control, and statistical dependence between them. Interestingly, we measure a slope of α ≈ 0.71, which is close to the theoretical value for phase-only modulation for Rayleigh speckles, and higher than typical experimentally measured slopes (e.g. α ≈ 0.57 in [30]). Naively, one could expect a lower slope in our system, since as discussed above, in our configuration the degrees of control are not independent. The large slope values we obtain may indicate that the bends change not only the relative phases between the guided modes (corresponding to phase modulation), but also their relative amplitudes (corresponding to amplitude modulation), via mode-mixing and polarization rotation.
To further study the linear scaling, we performed a set of numerical simulations. We used a simplified scalar model for the light propagating in a GRIN fiber, in which the fiber is composed of multiple sections, where each section is made of a curved and a straight segment. The curved segments simulate the bend induced by an actuator, and the straight segments simulate the propagation between actuators (see Supplementary Material for more details). As in the experiment, we use the PSO algorithm to focus the light at the distal end of the fiber. The numerical results exhibit a clear linear scaling, with slopes in the range of 0.57-0.64 (see Fig. S3 in Supplementary Material). Simulations for fibers supporting N = 280 modes, roughly the number of modes we excite in our experiment, exhibit a slope of α ≈ 0.64, slightly lower than the the experimentally measured slope.
As in experiments with SLMs, focusing is not limited to a single spot. To illustrate this, we used the optimization algorithm to simultaneously maximize the intensity at two target areas. Fig. 3(c) shows a typical result, exhibiting an enhancement which is half of the enhancement obtained when focusing to a single spot, as expected by theory [28]. In principle, it is possible to focus the light to an arbitrary number of spots, yet in practice we are limited by the number of available actuators.

Mode Conversion in a Few Mode Fiber
In the previous section, we demonstrate the possibility to use our system as an all-fiber SLM, i.e. to shape an output complex field by modifying the relative complex weight of the propagating modes. In the following, we show that we can go further by studying the feasibility of TM-shaping to tailor the output patterns in the fewmode regime, where the number of fiber modes is comparable with the number of actuators. Specifically, we are interested in converting an arbitrary superposition of guided modes to one of the linearly-polarized (LP) modes supported by the fiber. To this end, we utilize the PSO optimization algorithm to find the configuration of actuators that maximizes the overlap between the output intensity pattern and the desired LP mode. The target LP modes of the step-index fiber were computed numerically for the parameters of our fiber, and scaled to match the transverse dimensions of the fiber image. Fig. 4 presents a few examples of conversions between LP modes using 33 and 12 actuators. A mixture of LP 01 and LP 11 at two different polarizations can be converted to LP 11 in one polarization ( Fig. 4(a)). Alternatively, a horizontally polarized LP 11 mode can be converted to a superposition of horizontally LP 01 and a vertically polarized LP 11 (Fig. 4(b)). The Pearson correlation between the target and final patterns in these examples is 0.93. Similar results are obtained when we run the optimization with fewer active actuators, with a negligible reduction in the correlation between the target and final pattern. For example, with 12 actuators we observe correlations of 0.90 for the conversion presented in Fig. 4(c). Optimization with less than 12 actuator shows poorer performance, as the number of actuators becomes comparable with the number of guided modes.

Discussion
Controlling the transmission matrix of a multimode fiber, rather than the wavefront that is coupled to it, opens the door for an unprecedented control over the light at the output of the fiber. Since the number of degrees of control, the number of actuators in our implementation, is not limited by the number of fiber modes N , it can allow simultaneous control for orthogonal inputs and/or spectral components. In fact, if O(N 2 ) degrees of control are available, one can expect generating arbitrary N × N transformations between the input and output modes. Over the past two decades there is an ever-growing interest in realizing reconfigurable multimode transformations, for a wide range of applications, such as quantum photonic circuits [22,[31][32][33][34] optical communications [18,35], and nanophotonic processors [20,36]. These realizations require strong mixing of the input modes, as the output modes are arbitrary superpositions of the input modes. The mixing can be achieved, for example, by diffraction in free-space propagation between carefully designed phase plates [16][17][18][19], a mesh of Mach-Zehnder interferometers with integrated modulators [22], engineered scattering elements in multimode interferometers [20,21], or scattering by complex media [25,37]. In our implementation, we rely on the natural mode mixing and inter-modal interference in multimode fibers, allowing implementation using standard commercially available fibers.
The main limitation our current proof-of-concept suffers from is the achievable modulation rates. The piezobased implementation limits the achievable modulation rates. The response time of the system to abrupt changes of the piezos is approximately 30 ms (see Supplementary Material), allowing in principle for modulation rates as high as 30 Hz. In practice, our system works at slower rates (≈5 Hz), mainly due the latency of the piezoelectric actuators and the camera. The total optimization time for the focusing experiments is 50 minutes, and 12 − 15 minutes for the mode conversion experiments. Faster electronics and development of a stiffer and more efficient bending mechanism will allow higher modulation rates, limited by the resonance frequency of the piezo benders (≈ 300-500 Hz). To achieve even faster rates, a different technology should be used for applying perturbations to the fibers, e.g. utilizing all-fiber acousto-optical modulators [38] or the 'smart fibers' technology with integrated modulators [39]. Optical fibers with built-in modulators can also be utilized for a scalable implementation of our method.

Conclusions and Outlook
In this work we proposed a novel technique for controlling light in multimode optical fibres, by modulating its TM using controlled perturbations. We presented proof-of-principle demonstrations of focusing light at the distal end of the fiber, and conversion between guided modes, without utilizing any free-space components. Since our approach to modulate the TM of the fiber is general and not limited to mechanical perturbations, it could be directly transferred to other types of actuators, e.g. in-fiber electro-optical or acousto-optical modulators, to achieve all-fiber, loss-less, fast, and scalable implementations. The all-fiber configuration and the possibility to control more degrees of freedom than the number of guided modes, makes our method attractive for fiberbased applications that require control over multiple inputs and/or wavelengths. Moreover, the possibility to achieve high dimension complex operations opens the way to the implementation of optical neural networks. Our system can provide an important building block for linear reconfigurable transformations, which can be further used in combination with fibers and lasers that exhibit strong gain and/or nonlinearity for deep learning applications.

Response Time
To measure the typical response time of the system, we introduced abrupt changes to the voltages applied to a subset of the piezoelectric actuators, and recorded the speckle pattern obtained at the distal end of the fiber. We then calculated the 2D Pearson correlation coefficient between each of the captured frames and the first frame. The measurements were repeated using different subsets of piezos. Examples of a few of these measurements, for subsets that include between one and four actuators, are shown in Fig. S1(a). The abrupt voltage change causes a fast change to the recorded speckle pattern, yielding a sharp decrease in the computed correlation coefficient. As expected, the bigger the subset of the piezos, the stronger the correlation drop. This sharp decrease is the result of the change in the actuators configuration (the bend they pose), and manifests in a change to the captured speckle pattern. Once the actuator's position stabilized, the correlation stabilized on a lower value. To ensure that the patterns with lower correlation with regard to the first frame are correlated with one another (thus ensuring that the plateau is not a result of the statistical properties of speckles), we also calculated the 2D correlation coefficient of each frame from the last acquired frame. These results are shown in Fig. S1(b) for the same groups of actuators. The high correlation after the configuration change indeed verifies that the speckle pattern did not change further. Based on such measurements we were able to estimate the response time of the system at 30 ms, which corresponds to modulation rates of 33 Hz.

Decorrelation Time
To estimate the stability of the system, we calculated the 2D correlation coefficient of the speckle pattern at the distal end of the fiber over time when the system is idle, i.e. no changes are performed to the states of the actuators. This loss of correlation is known to be attributed to the sensitivity of bare optical fibers to thermal fluctuations in the room and changes of pressure due to air flow. With the GRIN MMF, we found that the system remained highly correlated (corr ≥ 0.99) for 10 minutes. The correlation decreased slowly and linearly for 55 minutes, reaching corr = 0.976. The correlation then decreased faster, reaching corr = 0.883 after two hours. With the SI FMF, the system remained stable and highly correlated (corr ≥ 0.996) over the course of 15 hours.

Rayleigh Statistics
The slopes of the linear scaling of the focusing enhancement factor η as a function of the number of degrees of control rely on the intensity statistics of the generated speckle patterns. The theoretical values reported in the main text are derived for Rayleigh intensity statistics [1]. It is therefore important to compare the intensity statistics of the speckle patterns we obtain in our system with the predictions of Rayleigh statistics. Such a comparison is depicted in S2, which shows excellent agreement with theory.

Optimization Technique
As described in the main text, the results were obtained by finding solutions to optimization problems. These problems used a feedback loop-at each iteration, the speckle pattern at the distal end of the fiber was recorded using the CMOS camera. This pattern was evaluated according to its similarity to a target pattern, and this score was given to the optimization algorithm as a cost, which it tried to minimize by changing the configuration of bends which are applied to the fiber segments. Lower costs were obtained for bend configurations which yielded patterns with high similarity to the target.
The optimization algorithm we chose to use is the Particle Swarm Optimization (PSO), which is a genetic algorithm. It randomly initializes a population of points (referred to as particles) in an M -dimensional search space, representing the voltages which are assigned to the M actuators. These positions are iteratively improved according to their local and global memory from previous iterations. Its stochastic nature helps avoiding local extrema in non-convex problems. An open-source implementation of PSO [2] was modified to fit our experimental setup and simulation. We defined a single run as a single instance of the optimization process, i.e. achieving a single optimized target speckle pattern, such as the example shown in Fig. 3(b) in the main text. With the GRIN MMF, each such run constituted of 80 iterations, with the following hyper-parameters: population size of 120, inertia weight of w = 1, inertia damping ratio of w damp = 0.99, personal learning coefficient of c 1 = 1.5, global learning coefficient of c 2 = 2. With the SI FMF, each run used 86-108 iterations, with a population of size 50. The values of the other hyper-parameters were not changed.

Simulation
Since our system is linear in the optical field, it is natural to describe the propagation of light in it with matrix formalism. We divided the fiber into multiple segments, calculated the transmission matrix (TM) of each segment and computed the total TM of the fiber by multiplying them. To represent our experimental system, we composed bent segments (which mimic the effect of actuators) and straight segments (for the propagation between actuators). A bent segment was approximated by a circular arc, with a defined curvature. To find the guided Figure S1: Response time of the experimental system. The 2D correlation coefficient of each frame with (a) the first and (b) the last of the acquired frames, when a configuration of actuators (the voltage which is applied on these piezos) is changed. Blue lines show a change of configuration of a single actuator, red of two actuators, yellow of three and purple of four. Figure S2: Intensity distribution of the speckle patterns at the end of the fiber. The natural logarithm of the probability distribution function (PDF) of the speckle patterns intensities, as a function of normalized pixel intensity (dots), and a linear fit (line). The data points correspond to experimental readouts, without background noise subtraction. modes and propagation constants of different segments, we used a numerical module [3] which solves the scalar Helmholtz equation under the weakly guiding approximation [4]. We used 10 radii of curvatures, to simulate 10 different vertical positions of the actuators, which impose 10 different perturbations. These radii were linearly spaced between a maximal and a minimal value, which we estimated from the experimental system.
Mode-mixing in short GRIN fibers mostly occurs within groups of degenerate modes. To mimic this phenomenon, we introduced unitary block matrices, whose block sizes were determined according to the mode degeneracy, as expressed in the propagation constants, allowing mixing between modes with the same propagation constants. It is note worthy that without introducing this feature, we were unable to achieve focusing.
We used the same discrete set of possible curvatures for all of the actuators in all runs, and the same optimization mechanism as the experimental setup to achieve a focus. The optimization process mapped one of the possible curvatures for each of the bent fiber segments. In runs where not all of the actuators where used, the remaining were set as straight segments (with no curvature) to maintain the same propagation distance in all runs. Fig. S3 shows the enhancement factor η as a function of the number of actuators whose curvatures were optimized for simulated fibers. It is noticeable that the enhancement factor scales linearly with the number of simulated actuators, where the slope ranges between 0.57-0.64 for the displayed fiber parameters, at a wavelength of λ = 632.8 nm. Figure S3: Simulation of focusing in a multimode optical fiber. The average enhancement factor that was achieved (circles) and the standard error (bars), as a function of the number of actuators whose curvatures where modified as part of the optimization process in the simulation. Several results obtained for different fiber parameters are shown in different colors, along with a linear curve (gray dashed line).