Chromato-axial memory effect in step index multimode fibers

Multimode fibers (MMF) are used in many applications from telecomunications to minimally invasive micro-endoscopic imaging. However, the numerous modes and their coupling make light-beam control and imaging a delicate task. To circumvent this difficulty, recent methods exploit priors about the transmission of the system, such as the so-called optical memory effect. Here, we quantitatively characterize a chromato-axial memory effect in step-index MMF, characterized through its slope $\delta z/\delta \lambda$ and its spectral and axial widths. We propose a theoretical model and numerical simulations in good agreement with experimental observations.

Multimode fibers (MMF) are used in many applications from telecomunications to minimally invasive micro-endoscopic imaging. However, the numerous modes and their coupling make lightbeam control and imaging a delicate task. To circumvent this difficulty, recent methods exploit priors about the transmission of the system, such as the so-called optical memory effect. Here, we quantitatively characterize a chromato-axial memory effect in step-index MMF, characterized through its slope δz/δλ and its spectral and axial widths. We propose a theoretical model and numerical simulations in good agreement with experimental observations. Multimode fibers (MMF) are ubiquitous tools having a key role in telecommunications [1], and driving research in many other related fields such as MMF lasers [2], fiber-based tunable optical cavities [3] and reconfigurable linear operators in quantum photonics applications [4,5]. MMF have also been shown to be useful as high-resolution spectrometers [6,7]. In addition, MMF have recently raised major hopes for developing noninvasive imaging techniques and minimally invasive endoscopes [8,9]. Due to the complex nature of multi-mode guided propagation, light field control through MMF has highly benefited from recent progresses in wavefront shaping techniques, initially developed for adaptive optics in astronomy [10,11] and then extended to complex media [12]. As for propagation through complex media, uncontrolled coherent light propagation through a MMF results in a speckle intensity pattern [13,14]. Light-field control at the distal tip of the fiber demands iterative optimisation [15] or transmission matrix (TM) measurement [12,16]. Not only this calibration is experimentally delicate, but it is also very sensitive to perturbations such as mechanical or thermal fluctuations. In this regard, significant progresses have been made for controlling short MMF eigenmodes, resulting in particular in an increased robustness to bending [17,18]. New strategies have also been developed to help TM measurement, using continuous optimisation [19] or exploiting priors [20].
A major prior information about the TM of complex media is the so-called "optical memory effect" (ME), referring to optical-field transforms which commute with the TM, or equivalently, which are diagonal in the eigenbasis of the TM [20]. The expression ME has been coined in the context of scattering media where for a thin enough diffuser, wavefront tilting has been demonstrated, both theoretically [21] and experimentally [22], to be preserved along propagation. This spatial ME has then been extended to scattering media with strong anisotropy factors [23,24]. Based on the cylindrical symmetry of MMF, rotational ME could be demonstrated [25,26]. Recently, the concept of ME has even been broadened to the spec-tral domain [27,28]. It has been shown experimentally that a chromatic shift induces an axial drift of a beam focused by wavefront shaping behind 1 mm-thick brain tissues, over spectral widths as large as ∼ 100 nm [27]. Theoretical modeling established that this broadband chromato-axial (χ-axial) ME is a characteristic of forward scattering media thinner than the transport mean free path [28], where the product λz (wavelength, axial plane), is specifically conserved when illuminated by a plane wave [28,29]. In essence, the χ-axial ME is due to the conservation of the transverse component of the wave-vector under spectral detuning. Interestingly, the transverse wave-vector conservation naturally holds in step-index (SI) MMF [16]. Furthermore, broadband light propagation through a MMF has recently been observed resulting in an axially extended focus allowing efficient volumetric imaging [30,31]. In this letter we experimentally characterize the χ-axial ME, so far unexplored, at the distal facets of SI-MMFs, and theoretically model it based on the study of the correlation function. The two eigen-axes of the correlation function, diagonal wherein the TM of the MMFs, are analytically derived, together with their respective widths.
As illustrated in Fig. 1a our experiment consists in a tunable mode-locked Ti:Sapphire laser emitting a linearly polarized monochromatic light coupled to a SI-MMF with a microscope objective. In all experiments, MMFs of 0.22 numerical aperture (NA) and of lengths L = 29 mm, 58 mm and 120 mm are tightly hold straight. The objectives' NA (0.25) are chosen slightly larger than the fibers ones to ensure that light is well coupled to all fiber modes. At the MMF output an objective lens mounted on a axial-translation stage collects the light and sends it to a CCD camera through a tube lens. This arrangement enables imaging different axial planes at the MMF output, while tuning the wavelength in the range 800 ± 5 nm, as illustrated in Fig. 1b. The longitudinal axis is denoted as z and the position z = 0 is defined as the MMF output facet. The input wavefront can be controlled by a spatial light modulator (SLM) conjugated with the back focal plane of the input objective in order to either refocus the beam after the MMF or generate random input fields.
In a first set of experiments, we engineered a focused 800 nm laser spot just after the MMF by input wavefront modulation, using the TM of the fiber. The TM is measured by sequentially displaying a basis of modes onto the SLM and gradually phase stepping them [12]. The focus is then generated by phase conjugation. When applying a spectral shift to the input beam, the focus is observed to be axially shifted at the MMF output (see Fig. 1c), similarly to what was observed through forward scattering media [27]. Noteworthy, over the explored spectral range, the contribution of the dispersion of the SLM and other optics to this effect is only a few percent [27] as discussed in the Supplemental Materials [32]. The χaxial ME is not only observed for an engineered focused spot, but also when coupling random light patterns to the fiber. In this case, the similarity of speckles can be observed by computing the zero-mean cross-correlation products between output speckles [28]: (1) withĨ = I − I , and where stands for spatial averaging over all speckle grains. The 0 index represents the reference speckle image at z = 0 and λ 0 = 800 nm, to which all the other ones are compared. Importantly, before computing the correlation product of experimental images we numerically corrected them from potential transverse drifts and intensity inhomogeneities [32]. The correlation products of so-corrected experimental intensity patterns are shown in Fig. 1d as a function of δz and δλ.
Similarly to what we obtain for a focused spot, a χaxial ME is obtained for output speckles: a correlation line appears when comparing speckle, proving the coupling between λ and z shifts. In addition, these plots illustrate that the χ-axial ME extends over limited spectral and axial ranges and that the longer the fiber the steeper the correlation slope (see . Noteworthy, the spectral and axial correlation ranges associated with the χ-axial coupling are both larger than the ones usually defined for fixed λ (Rayleigh length: l z = 2λ NA 2 [33]) and z (l λ = 2n1λ 2 LNA 2 [7, 34]). l z and l λ are illustrated in Fig. 1d and their experimental values for our fibers are given and discussed in [32]. We thus now present a simple model allowing us to discuss how the correlation slope scales with physical parameters in the ray-optics approximation.
Since the fiber diameters (50 µm and 105 µm) are large in comparison to the working wavelength λ 0 = 800 nm, we first consider a ray optics approximation model [35]. In a SI-fiber, the modulus of the transverse component of

Moveable
Step-Index NA = 0.22 core Ø: 50µm, 105µm FIG. 1. Experimental setup and procedure (a) A monochromatic tunable laser illuminates a MMF. The output filed at different axial positions is images on a CCD camera. A reflective phase-only SLM shapes the input wavefront (here represented as a transmission one for schematic convenience). (b) For each wavelength a microscope objective axial scan is performed and a set of images is acquired. (c) Observation of χ-axial effect. A TM is measured at λ0 = 800 nm and z0 = 0 and a focusing phase mask is applied on the SLM (central panel). The focus remains varying both the wavelength and the axial position. (d) Correlation plot of speckles acquired for the same fiber as in (c). A high-valued antidiagonal correlation region is visible. The colour dots make the correspondence with the images of (c). The widths of the profiles along the z axis (lz) and along the λ axis (l λ ) are presented on the side. the wavevector is conserved, and the optical path length along a multiply bouncing trajectory scales quadratically with the incident angle [26], similarly to the phase accumulated under free-space propagation in the Fresnel approximation. Fresnel diffraction formula implies that the intensity is unchanged when the product λz, appearing as the prefactor of the quadratic phase, is conserved, so resulting in the χ-axial effect [28]. In a SI-MMF, taking the differential of the equation λz = const. results in the following slope for the χ-axial ME at the fiber output: where L is the fiber length, and where n 1 , the refractive index contrast of the core with the output medium, is due to Snell-Descartes's law applying at the output facet of the MMF (see [32] for a detailed derivation of Eq. (2)). A spectral shift of the impinging beam thus results in a homothetic axial dilation of the intensity in the SI-MMF, with the origin of the homothety at the input facet of the fiber. This is noticeably different from the χ-axial ME in forward scattering slabs [27] wherein the origin of the dilation was found to be in a virtual plane located at 1/3 of the slab thickness [28,29]. The slope obtained from Eq. (2) is plotted in Fig. 2d, showing good agreement with the corresponding experimental measurements. Noteworthy, within the framework of the ray optics model, the slope does not depend on the fiber radius, in qualitative agreement with experimental observations for the 50 µm and the 105 µm fibers that exhibit similar slopes, since this model is justified in both cases. However, even though the ray model well predicts the measured slopes, it does not explain why the effect remains limited to finite χ-axial ranges, as observed experimentally ( Fig. 1d and Figs. 2a-c). A wave model is thus required to understand and estimate the limited spectral correlation width ∆λ and the limited axial range scanning ability ∆z that we observe in experiments. In the following, we thus propose to analyse our results in the framework of the SI-MMF eigen-modes, assuming cylindrical symmetry for the fibers.
Here, we express the fibers modes in terms of linearly polarized (LP) modes, to then calculate an analytical expression for the cross correlation product (Eq. (1)) as a function of z and δλ. The slope of the χ-axial effect as well as its spectral and axial bandwidths arise from this expression. Importantly, all calculations are carried out inside the fiber and do not consider free-space propagation outside, making these results not specific to our experimental system but revealing intrinsic spectro-axial properties of SI-MMF. Also, for the sake of simplicity of expressions, spectral components of the field are now described by their angular frequencies ω = 2πc/λ, with c the light speed in vacuum, rather than by their wavelength λ.
SI-MMF propagation eigenmodes are, in the weak guidance approximation, LP modes [36]. The field propagating in a fiber of radius a can then be expressed as: where l and m are the azimuthal and the radial number of the LP-modes, respectively, J l the Bessel function of the first kind of order l, and where β l,m and u l,m /a, the longitudinal and the transverse wavenumbers, are imposed by continuity equations at the core-cladding boundary.
Assuming that the fields in the fiber are random patterns with Gaussian statistics [37], the cross-correlation product used experimentally (Eq. (1)), may be written as: up to a normalization factor. As discussed in [32], intensity correlations at the fiber output are dominated by the phase delays accumulated along the fiber due to the chromatic dependence of β l,m : where, assuming ergodic hypothesis, statistical averaging is replaced by modal averaging over l and m values. Eq. (5) then appears as the characteristic function of variables δ [β l,m (ω)z] = β l,m (ω)z − β l,m (ω )z , which we approximate to the Normal one, hence giving: where the variance is thus calculated over all possible l and m values. The axial shift δz that maximizes the correlation coefficient (i.e. that minimizes the variance) for a given spectral shift δω is thus an extremum of C obtained by solving: where is a positive constant scaling as 1/v, with v = ωaNA/c, the normalized frequency. For large fiber cores (ωa/c 1), thus vanishes. Numerical simulation further reveals that the product v does not strongly depend either of NA or of the core diameter. On average v is found to be numerically on the order of 4.7, both for a 50 µm or a 105 µm core diameter [32]. In the large cores limit, Eq. (8) thus leads to Eq. (2) obtained in the ray optic modeel. The axial shift analytically found in Eq. (8) in the frame of LP-mode modeling is compared to the ones obtained experimentally and theoretically in the ray-optics framework in Fig. 2d.. To compare Eq. (8) to the experimental data, the slope value in Eq. (8) is divided by n 1 to take into account Snell-Descartes relation when exiting the fiber.
We now consider the correlation width of the χ-axial ME along the correlation line observed experimentally . Experimental spectral widths and axial displacement amplitudes are measured full width at half maximum (FWHM) by projecting maximal values of the cross-correlation function on the λ and z axes (see Fig. 3 top). The spectral correlation widths of fibers are observed to depend on their lengths. Conversely, the axial correlation width remains almost unchanged in relation to the fiber length. We measure spectral widths of 7 nm, 3 nm, 1 nm for the fibers of lengths 29 mm, 58 mm and 120 mm, respectively, and a mean value of 130 µm for the axial correlation width.
The analytical expression of the correlation width of the χ-axial ME along the correlation line may be obtained by plugging Eq. (8) into the expression of the correlation function Eq. (6). Replacing δz by δω in Eq. (6) according to Eq. (8), we get: where α is a prefactor close to one on average and whose exact value is discussed in [32]. The spectral width results are presented Fig. 3a (bottom), where we plotted together the analytical formula Eq. (9) (dashed lines) and the results of numerical simulations (solid lines). As expected and experimentally observed the spectral width depends on the fiber length. From the analytical LP model, we calculate expected spectral widths of 20 nm, 10 nm and 5 nm, for fibers of lengths 29 mm, 58 mm and 120 mm, respectively. These theoretical values are thus larger than the experimental ones by a factor ∼ 3 for the shortest fiber, up to 5 for the longest fiber but both the orders of magnitude and the global evolution trend with fiber lengths agree. Making use of Eq. (8), the correlation product can also be expressed as a function of the axial shift: with the axial correlation range which, in agreement with experimental observations, does not depend on the fiber length. The simulated axial width and the Gaussian analytical prediction are shown in Fig. 3b (bottom). The axial correlation range obtained with the LP model is 700 µm inside the fiber, yielding 480 µm outside the fiber, once Snell-Descartes law is taken into account. As for spectral correlation width, the theoretical value is a factor 3 above the one obtained experimentally (130 µm). In Fig. 3 it appears that the theoretical model matches the numerical simulations, except for the largest values of δλ and the smallest values of δz. We attribute this discrepancy to the limit of validity of our modeling hypothesis.
A detailed analysis of the contribution of the different modes to the total χ-axial ME presented in [32] reveals the key role of the less confined modes. Our model assumes that the transverse profiles of LP modes is achromatic, but this approximation becomes inaccurate when the spectral correlation width is too large, as implied for short fibers; hence the better matching of the model with numerical simulations for the longest fiber than for the shortest one in Fig. 3. The spectral sensitivity of the less confined modes also suggests that it is possible to broaden the spectral correlation width of the χ-axial ME by limiting the light coupling to the lowest NA-modes, or by filtering high-NA modes at the output (i.e. working at a numerical aperture slightly below the one of the fiber, still allowing high-resolution imaging if a in-purpose higher-NA fiber is used). Another way to improve the results could be to better fulfill the achromatic input field hypothesis by imaging the SLM in the fiber proximal end instead of imaging it on the microscope objective back aperture.
Finally, for applications where the χ-axial ME is of interest after exiting the fiber, the intrinsic correlation width of the fiber is not the only limitation to consider: a mere geometric limitation arises. Indeed when imaging planes far away, outside the fiber emission cone (defined by the fiber radius and its numerical aperture), the maximum spatial frequency decreases resulting in larger patterns that no-longer correlates with fields nearby the fiber output facet. This geometrical effect thus limits the axial scan to δz c = 2a NA 230 µm that translates into a spectral width limitation adding to the intrinsic correlation range of the fiber. Importantly for practical applications, both limitations are of same the order of magnitude: the FWHM of the intrinsic LP axial width is (using Eq. (11)) 2 ln (2)∆z 2.5δz c . None of the mentioned phenomena should then dominate.
In conclusion, we experimentally demonstrated and quantified the χ-axial ME in a SI-MMF. Two different theoretical approaches enable us to predict the value of the χ-axial shift. A simple ray optics model gives access to the shift only using easily accessible fiber parameters with a good degree of approximation. However, in the frame of this model, an infinite spectral width is predicted and slight discrepancies arise on the slope of the χ-axial ME when the field penetration inside the fiber cladding cannot be neglected. The framework of LP modes not only brings more accuracy on the slope of the χ-axial ME but also allows the derivation of its correlation width. We observe that studied χ-axial scanning ranges are not infinite and also obtain that it is possible to optimize the effect: analytical results demonstrate that working with the most confined modes would allow extending significantly the spectro-axial correlation range of this memory effect. For imaging purposes, when a large spectral bandwidth is required, working with short MMF should be preferred. Alternatively, avoiding the excitation of the less confined modes is expected to increase the spectral correlation range, in principle up to arbitrarily large values.
Characterization of the χ-axial ME in MMFs is of high importance to simplify light-beam manipulation by blind control of polychromatic wavefields at distal facets, especially for imaging applications. Thanks to the χ-axial ME, only one TM measurement would be needed to perform foci on different z planes, and thereby overcome the current inability of wavelength tuning while imaging [30]. Furthermore, the χ-axial ME knowledge enables a priori wavefront correction to achieve non-linear microscopy [38,39], and to extend the confocal microscopy technique of [40] for objects subjected to inelastic scattering or broadband fluorescence. This fast axial scan ability opens up the possibility of extending the "spot scanning" imaging technique in three dimensions [13], paving the way to non-invasive imaging (in biological media for instance).

ACKNOWLEDGMENTS
We would like to thank Michał Dąbrowski and Lorenzo Valzania for the careful reading of the manuscript and fruitful comments. This project was funding by the European Research Council under the grant agreement No. 724473 (SMARTIES). SG and MG are members of the Institut Universitaire de France.

MATERIAL
In the experimental setup presented Fig. 1a of the main text, the monochromatic tunable laser is a MaiTai HP (Spectra Physics). The CCD camera is a Basler ace (acA1300-30uc) and the SLM is from Meadowlarks (HSP512L-1064). Also the MMFs used for the experiments are purchased from Thorlabs (FG050LGA/FG105LCA).

SPECTRAL WIDTHS AT FIXED λ AND z
The axial correlation extent l z of a monochromatic speckle equals the Rayleigh length: l z = 2λ NA 2 [33]. From our experimental measurements presented in Fig. 1d and Figs. 2a-c of the main text, we obtain l z 31 µm FWHM, in agreement with the Rayleigh length expression (l z = 33 µm for a 0.22-NA fiber at 800 nm). Correspondingly, the spectral width l λ in a given transverse plane, is l λ = 2n1λ 2 LNA 2 [7, 34]. We measured spectral widths (FWHM) equal to 1.4, 0.8 and 0.4 nm for L = 29, 58 and 120 mm, respectively, in agreement with the analytical expression of l λ . Both l z and l λ lengths are represented in Fig. 1d.

EXPERIMENTAL PRECAUTIONS AND DATA PROCESSING
Alignment post-correction: To prevent any correlation drop due to unavoidable alignment issues during the z-scan we used a reference. For each experiment we performed a z-scan with a white light beam and tracked its center. The center displacement enabled us, by cropping the speckle data, to correctly correlate the images and minimize any misalignment induced correlation drop.
Correlation measurement: Before computing the correlations we selected on the data the central region with approximately 200 speckle grains. In this region we used the reference white light z-scan to measure an intensity profile for each axial position (see Fig. S1a). Indeed if the profile is flat for z = 0, imaging planes away from this position leads to intensities peaked at the centre. We then used the measured reference intensity profiles to correct the speckle images. The reason of these two precautions is the following: if, in addition to the speckle short range variations, long range variations (due to the spatial illumination non-flatness) are present, they impact the correlation value. One hence may observe a -relatively important -correlation background and the correlation value does not reach zero for uncorrelated speckles. It is however important to point out that the χ-axial effect is still well visible in absence of correction and that the later only affect the correlation background.
Speckle dilation and z-scan: For planes away from the fiber output (z = 0), the imaged speckle is dilated due to Fresnel diffraction. To evaluate the dilation we measured the speckle grain size by performing the autocorrelation of the intensity image and extracting its FWHM. The Fig. S1b shows the speckle grain size evolution along z for the L = 58 mm fiber. For each z position the grain size value results from the average for all wavelengths. For this experiment the relative variation of grain size is of R 58 = 0.45. The mean value for all three main experiments is R = 0.42. However because of free space propagation the speckle transverse dilatation on the axial scan is compensated by the wavelength detuning such that on the high correlation line of the effect no dilatation is observed as visible on Fig. S1c.

SLM impact:
We also checked that the chromaticity of the SLM had no impact on the experimental results by suppressing it on one experiment. It did not affect the χ-axial ME effect itself neither its quantitative value.

COMPARISON SI-MMF, GI-MMF AND LENS CHROMATICITY
To ensure that the χ-axial shift is not a trivial effect we evaluated the shift one could expect from the mere chromaticity of the experimental setup optical components. The microscope objectives used are treated against chromaticity in the visible range. At 800 nm it remains δz δλ ≈ 0.4 µm nm −1 . The lens used is plane-convex, made of N.BK7 glass. The chromaticity evaluated at 800 nm gives δz δλ ≈ 0.1 µm nm −1 . These values are far different from the effect observed both in order of magnitude and in sign. However it matches the slope obtained when performing the experiment with a GI fiber (see Fig. S2).
FIG. S1. (a) Evolution of the white light intensity profile along the z axis for a cut on the y dimension for the 58 mm long SI fiber. The bottom gray mark represents the speckle central position kept to compute the correlations. (b) Evolution of the speckle grain size along the z axis. The solid blue line represents the mean grain size (mean over all wavelength shifts) for the L = 58 mm fiber. The relative variation is R58 = 0.45. The red (resp. yellow) dotted lines is the mean grain size evolution for L = 29 mm (resp. 120 mm) fiber giving a relative variation of R29 = 0.34 (resp. R120 = 0.45). The mean relative variation for these three experiments is R = 0.42. (c) Speckles of raw images along the χ-axial high correlation line. No important transverse dilatation is observed.

CALCULATIONS FOR RAY OPTICS MODEL
To derive the χ-shift in the framework of ray optics, let us focus on the electric field on different planes perpendicular to the optical axis when the incident wavelength varies. We make the following assumptions: E ω (x, y, −L) does not depend on the angular frequency ω (save for a global phase factor); on the whole range of explored wavelengths, the refraction indices of the fiber (n 1 and n 2 ) do not depend on ω; the angles between optical rays and z-axis are small such that the paraxial approximation holds; with a step-index fiber, k 2 ⊥ and k z are conserved throughout the propagation inside the fiber, hence their input and output values are equal. Hence the light rays follow the same path inside the fiber whatever the angular frequency: varying ω just implements an homothetic transformation of the wave-vector distribution, enabling to write: where ∆φ = ∆φ(k x , k y , ω, ω 0 ) = φ ω − φ ω0 with φ the phase accumulated while propagating inside the fiber, andẼ indicates the spatial Fourier transform of E. Now expressing the electric field at a position z for an angular frequency ω gives where k z is derived from (ω, k x , k y ) through the dispersion relation : Using Eq. (S1) we obtain One needs to express ∆φ. We assume (with a posteriori confirmation) that ∆φ(k x , k y , ω, ω 0 ) = ∆φ 0 (ω, ω 0 ) + k 2 x +k 2 y (ω/c) 2 ∆ζ(ω, ω 0 ). This is equivalent to show that the phase φ is proportional/has a term in sin 2 θ, where θ is the incident angle of the input light ray. At this step the k ⊥ conservation is mandatory to identify the k 2 ⊥ = k 2 x + k 2 y term in φ (therefore after propagation in the fiber) with the input k 2 ⊥ = ω 2 c 2 sin 2 θ.
As a consequence, one has Thus, save for a spatial transverse magnification by a factor of ω0 ω , both speckle E ω resp. E ω0 are proportional when considered at abscissa z = 2∆ζ ω/c resp. z = 0. It is noteworthy that for λ 0 = 800 nm and ∆λ = 10 − 20 nm, the relative transverse magnification of 10 −2 , 2.10 −2 is hardly observable. If one limits to a short range of ω, Note that a local slope dz dω also exists if ω is far from ω 0 . It is just slightly more cumbersome to write down its analytical value. There is no ω-range limitation for this effect on this model.
To go forward, one needs to evaluate the phase φ 0 accumulated during the light propagation in the fiber. For a pure ray model the accumulated phase during propagation is Using Eq. (S7), one gets As commented in the main text the value of the slope is in global agreement with the experimental data.

Derivation of Eq. (5) of the main text
In this section we present the full derivation of the LP model used in the main text. We first recall the derivation central thread. After expressing field propagation in the fiber in terms of linearly polarized (LP) modes, we calculate an analytical expression for the cross correlation product (Eq. (1)) as a function of z and δλ and deduce from this expression the slope of the χ-axial effect as well as its spectral and axial bandwidths. Along the derivation we point out correspondences between the assumptions used and the ray optics model. Importantly, all the following calculations are carried out inside the fiber and do not consider free-space propagation outside, making these results not specific to our experimental system but revealing intrinsic spectro-axial properties of SI-MMF. For the sake of simplicity of expressions, spectral components of the field are described by their angular frequencies (independent of the refractive index of the fiber-core) ω = 2πc/λ, with c the light velocity in vacuum, rather than by their wavelength λ.
SI-MMF propagation eigenmodes are, in the weak guidance approximation, LP modes [36]. The field propagating in a fiber of radius a can then be expressed as: E ω (r, z) = l,m E l,m (ω)e iβ l,m z e ilϕ J l u l,m r a (S11) where l and m are the azimuthal and the radial number of the LP-modes, respectively, J l the Bessel function of the first kind of order l, and where β l,m and u l,m , the longitudinal and the transverse wavenumbers, are imposed by continuity equations at the core-cladding boundary. In particular, they satisfy the following equation: (S13) Assuming that the fields in the fiber are random patterns with Gaussian statistics [37], the cross-correlation product used experimentally (Eq. (1) of the main text), may be written as: up to a normalization factor. Using Eq. (S11), we get the expression of the two-wavelength mutual coherence function: To calculate analytically Eqs. (S14) and (S15), we need to make a few simplifications. Indeed, since u l,m depend on the wavelength, so do the LP-modes. When sending a field for which both amplitude and phase are assumed to be achromatic inside the fiber (for instance a pinhole illuminated with white light), the propagation of the different wavelengths should be written using different sets of eigenmodes. Alternatively, we shall consider that the transfer matrix from the eigen-basis at ω to the one at ω is not identity and thus, that the projection coefficients, for an achromatic impinging field, depend on the wavelength. However, in practice, the transverse profiles J l u l,m r a only weakly depend on the wavelength because u l,m is close to the m th zero of J l and spectral detuning dependence of u l,m is a second order perturbation. The following assumptions are then made: • Spectral fluctuations of the E l,m coefficients do not significantly contribute to the cross-correlation product (Eq. (S14)), since they just modify the weighting. Equivalently, although Bessel functions are not exactly orthogonal in Eq. (S15), their inner product will only contribute for a weighting change.
• In contrast, intensity correlations at the fiber output are dominated by the phase delays accumulated along the fiber due to the chromatic dependence of β l,m .
The two-wavelength mutual coherence function then becomes, It comes that C is as a weighted average of phasors e i(β l,m z−β l ,m z ) . Assuming that the coefficients E l,m (S17) As stated in the main text, we assume ergodic hypothesis making statistical averaging equivalent to modal averaging over l and m values. Eq. (S17) is the characteristic function of the variable δ [β l,m (ω)z] = β l,m (ω)z − β l,m (ω )z . In order to obtain a closed form expression, we approximate the characteristic function of the modal distribution of δ [β l,m (ω)z] by a Normal one: Differentiation of δ [β l,m (ω)z], relying on Eq. (S13) yields: with Eq. (S19) can then be written as a second order polynomial in ωδz and zδω: where A, B and C are given by: where indices l and m have been dropped for the sake of brevity. The variance in Eq. (S21) is thus the sum of two parabolas. The orthogonal eigen-axes of these parabola can be identified by writing Eq. (S21) in a matrix form: where it appears that eigen-axes are eigen-vectors of the involved 2 × 2 matrix. Eigen-values are trivially: Importantly, from Eqs. (S22), it appears that the determinant AC − B 2 appearing in Eq. (S24) vanishes if χ has a Dirac distribution. In practice, the distribution of χ is highly peaked and deviation from the peak value is due to the less confined modes as illustrated in Fig.S3. Since u scale as v and χ as 1/v (see Eq. (S41)), with v = ωaNA/c, the normalized frequency, it is easy to show that AC − B 2 scales as v 6 while (A + C) 2 scales as v 8 . Consequently, for highly multimode fibers, for which v 1, we may simplify the expressions of eigen-values s ± : resulting in a ratio s+ s− (A+C) 2 AC−B 2 scaling as v 2 1. Relying on the same approximation, the expression of the eigen-vectors are then In the eigen-basis (V + , V − ) of the system, the correlation function in Eq. (S18) can thus be simplified and expressed in the eigen-system of coordinates (ξ + , ξ − ) as: with: As a result, the sought-for correlation (Eq. (S18)) is a product of two Gaussians, one of which being much larger than the other because s + s − , so resulting in a cigar-like shape of the correlation function. It is noteworthy that the χ-axial ME is well established by diagonalizing the correlation function. This approach is similar to [26] and shows that not only diagonalizing the TM brings information on the field transformation but that the speckle correlation function, sometimes quite easily analytically accessible, also does. Indeed the intensity cross-correlation was already shown to allow interferometric measurements [41].

Expression of
The long axis of the cigar-shaped correlation function is a line satisfying ξ + = 0, so yielding Eq. (7) in the main text: where is: From numerical simulation, we observe that does not strongly depend either of NA or of the core diameter (for ωa/c 1) and we obtained that the product v may be considered as close to 4.7 on average, both for a 50 µm or a 105 µm core diameter (see Fig. S4).

Expression of α
Along the axis defined by ξ + = 0, the correlation function Eq. (S29) is: Because of successive approximations carried out in the expression of eigen-vectors and eigen-values, it is simpler to express the correlation coefficient starting from Eq. (S21) rather than from the expression of s − . The correlation function along the axis δξ + = 0 may then be re-written either as a function of δω or δz in a simplified form by plugging Eq.(S32) into Eq. (S21): Where the expression of α introduced in the Eq. (10) in the main text is: We introduced v together with the numerical aperture NA (v = kaNA) in the definition of ∆ω and ∆z in order that α be a dimension-less parameter. Numerical simulations show that α fluctuates between 0.45 and 1.4 as a function of the core diameter, because of abrupt changes in modal cut-offs (as for the product v ). However, assuming modal coupling in the fiber, we may consider that α is equal to unity on average, both for a 50 µm or a 105 µm core diameter as presented in Fig. S4.

SIMULATION OF THE LP MODES
For the simulation of the LP modes we used some functions of the Matlab library of Michael Hughes [42].
For different fibers (Ø = 50 µm or 105 µm), 0.22-NA, n core = n 1 = 1.4533 the modes of the fiber are numerically calculated. The prefactors α and v are determined using respectively the definition linked to Eq. (10) of the main text and Eq. (S34). The value of χ can be obtained by using the following equality: χ l,m = 2ω ∂ ∂ω [ln (u l,m )] = 2 v 2 − (u l−1,m ) 2 + l 2 + 1 , wherein for large enough core diameters (ka 1), u l,m can be replaced by the m th zero of the Bessel function of order l, J l [43]. The distribution of χ l,m given by Eq. (S41) is presented in Fig. S3a for the 50 µm diameter fiber and Fig. S3b for the 105 µm diameter one. It appears that χ l,m is almost uniform for the most confined modes (low l and m values) but diverges nearby the critical angle.
The values of α and v kept for the calculations are the mean values obtained from small variations of the fiber core (due to fabrication imperfections), see Fig. S4.
Fibers modes are then excited by an incident field. The field is propagated along the fiber on a distance L with the propagation constants of the modes previously determined. Performing these steps for different wavelengths and fiber lengths provides the data necessary to mimic the experience.

DISCUSSION OF THE ROLE OF THE LESS CONFINED MODES
For the calculation of the phase variance in Eq. (6) of the main text, differentiation of the product β l,m z for a given mode yields: wherein it appears that the specific phase for the mode (l, m), the second term, is unchanged for a χ-axial shift: which slightly differs from Eq. (8) in that is now replaced by χ l,m . When calculating the phase variance driving the limited correlation range of the χ-axial ME, a fine analysis of χ l,m appears very informative. As presented in Fig. S4, the χ l,m values are almost the same for all modes except for the less confined ones. If exciting only the most confined modes, we then get χ, with χ the common value of these confined modes, and α → +∞. As a result, if not injecting the less confined modes, responsible for field penetration in the fiber cladding, the χ-axial ME extends over an infinite range in full agreement with the ray optics approximation. Interestingly, it therefore means that the limited axial range of the χ-axial ME is dominated by the less confined modes.