# Characterization and Exploitation of the Rotational Memory Effect in Multimode Fibers

Rodrigo Gutiérrez-Cuevas, Arthur Goetschy, Yaron Bromberg, Guy Pelc, Esben Ravn Andresen, Laurent Bigot, Yves Quiquempois, Maroun Bsaibes, Pierre Sillard, Marianne Bigot, et al.

► To cite this version:

Rodrigo Gutiérrez-Cuevas, Arthur Goetschy, Yaron Bromberg, Guy Pelc, Esben Ravn Andresen, et al.. Characterization and Exploitation of the Rotational Memory Effect in Multimode Fibers. Physical Review X, 2024, 14 (3), pp.031046. 10.1103/PhysRevX.14.031046. hal-04813257

HAL Id: hal-04813257

https://hal.science/hal-04813257v1

Submitted on 1 Dec 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Characterization and Exploitation of the Rotational Memory Effect in Multimode Fibers

Rodrigo Gutiérrez-Cuevas,¹ Arthur Goetschy,¹ Yaron Bromberg,² Guy Pelc,² Esben Ravn Andresen,³

Laurent Bigot,³ Yves Quiquempois,³ Maroun Bsaibes,³ Pierre Sillard,⁴ Marianne Bigot,⁴ Ori Katz,²

Julien de Rosny,¹ and Sébastien M. Popoff¹,※

¹ESPCI Paris, PSL Univ., CNRS, Institut Langevin, Paris, France

²Racah Institute of Physics, The Hebrew Univ. of Jerusalem, Jerusalem, Israel

³Univ. Lille, CNRS, PhLAM-Physique des Lasers, Atomes et Molécules, F-59000 Lille, France

⁴Prysmian Group, Parc des Industries Artois Flandres, Douvrin, France

<figure><img src="image_2.png" /></figure>

(Received 30 October 2023; accepted 8 July 2024; published 16 September 2024)

In an ideal perfectly straight multimode fiber with a circular core, the symmetry ensures that rotating the input wave front leads to a corresponding rotation of the output wave front. This invariant property, known as the rotational memory effect (RME), remains independent of the typically unknown output profile. The RME thus offers significant potential for imaging and telecommunication applications. However, in real-life fibers, this effect is degraded by intrinsic imperfections and external perturbations, and is challenging to observe because of its acute sensitivity to misalignments and aberrations in the optical setup. Building on a previously established method for precisely estimating fiber transmission properties, we demonstrate an accurate extraction of RME properties. Additionally, we introduce a comprehensive theoretical framework for both qualitative and quantitative analysis, which specifically links the angular-dependent correlation of the RME to the core deformation's geometrical properties and the fiber's mode characteristics. Our theoretical predictions align well with experimental data and simulations for various amounts of fiber distorsion. Finally, we demonstrate the ability to engineer wave fronts with significantly enhanced correlation across all rotation angles. This work enables accurate characterization of distributed disorder from the fabrication process and facilitates calibration-free imaging in multimode fibers.

DOI: 10.1103/PhysRevX.14.031046

Subject Areas: Optics, Photonics

I. INTRODUCTION

Optical fibers present a unique opportunity for mini-
mally invasive imaging deep within the human body. Most
flexible medical endoscopes utilize multicore fibers or fiber
bundles [1]. Comparatively, multimode fibers (MMFs)
offer orders of magnitude higher information density,
allowing, in theory, an increase in image resolution or a
decrease in the endoscope footprint [2]. However,
dispersion distorts the input image, a phenomenon that
is exacerbated by mode coupling introduced by defects or
deformations within the fiber. For this reason, image
reconstruction techniques through multimode fibers hinge
on estimating [3] or measuring the transmission matrix
(TM) [4–8], i.e., the relationship between the input and
output fields of the optical system. Unfortunately, this TM
approach is prone to real-time changes due to dynamic fiber

bending and temperature fluctuations [9], which prevent
the direct use of a previously calibrated system.

Similar challenges are encountered in utilizing MMFs
for telecommunications, where they hold the potential to
significantly boost data rates compared to their single-mode
counterparts. Through mode-division multiplexing, it is in
principle possible to utilize different fiber modes as
independent channels, effectively multiplying data rates
by the number of modes employed, without substantially
increasing cost or footprint. However, the occurrence of
mode coupling, even in fibers with a limited number of
modes, i.e., MMFs with fewer than 10 modes, currently
hinders their application in long-haul communications.
Enhancing fiber design necessitates a deep understanding
of the disorder-induced effects that lead to crosstalk,
including defects arising from the fiber drawing process.
Addressing these challenges remains difficult, and numeri-
cal models frequently used in fiber design often overlook
these crucial factors [10].

Various fabrication techniques are employed based on the
fiber type and manufacturer, including modified chemical
vapor deposition, vapor axial deposition, outside vapor
deposition, and plasma-activated chemical vapor deposition
(PCVD). One well-recognized challenge across these tech-
niques, attributed to their inherent limitations in precision,

Published by the American Physical Society under the terms of
the Creative Commons Attribution 4.0 International license.
Further distribution of this work must maintain attribution to
the author(s) and the published article's title, journal citation,
and DOI.

<footnote>*Contact author: sebastien.popoff@espci.psl.eu</footnote>
lies in achieving a highly accurate radial index profile, particularly when dealing with sharp index changes. Imperfections in the fiber lead to a mixing of the information carried by the input field. This effect is analogous to the randomization of information when attempting to recover an image through a scattering medium [11,12]. To circumvent the necessity of measuring the TM, which implies an invasive procedure, an elegant solution for imaging is to exploit the invariant properties of the medium, specifically, the angular memory effect [13,14]. For a given illumination, even though the output random pattern remains unknown, the angular memory effect facilitates the shifting of this output speckle pattern in two directions with minimal to no deformation. Tilting the input wave front then allows scanning the object plane with the unknown pattern. Recording the reflected or fluorescent signal provides sufficient information to reconstruct the image of the hidden object [13]. While the range of such an effect is constrained, strategies have been proposed to recover images of objects beyond this limitation [15], making the memory effect highly attractive for noninvasive imaging applications.

Building on decades of research in scattering media, recent interest has surged in the study of coherent effects in disordered fibers [16,17]. Specifically, the TM approach [18–21] and random matrix theory [22,23] have emerged as particularly useful frameworks for these investigations. In particular, a close analogous effect to the angular memory effect in scattering media is observed in the special case of square-core fibers, where a translation of the input wave front results in a corresponding translation at the output [24], albeit with the noticeable presence of artifacts, which can nonetheless be exploited to recover images [25,26]. In more typical cylindrical-core fibers, a similar phenomenon, known as the rotational memory effect (RME), has been recently identified [27,28]. This effect is characterized by the rotation of an input wave front along the optical axis of the fiber leading to a corresponding rotation of the output pattern, even though the latter is unknown. In principle, this effect could be harnessed for imaging through a multimode fiber for which the TM has not been previously measured.

Nevertheless, since its initial observation, no prediction or quantitative description of the RME has been presented. Neither the angular range covered by the RME, nor its dependence on disorder, nor its potential robustness and modularity have been studied or elucidated. Furthermore, the manifestation of an angular revival effect, leading to secondary peaks in the correlation of the output pattern at the rotation angle $\pi$, has been observed but also remains unexplained. An important consideration is that the measurement of the RME is complicated by its high sensitivity to misalignments and aberrations in the optical system used to inject light into the fiber [27]. However, these adverse effects can be understood and compensated using a framework that some of us have introduced in Ref. [29].

The procedure involves learning the input and output aberrations by optimizing a model-based numerical model. This approach enables the retrieval of an accurate TM of the system, even when using imperfect measurements. Additionally, it provides the transformation needed to physically compensate for the input aberrations, which can be directly implemented using a spatial light modulator (SLM). The numerical compensation for aberrations is a crucial step, as it enables precise observation of the RME, which would otherwise be rapidly obscured by aberration effects.

In the present article, we first experimentally extract the RME properties with high accuracy building on a previously introduced method for precisely estimating the transmission matrix of MMFs in the presence of aberrations and misalignments [29]. We then introduce a theoretical framework based on precise disorder modeling, which yields analytical predictions for the shape of the RME correlation function, supported by numerical simulations of the microscopic wave equation in MMFs. Specifically, this framework properly predicts the angular range of the RME as a function of disorder and elucidates the origin of secondary peaks in the correlation function. These are attributed to the subtle interplay between mode symmetry and disturbance symmetry. The theory and simulations faithfully reproduce our experimental results obtained on various commercial MMFs. Finally, based on this analysis, we propose a method to identify input wave fronts that are only weakly sensitive to the symmetry breaking introduced by perturbations, hereby facilitating the exploitation of the RME for calibration-free imaging applications.

## II. MEASURING THE RME

A memory effect is defined in relation to a field transformation. A perfect memory effect exists when the application of this transformation before or after propagation through a given optical system produces the same effect. For RME to occur, the rotation operator $R(\theta)$, with $\theta$ the angle of rotation, must commute with the optical system's transmission matrix $T$ of the fiber [28]. In this study, we consider only one polarization of the field. The matrix $T$ therefore links the input field in a specific circular polarization channel to the output field in the same polarization channel.

### A. Experimental setup and measurement procedure

In principle, the measurement of the memory effect is straightforward, as it simply requires rotating the input wave front and measuring the corresponding output pattern. However, any factor that disrupts the rotational symmetry of the system leads to a degradation of the measurement. In particular, minute misalignments relative to the fiber axis and aberrations can alter the observed results [3,27]. This complicates quantitative characterization by making it impossible to separate the effects of fiber defects from
<figure><img src="image_3.png" /><figcaption>FIG. 1. The rotational memory effect in MMF. When the fiber is illuminated by a coherent wave front $|\psi_{in}\rangle$, a seemingly random transmitted field $T|\psi_{in}\rangle$ is observed at the output. In an ideal MMF with cylindrical symmetry, rotating the input wave front by an angle $\theta$ (i.e., sending $R(\theta)|\psi_{in}\rangle$) and measuring the transmitted field $TR(\theta)|\psi_{in}\rangle$ is equivalent to rotating the output field resulting from the propagation of $|\psi_{in}\rangle$ and measuring $R(\theta)T|\psi_{in}\rangle$. A local perturbation is then added by moving a tip in contact with the fiber over a distance $\Delta x$ transverse to the fiber axis.</figcaption></figure>

those introduced by the optical setup used for measurement. To address this issue, we use the approach we previously introduced in Ref. [29], which allows us to numerically estimate the detrimental effects of optical and experimental imperfections. We then experimentally compensate for these effects by correcting the input wave front impinging on the fiber. This not only ensures accurate, interpretable measurements, but also significantly reduces the time needed to change the sample and carry out the measurement, in just a few minutes.

We first measure the TM of a 24.5-cm segment of a straight 50-µm core graded-index fiber (BendBright OM4 [30]). The choice of fiber studied is guided by the fact that graded-index fibers are the standard for MMFs in telecommunications and are appealing for endoscopic imaging applications due to their relative robustness to bending [31]. Utilizing a fast digital micromirror modulator and an InGaAs camera, we follow the procedure detailed in Ref. [29], which enables us to identify and compensate for aberrations and misalignments. It also allows us to

accurately generate the input masks on the modulator that correspond to rotating the field with respect to the optical axis in the input facet plane of the fiber. The principle of the experiment is depicted in Fig. 1 and is further detailed in Appendix A 1.

## B. Measurement of the RME

To illustrate the effect of the RME, we first observe its impact on a focusing operation. We compute the mask that focuses light at a specific position in the output facet of the fiber using the TM [32,33]. It is noteworthy that the knowledge of the TM is not necessary for this step, nor for any of our measurements, as focusing can be achieved through methods like sequential optimization [34,35] or phase conjugation [6]. We then rotate the input wave front. We show in Fig. 2(a) the sum of the resulting output amplitude patterns for 10 values of the rotation angle. We can see that rotating the input masks allows the focusing spot to be rotated along the optical axis of the fiber, with limited degradation of focusing quality.

To further characterize the RME, we seek to quantify the similarity between a transmitted field $|\psi\rangle = T|\psi_{in}\rangle$ for a normalized input field $|\psi_{in}\rangle$ and the output field $|\psi_{\theta}\rangle = T_{\theta}|\psi_{in}\rangle$, where $T_{\theta} = R(-\theta)TR(\theta)$. The second field corresponds to a rotation of the input and output profiles by an angle of $\theta$ and $-\theta$, respectively. We define a correlation function for this purpose as

$$C(\theta) = \frac{|\langle\psi|\psi_{\theta}\rangle|}{\sqrt{\langle\psi|\psi\rangle}\sqrt{\langle\psi_{\theta}|\psi_{\theta}\rangle}} \quad (1)$$

In practice, we send a set of 100 random input wave fronts, rotate them in the plane of the input facet, and measure the output field. We then compute the average correlation function $\langle C(\theta) \rangle$. Figure 2(b) shows the experimentally measured $\langle C(\theta) \rangle$ for the unperturbed fiber (solid

<figure><img src="image_4.png" /><figcaption>FIG. 2. Experimental measurement of the RME. (a) Rotation of a focal spot. Light is focused at a given output position and the input phase mask is rotated along the axis of the fiber for 10 values of the rotation angle $\theta$. We sum all the resulting output amplitude patterns to reveal conservation of the focal spot, albeit with a variation in intensity, the latter being maximal for angle $\theta = 0$ and $\theta = \pi$. (b) Experimental measurement of the RME angular correlation Eq. (1) as a function of the level of perturbation $\Delta x$. The bright red curve shows results for the unperturbed fiber ($\Delta x = 0$ µm). The red to brown curves correspond to progressively increased disturbances, obtained by applying a local deformation using a translation stage. Results are obtained by averaging over 100 random inputs. Shaded areas correspond to the error estimated by the standard deviation of the experimental data.</figcaption></figure>
red line). A first observation is that the correlation deviates from 1 and that the rotated focal spot is degraded within the $2\pi$ range. This indicates the presence of imperfections in the fiber that break the cylindrical symmetry of the system. We emphasize that the TM measurement is not used to characterize the memory effect; only knowledge of input aberration and misalignment effects is exploited to accurately rotate the input field. Although the results shown here correspond to field correlations, we demonstrate in Appendix A 2 that very similar results are obtained for intensity correlation measurements. The latter can be expressed as $C_l(\theta) \approx C(\theta)^2$. As a result, the behavior of the angular correlation shown in Fig. 2(b) is qualitatively analogous to the result reported in Ref. [27], where the intensity RME correlation was measured. Furthermore, we show in Appendix A 3 that when the TM is known, the RME correlation can be accurately computed without the need for additional measurements.

To qualitatively observe the effect of perturbations on the RME, we gradually apply a controlled deformation to the fiber along an axis orthogonal to the propagation direction. The fiber is maintained on a V groove and we press locally on the fiber from the top with a spherical metallic tip using a motorized translation stage. The correlation $\langle C(\theta) \rangle$ as a function of the rotation angle $\theta$ is presented in Fig. 2(b), for different values of the displacement $\Delta x$ of the tip (red to brown curves).

We first observe that, even without applying a local perturbation to the fiber, the correlation decreases to approximately 60%. When the fiber is held straight, this effect can be attributed to the presence of defects in the fiber. This correlation curve exhibits a second maximum, close to 95%, at $\theta = \pi$, along with small local maxima at $\theta = \pi \pm \pi/2$. These features are indicative of the geometrical defects within the fiber. Upon applying the deformation, we observe that the correlation as a function of $\theta$ decreases globally, and the local maxima vanishes.

### III. MODELING THE RME FOR DEFECT CHARACTERIZATION

The design and manufacture of few-mode and multi-mode fibers with minimal mode coupling for mode division multiplexing in telecommunications pose significant challenges [36]. Achieving industry-compatible levels of crosstalk remains an elusive goal, even for fibers with a low mode count (<10 modes). The application of mode-coupling theory to specific types of disorder enables the prediction of certain adverse effects on telecommunications, such as mode-dependent losses and power mode coupling [37,38]. However, the precise characterization of perturbations in a given fiber is difficult to achieve. The common approach is to assume coupling solely between pairs of neighboring modes. Thus, with the goal to minimize mode-coupling-induced crosstalk, fiber manufacturers focus on solely optimizing the difference in

effective refractive index among the fiber modes [39,40]. By overlooking the interplay between the symmetries of the modes and the perturbations, this method reveals its limitations when attempting to achieve low levels of crosstalk [10], thus restricting the performance of current systems. An accurate disorder model with easily estimable parameters would be a crucial asset for designing multi-mode fibers with low mode coupling in the telecommunications industry.

The acute observed sensitivity of the RME to geometrical deformations, along with its dependence on the fiber's symmetry, indicates that one can derive information about the distributed disorder within MMFs. In this section, we present a theoretical model whose predictions are compared with experimental observations. This model proves capable of predicting all RME behaviors in the presence of the disturbances just described.

#### A. Model of disorder

In an ideal MMF, due to the axisymmetry of the system, perfect RME should be expected; i.e., the rotation of a given input wave front should result in a corresponding rotation of the output wave front. This corresponds to $C(\theta) = 1$ for all $\theta$. However, real fibers are rarely perfect, as demonstrated by the result shown in Fig. 2, resulting in mode coupling that is mainly influenced by the geometrical defects of the fibers [41]. Two main contributions can be identified: large radius bends, attributable to the geometrical conformation of the fiber, and minor distortions at the core-cladding interface, primarily due to fabrication inaccuracies [42–44].

We propose to model fiber disturbances by a deviation of the refractive index profile from a perfect axisymmetric function of the following form:

$$ \delta n(r, \phi, z) = g(z, r) \sum_q \Gamma_q \cos(q\phi + \varphi_q), \quad (2) $$

where $z$, $r$, and $\phi$ are the cylindrical coordinates corresponding, respectively, to the longitudinal (axis of the fiber), radial, and azimuthal directions. The longitudinal variations of $g(z, r)$ are characterized by random fluctuations with a correlation length $l_z$, which is typically of the order of 100 µm [44], while radial variations of $g(z, r)$ are discussed in detail in below. On the other hand, disorder in the azimuthal direction is decomposed into harmonics with orbital momentum $q$ and weight $\Gamma_q$ [23].

We attribute the radial fluctuations to variations between neighboring radial layers, stemming from inaccuracies in the deposition technique or interlayer diffusion of the doping elements. Specifically, in the case of the fiber under investigation, these inaccuracies are associated with the PCVD process. This leads us to approximate the fiber of length $L$ by a succession of $N_z = L/l_z$ segments, each of length $l_z$, in which the perturbation
term $\delta n$ is invariant along $z$. Specifically, for the $p$th segment in the interval $z \in [p l_z, (p+1)l_z]$, we write $\delta n_p(r, \phi) = g_p(r) \sum_q \Gamma_q \cos(q\phi + \phi_q)$. We model $g_p(r)$ as a Gaussian random variable with zero mean, characterized by a standard deviation $\sigma_g(r) = d_{\text{layer}}|dn_0(r)/dr|$ [45], where $n_0(r)$ is the radial profile of the unperturbed fiber, and $d_{\text{layer}} \approx 10$ nm is the typical thickness of each layer formed in the PCVD process [46]. For a gradient index fiber with a parabolic index profile (see Appendix B and Ref. [47]), the standard deviation can be put in the form

$$ \sigma_g(r) \approx \frac{r}{a^2} \frac{\text{NA}^2}{n_{\text{max}}} d_{\text{layer}} \mathcal{H}(r/a), \quad (3) $$

where NA and $a$ are, respectively, the numerical aperture and radius of the MMF, $n_{\text{max}}$ is the value of refractive index $n_0(r)$ at the center of the core, and $\mathcal{H}$ is the Heaviside function.

As detailed in Appendix B, the TM of the $p$th segment of length $l_z$ expressed in the unperturbed fiber mode basis can be written as

$$ \mathbf{T}_p = e^{-i(\mathbf{H}_0 + \mathbf{V})l_z}. \quad (4) $$

Here, $\mathbf{H}_0$ is the propagation operator in the absence of perturbation; it is a diagonal matrix containing the propagation constants $\beta_\mu$ of the modes of the unperturbed fiber, indexed by $\mu$. On the other hand, $\mathbf{V}$ represents the perturbation due to the index fluctuations, $\delta n_p(r, \phi)$, projected onto the mode basis (see Appendix B for further details). The complete TM is obtained by multiplying the TMs of all the segments.

## B. Theoretical predictions for $\langle C(\theta) \rangle$

In the limit of moderate disorder, we are able to find an analytical expression of the mean correlation function $\langle C(\theta) \rangle$, which involves the geometrical parameters of the fiber as well as the disorder strength. In Appendix C, we show that it can be put in the form $\langle C(\theta) \rangle = \tilde{C}(\theta)/\tilde{C}(0)$, with

$$ \tilde{C}(\theta) = 1 + A \sum_{q, \nu, \mu} \Gamma_q^2 \cos(q\theta) B_{\nu\mu}^q. \quad (5) $$

The prefactor $A = N_z(kl_z)^2/4N_{\text{modes}}$ is a coefficient that combines properties of the radial disorder with the number of propagating modes supported by the fiber. In addition, the coefficient $B_{\nu\mu}^q$ characterizes the energy coupling between eigenstates $\psi_\nu$ and $\psi_\mu$ of the unperturbed propagation operator $\mathbf{H}_0$. It is expressed as

$$ B_{\nu\mu}^q = \delta_{m_{\nu\mu}, q} \text{sinc}\left(\frac{\beta_{\nu} - \beta_{\mu}}{2} l_z\right) I_{\nu\mu}, \quad (6) $$

where $m_{\nu\mu} = |m_\mu - m_\nu|$ is the difference between orbital angular momentum of the eigenstates coupled by the azimuthal disorder, and

$$ I_{\nu\mu} = d_{\text{layer}} \int_0^\infty dr |\psi_\nu(r)|^2 |\psi_\mu(r)|^2 \sigma_g^2(r) r^2 \quad (7) $$

is the coupling term induced by disorder along the radial direction.

The expression (5) is a perturbative result, valid when photons scatter on average once over the disordered potential $\mathbf{V}$. As the extent of disorder increases, it becomes necessary to take higher-order perturbations into account. This means taking into account multiple interactions between photons and disorder. For all the results presented in this work, the single scattering contribution Eq. (5) is dominant, but we have also calculated the second-order perturbation contribution to obtain quantitative agreement with experimental results and simulations. The second-order contribution takes the following form:

$$ \tilde{C}^{(2)}(\theta) = \tilde{A} \sum_{q, q', \nu, \kappa, \mu} \Gamma_q^2 \Gamma_{q'}^2 \cos[(q + q')\theta] C_{\nu\kappa\mu}^{qq'}, \quad (8) $$

where $\tilde{A} = N_z(kl_z)^4/16N_{\text{modes}}$. Energy coupling is provided by the term

$$ C_{\nu\kappa\mu}^{qq'} = \frac{N_z - 1}{2} B_{\nu\kappa}^q B_{\kappa\mu}^{q'} + \delta_{m_{\nu\kappa}, q} \delta_{m_{\kappa\nu}, q'} Q_{\mu\kappa\nu} I_{\nu\kappa} I_{\kappa\mu}, \quad (9) $$

where the explicit but lengthy expression of the coefficient $Q_{\mu\kappa\nu}$ in terms of the propagation constants $\beta_\mu, \beta_\kappa$, and $\beta_\nu$ is given in Appendix C.

## C. Validation of the model of disorder and theory

To validate the model of disorder as well our theoretical predictions based on Eqs. (5) and (8), we first find the values of the coefficients $\Gamma_q$ that match best the experimental profile of the mean correlation function $\langle C(\theta) \rangle$. For the graded-index fibers used in our experiments, we find that it is sufficient to use only four nonzero coefficients corresponding to $q \in [1, 2, 3, 4]$. The experimental results shown in Fig. 3 (blue solid lines) are virtually indistinguishable from the analytical curves (black solid lines). To assess the physical relevance of these coefficients, we then perform simulations using the same values of $\Gamma_q$, without adding any fitting parameter. The simulation consists in dividing the fiber into segments of length $l_z$. For each segment, we add to the index profile matching the specifications of the fiber a perturbation of the form given by Eq. (2). We then estimate its TM using a custom fiber solver [48]. Finally, the complete TM of the fiber is obtained by multiplying the TMs of all segments, each segment corresponding to a different realization of the radial disorder (see Appendix B for more details). We calculate
<figure><img src="image_5.png" /></figure>

FIG. 3. Comparison between experiment, simulations, and theory. Mean angular correlation function of the RME, as defined in Eq. (1), for various levels of deformation. The fiber used is a typical graded-index fiber (Prysmian BendBright OM4 [30]), with radius $a = 50$ µm, NA = 0.2, and $N_{modes} = 55$. The correlation length in the model and simulations is set at $l_z = 100$ µm. Experimental data (blue lines) are compared with theoretical predictions based on Eqs. (5) and (8) (black lines, and simulation results for wave propagation inside disordered MMFs (red lines). The parameters $\Gamma_q$ of the model are found by fitting to the experimental results, and simulations are obtained with the same parameters. Error bars represent the standard deviation computed over 100 random input wave fronts for the simulations and experiments, as well as 20 disorder realizations for the simulations.

the RME correlation as a function of angle and average over 20 realizations of the fiber. The results, presented in Fig. 3, show good agreement between simulations and theory, with no adjustment parameters required.

## D. Discussion and interpretation

The different values of the angular momenta $q$ of the deformation have different origins and impacts on the RME. Global radial index variations, corresponding to $q = 0$, alter the shape of the modes but do not break the axisymmetry of the system. In this case, their effect can simply be equated to a change in the length of the fiber up to first order [31] that impacts only the relative phase between the modes. The TM remains diagonal in the mode basis and commutes with the rotation operator $\mathbf{R}(\theta)$ for any angle $\theta$. As a result, the $q = 0$ component does not impact the RME.

In the absence of external perturbation ($\Delta x = 0$), thanks to prior compensation for aberrations, the deviation of the correlation curve from a perfect RME ($C(\theta) = 1$) is due to intrinsic fiber defects caused by the fabrication process, which give rise to nonzero $\Gamma_q$, for $q > 0$. In this regime, we find that the correlation function is dominated by the contributions of even values of $q$. The contribution $q = 2$

is responsible for the valleys found at $\theta = \pi \pm \pi/2$, and the contribution $q = 4$ for the valleys observed at $\theta = \pi/2 \pm \pi/4$ and $\theta = 3\pi/2 \pm \pi/4$. Even contributions have no impact on the value of the correlation at $\theta = \pi$, simply because they correspond to $\pi$-symmetric deformations (see inset of Fig. 4). Consequently, the slight decrease in correlation at $\theta = \pi$ is entirely controlled by the odd deformations.

Although all $\Gamma_q$ terms are of the same order of magnitude, we observe that the effect of odd contributions, which couple modes of different parity to the orbital angular momentum, is much less pronounced than even contributions, which couple modes of the same parity. This is explained by the modal properties of the fiber. Indeed, in ideal graded-index fibers, modes in quasidegenerate groups, i.e., with similar propagation constants $\beta_\mu$, have the same parity of the angular orbital momentum. This property is inherited from the modes of the two-dimensional isotropic harmonic oscillator which represents the idealized parabolic graded-index fiber with no boundary [49–52]. Consequently, for pairs of modes for which $m_{\mu\nu}$ is odd, the difference of propagation constants $\beta_\mu - \beta_\nu$ is non-negligible, leading to weak contributions of $B_{\nu\mu}^q$ appearing in Eq. (5). This effect has the same origin as the observation that disorder preferentially induces coupling between degenerate modes [29].

The previous analysis fully explains the robustness of the correlation observed at $\theta = \pi$ for small deformations. We note that this correlation revival can equivalently be interpreted in terms of commutation between the matrix $\mathbf{T}$ and the rotation matrix $\mathbf{R}(\theta)$. At small deformations, $\mathbf{T}$ is block diagonal, with blocks corresponding to quasidegenerate modes. Since, within each block, the different angular momenta $m_\mu$ share the same parity (corresponding to constant values of $|m_\mu| + 2p_\mu$, where $p_\mu$ is the radial index [38]), the expression of $\mathbf{R}(\theta)$ restricted to each block necessarily satisfies $\mathbf{R}(\theta = \pi) = \pm \mathbf{1}$, where $\mathbf{1}$ is the identity matrix. As a result $\mathbf{R}(\theta = \pi)$ and $\mathbf{T}$ commute, regardless of the coupling complexity within each group of modes.

<figure><img src="image_6.png" /></figure>

FIG. 4. Influence of deformation on the perturbation contributions. Values of the normalized deformation parameters $\tilde{\Gamma}_q = kl_z\sigma_q(r=a)\Gamma_q$. The values of $\Gamma_q$ are found by fitting the theoretical model [Eqs. (5) and (8)] to the experimental data as a function of the deformation. Inset: symmetry corresponding to the perturbation associated with each value of $q$.
When the external mechanical deformation is introduced, we find that the value of $\Gamma_1$, corresponding to a flattening of the fiber, gradually increases (see Fig. 4). This means that modes of different propagation constant become more and more coupled, and that the TM progressively loses its block-diagonal structure in the mode basis. This explains the disappearance of the dominant revival effect at $\theta = \pi$, as well as the loss of the local maxima at $\theta = \pi/2$ and $\theta = 3\pi/2$. As $\Delta x$ is further increased, we also observe that higher harmonics start to play a role (see $\Gamma_3$ in Fig. 4), and the width of the correlation function $\langle C(\theta) \rangle$ starts to decrease. This indicates that the width of the correlation function at large deformation depends on the disorder strength and is intimately connected to the loss of the block-diagonal structure of the TM.

In Appendix D, we present measurements of $\langle C(\theta) \rangle$ obtained for various graded-index fibers, which exhibit advertised properties similar to those of the fiber used in Figs. 3 and 4. Although we obtain qualitatively similar results, we do find some quantitative reproducible differences, expressed in terms of different values for the $\Gamma_q$ weights. This demonstrates that RME is a very good indicator for probing the small variations in disturbances that occur during the MMF manufacturing process.

The possibility to characterize the distributed perturbation along a fiber using intensity measurements from the input and output facets represents important prospects for telecommunication applications. Indeed, once the parameters of the perturbations ($l_z, \Gamma_q$) are found from a fit of the measured RME correlation, one has access to a more accurate model of the fiber. This model can then be exploited through simulations or mode-coupling theory computations to predict the crosstalk and losses given the estimated index profile. This can be used for improving the design and the characterization of new optical fibers with desired properties.

## IV. IMPROVING THE RME FOR IMAGING PROSPECTS

The potential application of the RME extends to the possibility of facilitating blind imaging through an unknown MMF, in a manner similar to demonstrations in scattering media [13,14]. For successful image reconstruction, it is essential to collect information from the output facet, which necessitates maintaining a high RME correlation throughout the entire $2\pi$ range. As shown in the previous section, the memory effect is affected by the presence of disorder stemming from the fabrication process. Notably, this disorder degrades the RME correlation, particularly reducing the range over which the correlation remains close to 1. However, this correlation was obtained by averaging over random input wave fronts. We now ask whether it is possible to find specific input wave fronts for which the correlation is significantly higher than the mean value for one given angle or for a wide range of angles. As with the approaches to

tailoring the angular memory effect in scattering media [53], we can build operators whose eigenstates optimize the memory effect at a given angle. Since losses are low in the fiber, TM is close to unitary and the lower part of Eq. (1) is approximately constant. Then, an interesting operator is the one involved in the upper part of the correlation function:

$$ \mathbf{O}(\theta_t) = \mathbf{T}^{\dagger}\mathbf{R}(\theta_t)^{\dagger}\mathbf{T}\mathbf{R}(\theta_t). \quad (10) $$

This operator can thus be used to improve the RME for a specific value $\theta_t$ of $\theta$, as shown in Appendix E 1. In order to improve the correlation over the entire $2\pi$ range, we can also study the operator built using the sum of operators describing the correlation at different angles:

$$ \mathbf{O}_{\text{sum}} = \sum_{t} \mathbf{T}^{\dagger}\mathbf{R}(\theta_t)^{\dagger}\mathbf{T}\mathbf{R}(\theta_t). \quad (11) $$

To optimize the RME correlation, we construct this operator using the experimentally measured TMs with $\theta_t = t\pi/4$, where $t \in [0, 7]$. We then compute the singular vectors of this operator corresponding to the singular values with the largest modulus. We present in Fig. 5 the resulting correlation $C(\theta)$ for the two first eigenvectors in the case of no deformation and strong deformations ($\Delta x = 60$ µm) and the corresponding output field profiles. Our results demonstrate that it is possible to find input wave fronts for which the output profiles remain highly correlated across the entire $2\pi$ range.

Compared to previous works that aim to enhance the robustness of an input channel against the level of disorder using a Wigner-Smith operator [29], our approach does not depend on differentiating the TM, allowing us to study a

<figure><img src="image_7.png" /><figcaption>FIG. 5. Tailoring the rotational memory effect. The angular correlation function $C(\theta)$ is constructed using experimentally measured input channels with improved RME range, for two values of the deformation ($\Delta x = 0$ µm and $\Delta x = 60$ µm). The results for the first two singular vectors of the operator defined in Eq. (11) are compared with the average results for random input profiles (dashed line). Insets: output spatial transverse profiles of the corresponding singular vectors.</figcaption></figure>
fixed realization of disorder. Furthermore, the projection of
this operator onto a given input wave front directly
represents the targeted quantity, i.e., the correlation of
the RME for this specific wave front.

An interesting application of memory effects is the
ability to retrieve information from the distal side, where
the field for a given input wave front is a *priori* unknown.
For imaging applications, the range of the memory effect
must be wide enough to cover the size of the object to be
imaged, and the output excitation must have a pronounced
peak autocorrelation function [13]. Indeed, even if the
precise shape of the output pattern is unknown, approxi-
mating its autocorrelation by a Dirac function in space
allows for accessing the autocorrelation of the hidden
object. With this information, the image of the object
can then be estimated using a numerical process.

This condition is guaranteed in multiple scattering media
by the presence of strong disorder that randomizes the field
for any given input wave front. However, this is not the case
in MMFs, where the disorder does not affect all modes in
the same way [33]. A trivial solution for maximizing the
RME range would be to use the fundamental mode, which
is less affected by external perturbations [29]. But, due to
its rotational symmetry, the autocorrelation of this mode
with respect to angular rotation is close to one. So, even
though the field profile remains correlated at the fiber's
output when the input profile is rotated, this mode cannot
be used to provide information about the distal end of the
fiber. As shown in Fig. 5, the first singular vector of the
operator Eq. (11) is very close to the fundamental mode for
any $\Delta x$ and is therefore not useful for imaging. However,
the second combines the properties of a large-range RME
and an output pattern with a peaked autocorrelation
function (see Appendix E2 for details). It is therefore a
good candidate for recovering information about the fiber
distal end.

V. CONCLUSION

In this article, we first perform precise measurements of
the RME properties in MMFs. Accurate estimation of mode
propagation properties is essential, as the RME is extremely
sensitive to the effects of aberrations and misalignments. We
then propose a model of disorder and provide a theoretical
calculation of the RME correlation function, which is shown
to be in good agreement with both experimental data and
realistic wave propagation simulations. In particular, our
analysis makes it possible to estimate geometric perturba-
tions in the fiber, whether due to fabrication imperfections or
mechanical deformations. From a fundamental perspective,
this approach can serve as a powerful tool for the study of
MMF defects resulting from the breaking of fiber symmetry.
Moreover, the unknown parameters of the disorder along the
azimuthal and longitudinal directions can be determined
through simple input-output measurements. This provides
promising opportunities for designing and characterizing

fibers for telecommunications applications, aiming to reduce
modal crosstalk that currently restricts their practical utility.

Finally, we tackle the issue of the robustness of the RME
with respect to the rotation angle for a given disorder. This
is an important issue for the prospect of harnessing the
RME for imaging applications. We demonstrate the pos-
sibility of generating channels that exhibit a drastic
improvement in the RME. In particular, we can create
channels that are more robust to deformations compared to
random inputs or standard fiber modes, and that also exhibit
a random profile with high spatial frequencies.

Raw and processed data, sources to regenerate all the
figures, and sample codes for the treatment preprocessing and
postprocessing are available in the dedicated repository [54].

ACKNOWLEDGMENTS

R. G.-C., A. G., J. d. R., and S. M. P. acknowledge the French Agence Nationale pour la Recherche Grant No. ANR-23-CE42-0010-01 MUFFIN and the Labex WIFI Grants No. ANR-10-LABX-24 and No. ANR-10-IDEX-0001-02 PSL*. R. G.-C., A. G., E. R. A., L. B., Y. Q., M. B., P. S., M. B., J. d. R., and S. M. P acknowledge the French Agence Nationale pour la Recherche Grant No. ANR-20-CE24-0016 MUPHTA. O. K. acknowledges funding from the European Research Council under the European Union's Horizon 2020 Research and Innovation Program Grant No. 101002406.

APPENDIX A: TM AND RME MEASUREMENTS

1. Aberration compensation and TM measurement

To decouple the effects of the RME and measurement
inaccuracies, we use the approach we developed to learn
and compensate for aberrations and misalignments in
Ref. [29]. The idea is to first measure the TM on a pixel
basis and then project it onto the mode basis. Without
aberrations, this projection into the mode basis should
conserve energy, since all the energy must be conveyed by
those modes. Using a model-based algorithm, constructed
with the deep-learning framework PyTorch [55], we identify
the aberrations and misalignments of the system that
minimize the loss when projecting onto the mode basis.
First, this process enables us to accurately recover the TM
in the mode basis, which contains all the information about
light propagation in the MMF. Second, it facilitates the
identification of the aberrations that need to be corrected in
order to obtain a desired pattern in the input facet plane of
the fiber. The correction can then be implemented onto the
SLM [56,57].

2. Intensity versus field correlation

In the main text, we studied the RME using the field
correlation function Eq. (1). Another way to characterize
the RME amounts to estimating the correlation between the
<figure><img src="image_8.png" /><figcaption>FIG. 6. Comparison between the field and the intensity correlations as a function of the rotation angle $\theta$. For different values of the deformation $\Delta x$, we show the field correlation as defined in Eq. (1) (blue curve) and the correlation obtained using the TM (red curve). Error bars represent the standard deviation computed over 100 random input wave fronts.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>FIG. 7. Comparison between the correlation $C(\theta)$ based on the measurement of the output fields and the one estimated using the TM. For different values of the deformation $\Delta x$, we show the field correlation as defined in Eq. (1) (blue curve) and the one obtained using the TM (red curve). Error bars represent the standard deviation computed over 100 random input wave fronts.</figcaption></figure>

output intensity patterns associated to the fields $\psi(\mathbf{r})$ and $\psi_{\theta}(\mathbf{r})$ involved in Eq. (1) [27]. It corresponds the output intensity pattern $I(\mathbf{r}) = |\psi(\mathbf{r})|^2$ for a given input, and the intensity $I_{\theta}(\mathbf{r}) = |\psi_{\theta}(\mathbf{r})|^2$ obtained when rotating the input and output profiles by an angle $\theta$. The intensity correlation reads

$$C_I(\theta) = \frac{\int d\mathbf{r} I(\mathbf{r}) I_{\theta}(\mathbf{r})}{\sqrt{\int d\mathbf{r} I(\mathbf{r})^2 \int d\mathbf{r} I_{\theta}(\mathbf{r})^2}}. \quad (A1)$$

Such correlation function was originally used in Ref. [27]. To estimate it, we use the same procedure and apparatus as the ones discussed in Sec. II B, where intensities can be evaluated from the field measurements. We then compare the mean field correlation Eq. (1) with the square root of the mean intensity correlation Eq. (A1). As illustrated in Fig. 6, the mean values of the two correlation functions are in excellent agreement. This demonstrates that the observables $\langle C(\theta) \rangle$ and $\langle C_I(\theta) \rangle$ are equivalent.

### 3. Estimation of the correlation using the measured TM

In the present study, we use the compensation of aberrations, facilitated by the transmission matrix measurement, but we do not directly employ the knowledge of the TM itself. However, the TM provides access to the output field $|\psi_{\text{out}}\rangle$ for any given input field $|\psi_{\text{in}}\rangle$. We can thus use the TM to estimate the output of a rotated wave front and compute the correlation function defined in Eq. (1). For each deformation, we compute the mean correlation for 100 random input wave fronts. We show in Fig. 7 a good agreement between the estimation based on

the TM and the one based on the explicit measurement of the output field. This demonstrates that the measurement of the TM can drastically reduce the time needed for characterizing the RME, as it does not require any additional measurement. In comparison, the explicit measurements presented in the main text necessitate, after the initial calibration, to average over 100 random input wave fronts for 50 different angles.

## APPENDIX B: EFFECTIVE HAMILTONIAN AND TRANSMISSION MATRIX

In the situation where the coupling between different polarization channels can be neglected, and in the weakly guiding approximation [47] (i.e., for variations of the index of refraction small compared to the average index value), each polarization of the transverse part of the field at frequency $\omega$ satisfies the scalar wave equation,

$$[\nabla_{\perp}^2 + \partial_z^2 + k^2 n(\mathbf{r}, z)^2]\psi(\mathbf{r}, z) = 0, \quad (B1)$$

where $k = \omega/c$ and $\mathbf{r} = (r, \phi)$ labels the position in the transverse plane. The refractive index is further decomposed into an unperturbed axisymmetric component $n_0$ and a perturbation $\delta n \ll n_0$,

$$n(\mathbf{r}, z) = n_0(r) + \delta n(\mathbf{r}, z). \quad (B2)$$

To identify the effective Hamiltonian that controls the dynamics in the presence of disorder, it is convenient to rewrite the wave equation in the operator form,

$$\partial_z^2 |\psi(z)\rangle = -\hat{H}(z)^2 |\psi(z)\rangle, \quad (B3)$$
where

$$ 
\begin{aligned}
\hat{H}(z) &= [\hat{\nabla}_{\perp}^2 + k^2 \hat{n}(\mathbf{r}, z)^2]^{1/2} \\
&\approx [\hat{H}_0^2 + 2k^2 \hat{n}_0(r) \delta\hat{n}(\mathbf{r}, z)]^{1/2} \\
&\approx \hat{H}_0 + k^2 \hat{H}_0^{-1} \hat{n}_0(r) \delta\hat{n}(\mathbf{r}, z).
\end{aligned} 
\quad (B4) $$

The Hamiltonian of the unperturbed problem reads $\hat{H}_0 = [\hat{\nabla}_{\perp}^2 + k^2 \hat{n}_0(r)^2]^{1/2}$. Since in a realistic MMF the relative variations of $n_0(r)$ in the radial direction are small, the eigenvalues $\beta_{\mu}$ of $\hat{H}_0$ are close to $kn_0$, where $n_0$ is the typical refractive index of the core. Therefore, a good approximation of $\hat{H}$ is

$$ \hat{H}(z) \approx \hat{H}_0 + k \delta \hat{n}(\mathbf{r}, z). \quad (B5) $$

In the present work, back reflections can be neglected, and Eq. (B3) is equivalent to

$$ \partial_z |\psi(z)\rangle = -i \hat{H}(z) |\psi(z)\rangle. \quad (B6) $$

This shows that the field transmitted through the fiber of length $L$ can be expressed in terms of a unitary transmission matrix $\mathbf{T}$ as $|\psi(L)\rangle = \mathbf{T}|\psi(0)\rangle$. The matrix $\mathbf{T}$ reads

$$ \mathbf{T} = \mathcal{T} e^{-i \int_0^L d\zeta' \hat{H}(\zeta')}, \quad (B7) $$

where $\mathcal{T}$ is the time-ordering operator ($z$ plays the role of time here). In the following, we model the disorder along the propagation direction $z$ as a succession of $N_z = L/l_z$ independent segments of length $l_z$, where the refractive index depends only on the transverse coordinate $\mathbf{r}$. In that case, the transmission matrix takes the form

$$ \mathbf{T} = \prod_{p=1}^{N_z} \mathbf{T}^{(p)}, \quad (B8) $$

with

$$ \mathbf{T}^{(p)} = e^{-i[\hat{H}_0 + k\delta\hat{n}_p(\mathbf{r})]l_z}. \quad (B9) $$

The index fluctuations of each sector $p$ are expressed as the product of a random function along the radial direction and a random function decomposed on the azimuthal harmonics:

$$ \delta n_p(r, \phi) = g_p(r) \sum_q \Gamma_q \cos(q\phi + \varphi_q). \quad (B10) $$

Here, $g_p(r)$ is a Gaussian random variable with zero mean and variance $\langle g_p(r)g_p(r') \rangle = \sigma_g(r)^2 d_{\text{layer}}\delta(r-r')$, where $d_{\text{layer}}$ is the thickness of each layer obtained in the chemical vapor deposition process. In addition, the phases $\varphi_q$ are

random independent variables with uniform distribution, added to mitigate the effect of the orientation of the perturbation.

In this work, we focus on the properties of graded-index fibers, where the refractive index $n_0(r)$ takes the form

$$ n_0(r)^2 = n_{\max}^2 \left( 1 - 2\Delta \frac{r^2}{a^2} \right) \quad (B11) $$

in the core of the fiber of radius $a$. Here $\Delta = (n_{\max} - n_{\text{cl}})/n_{\max}$, where $n_{\text{cl}}$ is the refractive index in the cladding, i.e., for $r > a$. In the weakly guiding approximation (NA $\ll$ 1), $\Delta \approx \text{NA}^2/2n_{\max}^2$ and the refractive index profile in the core is well approximated by a parabolic function, $n_0(r) \approx n_{\max}(1 - \Delta r^2/a^2)$. This yields the explicit expression (3) for the amplitude of the radial disorder $\sigma_g(r)$.

The expressions (B8), (B9), and (B10) are used both in the theoretical treatment developed in Appendix A 2 and in the numerical simulations. For simulation purposes, the modes profiles $\psi_{\mu}$ and propagation constants $\beta_{\mu}$ of the unperturbed fiber (which are the eigenstates and eigenvalues of $\hat{H}_0$) are computed using the pyMMF package [29,48]. The Hamiltonian Eq. (B5) and transmission matrix Eq. (B9) of each sector $p$ is then computed in the basis $\{\psi_{\mu}\}$. Finally, the total TM is found by multiplying the TMs of all the segments, as in Eq. (B8). Details of the simulations, performed in Python, are available in the dedicated repository [54].

## APPENDIX C: ANALYTICAL PREDICTIONS FOR THE RME

In this appendix, we evaluate the mean correlator $\langle C(\theta) \rangle = \tilde{C}(\theta)/\tilde{C}(0)$, where

$$ \tilde{C}(\theta) = \overline{\langle \psi | \psi_{\theta} \rangle} = \overline{\langle \psi_{\text{in}} | \mathbf{T}^{\dagger} \mathbf{T}_{\theta} | \psi_{\text{in}} \rangle}, \quad (C1) $$

and $\overline{\cdot \cdot \cdot} = \langle \cdots \rangle$ stands for the average over different configurations of the disorder. We first decompose the input field in the mode basis $\{\psi_{\mu}\}$ of the unperturbed Hamiltonian $\hat{H}_0$,

$$ |\psi_{\text{in}}\rangle = \sum_{\mu=1}^{N} c_{\mu} |\psi_{\mu}\rangle, \quad (C2) $$

where $\sum_{\mu=1}^{N} |c_{\mu}|^2 = 1$. In the following, we write the unperturbed eigenmodes in the form

$$ \psi_{\mu}(r, \phi) = \frac{1}{\sqrt{2\pi}} \varphi_{\mu}(r) e^{im_{\mu}\phi}, \quad (C3) $$

so that the normalization condition $\langle \psi_{\mu} | \psi_{\mu} \rangle = 1$ reads
$$ \int_0^\infty dr r |\varphi_\mu(r)|^2 = 1. \qquad (C4) $$

In addition, we consider random input wave fronts, uniformly distributed over the $N_{modes}$ modes of the MMF. Using $\langle c_\mu c_{\mu'} \rangle = \delta_{\mu,\mu'}/N_{modes}$, we express the correlator Eq. (C1) as

$$ \tilde{C}(\theta) = \frac{1}{N_{\text{modes}}} \sum_{\nu, \mu} e^{-i(m_{\nu} - m_{\mu})\theta} \langle |T_{\nu\mu}|^2 \rangle. \quad (C5) $$

We then use the decomposition Eq. (B8), where TMs $\mathbf{T}^{(p)}$ are independent of each other, and satisfy $\langle \mathbf{T}^{(p)} \rangle = 0$. This gives

$$ \langle |T_{\nu\mu}|^2 \rangle = \left( \prod_{p=1}^{N_z} \langle \mathbf{T}^{(p)} \otimes \mathbf{T}^{(p)\dagger} \rangle \right)_{\nu\mu}. \quad (C6) $$

In the case of weak disorder ($N_z(k_lz)^2 \langle \delta n^2 \rangle \lesssim 1$), we can evaluate the previous correlator using a perturbative expansion of each matrix $\mathbf{T}^{(p)}$. To obtain an explicit form of the latter, it is more convenient to work with the interaction representation $\mathbf{T}_I(z) = e^{i\hat{H}_{0z}}\mathbf{T}(z)$ than directly using the expansion of Eq. (B9). As the matrix $\mathbf{T}_I(z)$ obeys the equation $\partial_z \mathbf{T}_I(z) = -i\hat{V}_I(z)\mathbf{T}_I(z)$, where $\hat{V}_I(z) = e^{i\hat{H}_{0z}}\hat{V}e^{-i\hat{H}_{0z}}$ and $\hat{V}(z) = k\delta\hat{n}(z)$, it can be expanded, up to the second order in $\hat{V}_I$, in the form

$$ \begin{aligned} \mathbf{T}_I(z) &= 1 - i \int_0^z dz' \hat{V}_I(z') \mathbf{T}_I(z') \\ &\approx 1 - i \int_0^z dz' \hat{V}_I(z') - \int_0^z dz' \int_0^z dz'' \hat{V}_I(z') \hat{V}_I(z''). \end{aligned} \quad (C7) $$

Physically, this expansion corresponds to a situation where photons interact at most twice with the disordered potential located in a section of the fiber of length $z$. Inside each sector of length $l_z$, the potential $\hat{V}(z)$ is invariant along $z$, so that integrals in Eq. (C7) can be evaluated explicitly. This allows us to find the expression of $\mathbf{T}^{(p)} = e^{-i\hat{H}_{0z}}\mathbf{T}_I(l_z)$, up to second order in $\hat{V} = k\delta\hat{n}_p$,

$$ T_{\nu\mu}^{(p)} \approx e^{-i\hat{\beta}_{\nu\mu}l_z} (\delta_{\nu\mu} + T_{\nu\mu}^{(p,1)} + T_{\nu\mu}^{(p,2)}), \quad (C8) $$

where

$$ T_{\nu\mu}^{(p,1)} = -il_z e^{i\hat{\beta}_{\nu\mu}l_z/2} \mathrm{sinc}(\beta_{\nu\mu}l_z/2) V_{\nu\mu}, \quad (C9) $$

$$ T_{\nu\mu}^{(p,2)} = -il_z \sum_{\kappa} \frac{e^{i\hat{\beta}_{\nu\kappa}}}{\beta_{\nu\mu}} \left[ \begin{aligned} & e^{i\hat{\beta}_{\kappa\mu}l_z/2} \mathrm{sinc}(\beta_{\kappa\mu}l_z/2) \\ & + e^{-i\hat{\beta}_{\kappa\mu}l_z/2} \mathrm{sinc}(\beta_{\kappa\mu}l_z/2) \end{aligned} \right] V_{\nu\kappa} V_{\kappa\mu}, \quad (C10) $$

with $\beta_{\nu\mu} = \beta_\nu - \beta_\mu$. Inserting the expansion Eq. (C8) into Eq. (C6) and keeping terms up to second order in $V$, we obtain

$$ \begin{aligned} \langle |T_{\nu\mu}|^2 \rangle \approx & \delta_{\nu\mu} + N_z \langle |T_{\nu\mu}^{(p,1)}|^2 \rangle + N_z \langle |T_{\nu\mu}^{(p,2)}|^2 \rangle \\ & + \frac{N_z(N_z - 1)}{2} \sum_{\kappa} \langle |T_{\nu\kappa}^{(p,1)}|^2 \rangle \langle |T_{\kappa\mu}^{(p,1)}|^2 \rangle. \end{aligned} \quad (C11) $$

First-order contributions are of the form

$$ \langle |T_{\nu\mu}^{(p,1)}|^2 \rangle = I_z^2 \mathrm{sinc}(\beta_{\nu\mu} l_z/2)^2 \langle |V_{\nu\mu}|^2 \rangle, \quad (C12) $$

where $V_{\nu\mu} = k \langle \psi_\nu | \delta \hat{n}_p | \psi_\mu \rangle$. For the model of disorder given by Eq. (B10), we find

$$ \langle |V_{\nu\mu}|^2 \rangle = \frac{k^2}{4} I_{\nu\mu} \sum_q \Gamma_q^2 \delta_{q, |m_\nu - m_\mu|}, \quad (C13) $$

where

$$ I_{\nu\mu} = d_{\text{layer}} \int dr |\psi_{\nu}(r)|^2 |\psi_{\mu}(r)|^2 \sigma_g(r)^2 r^2. \quad (C14) $$

Second-order contributions $\langle |T_{\nu\mu}^{(p,2)}|^2 \rangle$ involve averages of the form $C_{\nu\kappa\mu}^{\nu\kappa\mu} = \langle V_{\nu\kappa} V_{\kappa\mu} V_{\nu\kappa'}^* V_{\kappa'\mu}^* \rangle$, which we can contract as

$$ \begin{aligned} C_{\nu\kappa\mu}^{\nu\kappa\mu} &= \langle V_{\nu\kappa} V_{\nu\kappa'}^* \rangle \langle V_{\kappa\mu} V_{\kappa'\mu}^* \rangle + \langle V_{\nu\kappa} V_{\kappa\mu}^* \rangle \langle V_{\kappa\mu} V_{\nu\kappa'}^* \rangle \\ &\approx \langle |V_{\nu\kappa}|^2 \rangle \langle |V_{\kappa\mu}|^2 \rangle \delta_{\kappa\kappa'} + \langle |V_{\nu\mu}|^2 \rangle^2 \delta_{\kappa\kappa'} \delta_{\nu\mu} \delta_{\nu\kappa}. \end{aligned} \quad (C15) $$

Combining the expression (C10) with the previous result, we obtain

$$ \langle |T_{\nu\mu}^{(p,2)}|^2 \rangle \approx l_z^4 \sum_{\kappa} Q_{\nu\kappa\mu} \langle |V_{\nu\kappa}|^2 \rangle \langle |V_{\kappa\mu}|^2 \rangle, \quad (C16) $$

where $Q_{\nu\kappa\mu}$ is a coupling weight between different energy subspaces:

$$ \begin{aligned} Q_{\nu\kappa\mu} = & \frac{1}{\beta_{\nu\mu}^2 l_z^2} [\mathrm{sinc}(\beta_{\nu\kappa} l_z/2)^2 + \mathrm{sinc}(\beta_{\kappa\mu} l_z/2)^2 \\ & -2\mathrm{sinc}(\beta_{\nu\kappa} l_z/2)\mathrm{sinc}(\beta_{\kappa\mu} l_z/2) \cos(\beta_{\nu\mu} l_z/2)] \\ & + \frac{1}{4} \delta_{\nu\kappa} \delta_{\kappa\mu}. \end{aligned} \quad (C17) $$

Finally, we insert the result Eq. (C11) into the expression (C5) of the correlator, to get an expansion of the form

$$ \tilde{C}(\theta) = 1 + \tilde{C}^{(1)}(\theta) + \tilde{C}^{(2)}(\theta). \quad (C18) $$

The first order in $V^2$ reads

$$ \tilde{C}^{(1)}(\theta) = A_1 \sum_{\substack{q, v, \mu \\ m_v - m_\mu = \pm q}} \Gamma_q^2 \cos(q\theta) I_{\nu\mu} \mathrm{sinc}\left(\frac{\beta_{\nu\mu} l_z}{2}\right), \quad (C19) $$
where $A_1 = N_z(kl_z)^2/4N_{\text{modes}}$. On the other hand, the second order in $V^2$ reads

$$ \tilde{C}^{(2)}(\theta) = A_2 \sum_{\substack{q, q', \nu, \kappa, \mu \\ m_{\nu} - m_{\kappa} = q \\ m_{\kappa} - m_{\mu} = \pm q'}} \Gamma_q^2 \Gamma_{q'}^2 \cos[(q \pm q')\theta] I_{\nu\kappa} I_{\kappa\mu} \tilde{Q}_{\nu\kappa\mu}, \quad (\text{C20}) $$

where $A_2 = N_z(kl_z)^4/8N_{\text{modes}}$, and

$$ \tilde{Q}_{\nu\kappa\mu} = Q_{\nu\kappa\mu} + \frac{N_z - 1}{2} \text{sinc}\left(\frac{\beta_{\nu\kappa} l_z}{2}\right)^2 \text{sinc}\left(\frac{\beta_{\kappa\mu} l_z}{2}\right)^2. \quad (\text{C21}) $$

Equations (C18)–(C20) are equivalent to Eqs. (5) and (8) of the main text. They are used to generate the theoretical predictions shown in Figs. 3 and 8.

<figure><img src="image_10.png" /><figcaption>FIG. 8. RME correlation results for a batch of GIF50C. (a) Angular correlation function of the RME, as defined in Eq. (1), for various levels of deformation $\Delta x$. Experimental data (blue) are compared to the theoretical prediction based on Eqs. (5) and (8) (black) and to simulation results obtained with the same parameters as those used in the theoretical model (red). Error bars represent the standard deviation computed over 100 random input wave fronts for the simulations and experiments, as well as 20 disorder realizations for the simulations. (b) Values of the normalized deformation parameters $\tilde{\Gamma}_q = kl_z\sigma_g(r=a)\Gamma_q$. The values of $\Gamma_q$ are found by fitting the theoretical model [Eq. (5)] to the experimental data as a function of the deformation. Inset: symmetry corresponding to the perturbation associated with each value of $q$.</figcaption></figure>

## APPENDIX D: RME CHARACTERIZATION OF DIFFERENT GRADED-INDEX FIBERS

In this appendix, we report the measurements for different fiber segments of the same length ($L = 24.5$ cm), and with advertised properties similar to those of the fiber used in the main text. Specifically, we use samples from a Thorlabs 50-µm core OM2 graded-index fiber (GIF50C, NA = 0.2).

Results for different fiber segments of the same spool are reproducible. We present typical results for one sample in Fig. 8. We observe different contributions of the $\Gamma_q$ terms as the ones reported in Fig. 3, where a Prysmian BendBright OM4 fiber was used [30]. In particular, $\Gamma_4$ is much smaller, leading to the absence of observed local maxima of the correlation at $\pi/2$ and $3\pi/2$.

## APPENDIX E: PROPERTIES OF THE RME CHANNELS

### 1. RME Operator for one angle value

We consider here operator $O(\theta_t)$ defined in Eq. (10), which represents the upper part of the correlation function Eq. (1). Computing the singular values of this operator enables the identification of input wave fronts that maximize the angular correlation for a specific value $\theta_t$ of $\theta$. We present in Fig. 9 the resulting correlation $C(\theta)$ of the first two singular vectors for $\theta_t = \pi/2$, in the cases of no deformation and strong deformations ($\Delta x = 60$ µm), along with the corresponding output field profiles. As with the results presented for the sum operator in Fig. 5, the first

<figure><img src="image_11.png" /><figcaption>FIG. 9. Tailoring the rotational memory effect. The angular correlation function $C(\theta)$ is constructed using experimentally measured input channels with improved RME range, for two values of the deformation ($\Delta x = 0$ µm and $\Delta x = 60$ µm). The results for the first two singular vectors of the operator defined in Eq. (11) are compared with the average results for random input profiles (dashed line). Insets: output spatial transverse profiles of the corresponding singular vectors.</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>FIG. 10. Autocorrelation of the RME channels. Experimental rotational autocorrelation of the output field when injecting the wave front associated with the first (solid dark blue line) and second (dashed light blue line) highest singular modes, without and with external deformations. For comparison, we also show the average autocorrelation of the output field for random input wave front (dashed black line), that represents the optimal autocorrelation.</figcaption></figure>

singular vector of the operator Eq. (10) closely resembles the fundamental mode for any $\Delta x$. The second singular vector exhibits higher spatial frequencies and achieves a maximum in the correlation at the target angular value $\theta_t$. However, although the correlation curve is consistently higher than the average correlation, it displays significant fluctuations over the $2\pi$ range.

## 2. Autocorrelation of the RME channels

As stated in the main text, for efficient information retrieval from the hidden side of a complex medium, the output pattern used should possess both a substantial memory effect range and a narrow autocorrelation function. This implies that the correlator $C(\theta)$ defined in Eq. (1) displays a broad width, and that the angular autocorrelation,

$$C_0(\theta) = \frac{|\langle\psi|\mathbf{R}(\theta)|\psi\rangle|}{\langle\psi|\psi\rangle}, \quad (E1)$$

of the transmitted field $|\psi\rangle = \mathbf{T}|\psi_{in}\rangle$ is sharply peaked.

We present in Fig. 10 the autocorrelation of the first two RME eigenchannels of the operator $\mathbf{O}_{\text{sum}}$ defined in Eq. (11), under different conditions—absence of external deformations and presence of strong deformations

($\Delta x = 60 \mu\text{m}$). As anticipated in the main text, the first singular mode, being less susceptible to perturbations, exhibits a rotational symmetry with a strong correlation across the entire angular range, even under significant deformations. This indicates that this mode, despite its wide RME range, does not meet the criteria defined in the main text to be a viable candidate for image or information retrieval. However, the second singular mode does exhibit a pronounced RME peak around $\theta = 0$. The combination of this property with the wide angular range of the RME makes this mode an excellent candidate for effective use in discerning the angular properties of an object or signal obscured at the fiber distal end.

[1] A. Boese, C. Wex, R. Croner, U. B. Liehr, J. J. Wendler, J. Weigt, T. Walles, U. Vorwerk, C. H. Lohmann, M. Friebe, and A. Illanes, *Endoscopic imaging technology today*, Diagnostics **12**, 1262 (2022).

[2] S. Turtaev, I. T. Leite, and T. Čižmár, *Multimode fibres for micro-endoscopy*, Optofluid. Microfluid. Nanofluid. **2**, 31 (2015).

[3] M. Plöschner, T. Tyc, and T. Čižmár, *Seeing through chaos in multimode fibres*, Nat. Photonics **9**, 529 (2015).

[4] T. Čižmár and K. Dholakia, *Exploiting multimode waveguides for pure fibre-based imaging*, Nat. Commun. **3**, 1 (2012).

[5] Y. Choi, C. Yoon, M. Kim, T. D. Yang, C. Fang-Yen, R. R. Dasari, K. J. Lee, and W. Choi, *Scanner-free and wide-field endoscopic imaging by using a single multimode optical fiber*, Phys. Rev. Lett. **109**, 203901 (2012).

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

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

[8] P. Jákl, M. Šiler, J. Ježek, A. Cifuentes, J. Trägårdh, P. Zemánek, and T. Čižmár, *Endoscopic imaging using a multimode optical fibre calibrated with multiple internal references*, Photonics **9**, 37 (2022).

[9] J. Yammine, A. Tandjè, M. Dossou, L. Bigot, and E. R. Andresen, *Time-dependence of the transmission matrix of a specialty few-mode fiber*, APL Photonics **4**, 022904 (2019).

[10] R. Maruyama, N. Kuwaki, S. Matsuo, and M. Ohashi, *Relationship between mode coupling and fiber characteristics in few-mode fibers analyzed using impulse response measurements technique*, J. Lightwave Technol. **35**, 650 (2017).

[11] H. Yu, J. Park, K. Lee, J. Yoon, K. Kim, S. Lee, and Y. Park, *Recent advances in wavefront shaping techniques for biomedical applications*, Curr. Appl. Phys. **15**, 632 (2015).

[12] S. Yoon, S. Y. Cheon, S. Park, D. Lee, Y. Lee, S. Han, M. Kim, and H. Koo, *Recent advances in optical imaging through deep tissue: Imaging probes and techniques*, Biomater. Res. **26**, 57 (2022).
[13] J. Bertolotti, E. G. van Putten, C. Blum, A. Lagendijk, W. L. Vos, and A. P. Mosk, *Non-invasive imaging through opaque scattering layers*, Nature (London) **491**, 232 (2012).

[14] O. Katz, P. Heidmann, M. Fink, and S. Gigan, *Non-invasive single-shot imaging through scattering layers and around corners via speckle correlations*, Nat. Photonics **8**, 784 (2014).

[15] T. Yeminy and O. Katz, *Guidestar-free image-guided wavefront shaping*, Sci. Adv. **7**, eabf5364 (2021).

[16] Y. Bromberg, B. Redding, S. M. Popoff, and H. Cao, *Control of coherent backscattering by breaking optical reciprocity*, Phys. Rev. A **93**, 023826 (2016).

[17] L. Devaud, M. Guillon, I. Gusachenko, and S. Gigan, *Chromato-axial memory effect in step-index multimode fibers*, APL Photonics **6**, 126105 (2021).

[18] J. Carpenter, B. J. Eggleton, and J. Schröder, 110 × 110 optical mode transfer matrix inversion, Opt. Express **22**, 96 (2014).

[19] J. Carpenter, B. J. Eggleton, and J. Schröder, *Observation of Eisenbud–Wigner–Smith states as principal modes in multi-mode fibre*, Nat. Photonics **9**, 751 (2015).

[20] W. Xiong, P. Ambichl, Y. Bromberg, B. Redding, S. Rotter, and H. Cao, *Spatiotemporal control of light transmission through a multimode fiber with strong mode coupling*, Phys. Rev. Lett. **117**, 053901 (2016).

[21] P. Chiarawongse, H. Li, W. Xiong, C. W. Hsu, H. Cao, and T. Kottos, *Statistical description of transport in multimode fibers with mode-dependent loss*, New J. Phys. **20**, 113028 (2018).

[22] Y. Li, D. Cohen, and T. Kottos, *Coherent wave propagation in multimode systems with correlated noise*, Phys. Rev. Lett. **122**, 153903 (2019).

[23] Y. Li, D. Cohen, and T. Kottos, *Enforcing Levy relaxation for multi-mode fibers with correlated disorder*, New J. Phys. **24**, 053012 (2022).

[24] A. M. Caravaca-Aguirre, A. Carron, S. Mezil, I. Wang, and E. Bossy, *Optical memory effect in square multimode fibers*, Opt. Lett. **46**, 4924 (2021).

[25] S. Mezil, I. Wang, and E. Bossy, *Imaging through a square multimode fiber by scanning focused spots with the memory effect*, Opt. Lett. **48**, 4701 (2023).

[26] D. Bouchet, A. M. Caravaca-Aguirre, G. Godefroy, P. Moreau, I. Wang, and E. Bossy, *Speckle-correlation imaging through a kaleidoscopic multimode fiber*, Proc. Natl. Acad. Sci. U.S.A. **120**, e2221407120 (2023).

[27] L. V. Amitonova, A. P. Mosk, and P. W. H. Pinkse, *Rotational memory effect of a multimode fiber*, Opt. Express **23**, 20569 (2015).

[28] S. Li, S. A. R. Horsley, T. Tyc, T. Čižmár, and D. B. Phillips, *Memory effect assisted imaging through multimode optical fibres*, Nat. Commun. **12**, 3751 (2021).

[29] M. W. Matthès, Y. Bromberg, J. de Rosny, and S. M. Popoff, *Learning and avoiding disorder in multimode fibers*, Phys. Rev. X **11**, 021060 (2021).

[30] Datasheet for BendBright™ OM4, https://www.prysmiangroup.com/en/bendbright-tm-om4.

[31] D. E. Boonzajer Flaes, J. Stopka, S. Turtaev, J. F. de Boer, T. Tyc, and T. Čižmár, *Robustness of light-transport processes to bending deformations in graded-index multimode waveguides*, Phys. Rev. Lett. **120**, 233901 (2018).

[32] S. M. Popoff, G. Lerosey, R. Carminati, M. Fink, A. C. Boccara, and S. Gigan, *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).

[33] T. Čižmár and K. Dholakia, *Shaping the light transmission through a multimode optical fibre: Complex transformation analysis and applications in biophotonics*, Opt. Express **19**, 18871 (2011).

[34] I. M. Vellekoop and A. P. Mosk, *Focusing coherent light through opaque strongly scattering media*, Opt. Lett. **32**, 2309 (2007).

[35] T. Čižmár, M. Mazilu, and K. Dholakia, *In situ wavefront correction and its application to micromanipulation*, Nat. Photonics **4**, 388 (2010).

[36] B. J. Puttnam, G. Rademacher, and R. S. Luís, *Space-division multiplexing for optical fiber communications*, Optica **8**, 1186 (2021).

[37] D. Marcuse and H. M. Presby, *Mode coupling in an optical fiber with core distortions*, Bell Syst. Tech. J. **54**, 3 (1975).

[38] R. Olshansky, *Mode coupling effects in graded-index optical fibers*, Appl. Opt. **14**, 935 (1975).

[39] L. Gruner-Nielsen, Y. Sun, J. W. Nicholson, D. Jakobsen, K. G. Jespersen, R. Lingle, and B. Palsdottir, *Few mode transmission fiber with low DGD, low mode coupling, and low loss*, J. Lightwave Technol. **30**, 3693 (2012).

[40] P. Sillard, M. Bigot-Astruc, D. Boivin, H. Maerten, and L. Provost, *Few-mode fiber for uncoupled mode-division multiplexing transmissions*, in *Proceedings of the 37th European Conference and Exposition on Optical Communications, OSA Technical Digest (CD)* (Optica Publishing Group, 2011), paper Tu.5.LeCervin.7, 10.1364/ECOC.2011.Tu.5.LeCervin.7.

[41] K.-P. Ho, J. M. Kahn, I. Kaminow, T. Li, and A. Willner, *Mode coupling and its impact on spatially multiplexed systems*, Opt. Fiber Telecommun. VI **17**, 1386 (2013).

[42] D. Marcuse, *Coupled mode theory of round optical fibers*, Bell Syst. Tech. J. **52**, 817 (1973).

[43] P. Mazumder, S. L. Logunov, and S. Raghavan, *Analysis of excess scattering in optical fibers*, J. Appl. Phys. **96**, 4042 (2004).

[44] M. Bsaibes, Y. Quiquempois, E.-R. Andresen, M. Bigot, J.-B. Trinel, P. Sillard, and L. Bigot, *Mode coupling in few-mode optical fibers: The impact of light scattering (unpublished)*.

[45] H. Lydtin, *PCVD: A technique suitable for large-scale fabrication of optical fibers*, J. Lightwave Technol. **4**, 1034 (1986).

[46] P. Geittner and H. Lydtin, *Manufacturing optical fibres by the PCVD process*, Philips Tech. Rev. **44**, 241 (1989), https://pearl-hifi.com/06_Lit_Archive/02_PEARL_Arch/Vol_16/Sec_53/Philips_Tech_Review/PTechReview-44-1988_89-241.pdf.

[47] D. Marcuse, *Theory of Dielectric Optical Waveguides* (Elsevier, New York, 2013).

[48] S. M. Popoff and P. Gostev, *wavefrontshaping/pyMMF: 0.5, 10.5281/zenodo.5524192*.

[49] S. van Enk, *Geometric phase, transformations of Gaussian light beams and angular momentum transfer*, Opt. Commun. **102**, 59 (1993).
[50] R. Gutiérrez-Cuevas, M. R. Dennis, and M. A. Alonso, *Generalized Gaussian beams in terms of Jones vectors*, J. Opt. **21**, 084001 (2019).

[51] R. Gutiérrez-Cuevas, S. A. Wadood, A. N. Vamivakas, and M. A. Alonso, *Modal Majorana sphere and hidden symmetries of structured-Gaussian beams*, Phys. Rev. Lett. **125**, 123903 (2020).

[52] R. Gutiérrez-Cuevas, D. H. J. O'Dell, M. R. Dennis, and M. A. Alonso, *Exactly solvable model behind Bose-Hubbard dimers, Ince-Gauss beams, and aberrated optical cavities*, Phys. Rev. A **107**, L031502 (2023).

[53] H. Yilmaz, M. Kühmayer, C. W. Hsu, S. Rotter, and H. Cao, *Customizing the angular memory effect for scattering media*, Phys. Rev. X **11**, 031010 (2021).

[54] R. Gutiérrez-Cuevas and S. M. Popoff, *Code examples and supporting data for the article “Tailoring the rotational memory effect in multimode fibers,” https://github.com/wavefrontshaping/article_RME_2023, accessed Sept. 1, 2023.*

[55] A. Paszke, S. Gross, F. Massa, A. Lerer, J. Bradbury, G. Chanan, T. Killeen, Z. Lin, N. Gimelshein, L. Antiga et al., *PyTorch: An imperative style, high-performance deep learning library*, in *Advances in Neural Information Processing Systems* 32, edited by H. Wallach, H. Larochelle, A. Beygelzimer, F. d'Alché Buc, E. Fox, and R. Garnett (Curran Associates, Inc., Red Hook, NY, 2019), pp. 8024–8035.

[56] S. M. Popoff, Y. Bromberg, M. W. Matthès, and R. Gutiérrez-Cuevas, *A practical guide to digital micro-mirror devices (DMDs) for wavefront shaping*, J. Phys. Photonics **6**, 043001 (2024).

[57] R. Gutiérrez-Cuevas and S. M. Popoff, *Binary holograms for shaping light with digital micromirror devices*, arXiv: 2311.16685.