# Characterising optical fibre transmission matrices using metasurface
reflector stacks for lensless imaging without distal access

George S.D. Gordon  Department of Electrical and Electronic Engineering,
University of Nottingham<EMAIL_ADDRESS>Milana Gataric
Department of Pure Mathematics and Mathematical Statistics, University of
Cambridge  Alberto Gil C. P. Ramos  Nokia Bell Labs, Cambridge, UK  Ralf
Mouthaan  Electrical Engineering Division, Dept. of Engineering, University of
Cambridge  Department of Physics, University of Cambridge  Calum Williams
Electrical Engineering Division, Dept. of Engineering, University of Cambridge
Jonghee Yoon  Department of Physics, University of Cambridge  Timothy D.
Wilkinson  Electrical Engineering Division, Dept. of Engineering, University
of Cambridge  Sarah E. Bohndiek  Department of Physics, University of
Cambridge  CRUK Cambridge Institute, University of Cambridge

###### Abstract

The ability to form images through hair-thin optical fibres promises to open
up new applications from biomedical imaging to industrial inspection.
Unfortunately, their deployment has been limited because small changes in
mechanical deformation (e.g. bending) and temperature can completely scramble
optical information, which distorts the resulting images. Since such changes
are dynamic, correcting them requires measurement of the fibre transmission
matrix _in situ_ immediately before imaging. Transmission matrix calibration
typically requires access to both the proximal and distal facets of the fibre
simultaneously, which is not feasible during most realistic usage scenarios
without compromising the thin form factor with bulky distal optics. Here, we
introduce a new approach to determine the transmission matrix of multi-mode or
multi-core optical fibre in a reflection-mode configuration without requiring
access to the distal facet. A thin stack of structured metasurface reflectors
is used at the distal facet of the fibre to introduce wavelength-dependent,
spatially heterogeneous reflectance profiles. We derive a first-order fibre
model that compensates these wavelength-dependent changes in the fibre
transmission matrix and show that, consequently, the reflected data at 3
wavelengths can be used to unambiguously reconstruct the full transmission
matrix by an iterative optimisation algorithm. We then present a method for
sample illumination and imaging following reconstruction of the transmission
matrix. Unlike previous approaches, our method does not require the fibre
matrix to be unitary making it applicable to physically realistic fibre
systems that have non-negligible power loss. We demonstrate the transmission
matrix reconstruction and imaging method first using simulated non-unitary
fibres and noisy reflection matrices, then using much larger experimentally-
measured transmission matrices of a densely-packed multi-core fibre. Finally,
we demonstrate the method on an experimentally-measured multi-wavelength set
of transmission matrices recorded from a step-index multimode fibre. Our
findings pave the way for online transmission matrix calibration _in situ_ in
hair-thin optical fibres.

optical fibres $|$ imaging $|$ endoscopy $|$ metasurfaces $|$ nanostructures
$|$ polarisation

††preprint: APS/123-QED

##  I Introduction

Lensless imaging through hair-thin optical fibres is a technique that promises
to open up new areas in biomedical imaging such as: _in vivo_ bright-field,
dark-field, and fluorescence microscopy Cižmár and Dholakia (2012);
quantitative phase and polarimetric imaging with applications in early
detection of cancer Gataric  _et al._ (2019); Gordon  _et al._ (2018); and
endoscopic confocal microscopy for high-resolution imaging Loterie  _et al._
(2015). For all these applications, accurate characterisation of the
deterministic propagation of light through the fibre, discretised as the
transmission matrix (TM) Rotter and Gigan (2017), is essential for imaging,
whether used directly Gordon  _et al._ (2018), indirectly via phase
conjugation Dunning and Lind (1982) or via optimisation Di Leonardo and
Bianchi (2011).

Unfortunately, the TM is highly sensitive to small changes in mechanical
deformation and temperature that are unavoidable in living subjects, so the TM
should be calibrated immediately before imaging. The most direct way of
measuring a TM is to send pre-defined optical fields into the distal facet of
the fibre (furthest from the operator) and measure the resulting fields at the
proximal facet (nearest the operator). This widely-used approach has a
critical limitation: the distal facet of the fibre is deep within the subject
during imaging so the pre-defined optical fields cannot be reliably generated
without additional bulky distal optics that would compromise the ultra-thin
form factor.

Several methods have been proposed to overcome this limitation and open new
avenues for minimally invasive optical imaging in inaccessible areas of the
body, e.g. deep in the brain where the fibre may need to curve for safe access
Vasquez-Lopez  _et al._ (2018). Placing a holographic plate on the distal
facet of the fibre, which is illuminated via a physically separate single-mode
fibre, can be used to create a ‘virtual beacon’ Farahi  _et al._ (2013). A
phase-conjugated version of the resultant proximal field is used to recreate a
single focussed spot at the distal facet. This spot is typically fixed in
position but with a multi-core fibre (MCF), the ‘memory effect’ could be
exploited to enable scanning over a small distances Stasio  _et al._ (2015).
This method is only able to partially recover the TM and so while confocal
imaging is possible, wide-field imaging (e.g. quantitative phase) is not
achievable. Further, the addition of a single-mode fibre to illuminate the
beacon adds bulk to the system.

Highly accurate modelling of TM perturbations has also been shown to be
feasible, but requires precision manufactured fibres such as
telecommunications OM4 grade multi-mode fibre (MMF), is limited to short
distances ($<$10cm) Plöschner  _et al._ (2015) and needs very precise bending
characterisation, for example through the addition of a Bragg-grating shape
sensing fibre Spillman and Udd (2014); Udd (2011). Another recently
demonstrated method uses reflectors at the distal end to introduce time delays
to resolve the fibre TM (Chen  _et al._ , 2018). However, this adds
undesirable complexity to the system and also requires a unique reflector for
each characterised propagation mode, which does not scale well for high-
resolution imaging applications.

Gu et al. propose a more generic approach to _in situ_ TM characterisation: to
infer the forward TM based on the reflection matrix (RM) i.e. light that has
made a round-trip into the proximal facet, off a distal reflector and back out
the proximal facet Gu  _et al._ (2015). Applied naively, this approach suffers
two major limitations. The first arises due to the transpose symmetry of fibre
TMs, which means that measured RMs exhibit a quadratic relationship with their
respective TMs, resulting in ambiguities during recovery Lee  _et al._ (2018);
Gu  _et al._ (2015). In some cases this ambiguity reduces to an
$N$-dimensional vector of sign errors (i.e. a vector $\in\\{-1,1\\}^{N}$) for
an $N$-pixel image, for example, when using MCF with low core-to-core coupling
Warren  _et al._ (2016) or using an MMF that is near perfectly unitary (via
the Autonne-Takagi matrix factorisation) Gu  _et al._ (2015). In such cases,
the ambiguity may be resolved by _in situ_ optimisation assuming prior
knowledge of the sample Warren  _et al._ (2016) or by treating the ambiguities
as invariants of the fibre that can be measured in advance Gu  _et al._
(2015). In practical cases, however, TMs are not unitary even in precision-
manufactured MMF Gu (2017); Carpenter  _et al._ (2014) because high-order
modes, with more power concentrated in the fibre cladding, tend to exhibit
power leakage at bends Gordon  _et al._ (2014).

The second limitation arises because light must be able to pass through the
distal facet of the fibre to enable imaging. Gu et al. address this by
proposing to include a shutter at the distal facet, but again, this added bulk
would compromise the ultra-thin form factor Gu  _et al._ (2015).

Seeking to overcome the need for precision-manufactured MMFs, while
maintaining the hair-thin property of the fibres, we demonstrate here a new
method that enables measurement of non-unitary fibre TMs of arbitrary size
without access to the distal facet and without adding distal bulk, making it
suitable for ultra-thin flexible imaging devices. The key innovation is a
distal reflector comprising a multi-layer stack of spatially heterogeneous
partially reflecting wire-grid polarisers (termed ‘metasurfaces’) and long-
pass optical filters. The filter stack creates a spatially and spectrally
varying reflector, whose properties can be characterised prior to use and will
remain fixed throughout any imaging experiment. By combining multiple RMs
recorded at different wavelengths it is possible to estimate the instantaneous
TM and then perform imaging at another wavelength.

We first introduce a new method that analytically extracts the TM based on 3
measured RMs at different wavelengths and prior knowledge of the reflector
properties (the ‘zeroth-order method’). We then extend the physical model to
account for realistic path-length differences and present an iterative
numerical method that is able to correct analytically derived TMs to satisfy
this extended model (the ‘first-order method’). Neither method requires that
the TM be unitary, only that it is invertible, meaning they could be applied
to lossy fibres or general scattering media. Following TM reconstruction, we
then show how imaging can be performed at a fourth wavelength.

Next, we computationally demonstrate the first-order method using simulated
32$\times$32 non-unitary fibre matrices with noise. We then demonstrate the
recovery algorithm using measured MCF TMs ($1648\times 1648$) and
realistically simulated reflector stacks, showing the method can be scaled up
effectively. Finally, we demonstrate TM recovery using a multi-wavelength set
of TMs measured from a step-index MMF validating both that the method works
for real MMF and that the first-order physical model used to predict changes
in TM with wavelength is physically realistic. Our findings underpin a
compelling new approach for enabling ultra-thin flexible lensless imaging via
optical fibres.

##  II Theory

###  II.1 Formulation of the physical model

The physical model used for reference in this paper is shown in Figure 1. We
consider a well-established model where a field containing $M$ pixels, each
complex-valued to encode amplitude and phase, in two orthogonal polarisations
is recorded. The field at some input plane (e.g. the proximal facet), is
represented by a vector $\mathbf{x}\in\mathbb{C}^{2M}$, and the field at an
output plane (e.g. the distal facet) is represented by a vector
$\mathbf{y}^{\prime}\in\mathbb{C}^{2M}$. In the forward propagation direction,
these vectors are related by the fibre TM, $\mathbf{A}\in\mathbb{C}^{2M\times
2M}$:

$\mathbf{y}^{\prime}=\mathbf{A}\mathbf{x}$ (1)

This equation can be derived by discretisation of the linear propagation
relationships between the two planes Rotter and Gigan (2017); Fan and Kahn
(2005). For example, a continuous optical field might be sampled using delta
functions at $M$ spatial locations in two orthogonal polarisation states using
a Jones vector formalism Goldstein (2010); Gataric  _et al._ (2019). The
complex-valued samples, representing amplitude and phase across two
polarisations, can then be ordered to produce the vector $\mathbf{x}$ with
$2M$ elements Gordon  _et al._ (2018).

Figure 1: Physical model used for fibre TM characterisation in reflection mode
and for imaging. The proximal facet of the fibre is nearest to the imaging
system (i.e. outside of the subject) while the distal facet is nearest the
object being imaged (i.e. inside the subject). Vectors at the distal facet are
denoted using the prime symbol (′), vectors at the exit facet of the reflector
stack are denoted with double prime (′′), and vectors at the imaging plane are
denoted with triple prime (′′′). $\mathbf{A}$ represents the fibre TM,
$\mathbf{R}$ the PM of the reflector stack, $\mathbf{C}$ the round-trip RM
(down the fibre, off the reflector, back up the fibre), $\mathbf{A}_{refl}$
the TM of the reflector stack, $\mathbf{G}$ a free-space propagation operator
(parameterised by the distance, $d$, between the reflector structure and the
imaging plane), and $\mathbf{r}_{target}$ the reflectance profile of the
sample.

Similarly, if we consider a field, $\mathbf{x}^{\prime}\in\mathbb{C}^{2M}$, at
the distal facet propagating in the reverse direction to become a field,
$\mathbf{y}\in\mathbb{C}^{2M}$, at the proximal facet, these are related by
the transpose of the fibre TM Gu  _et al._ (2015):

$\mathbf{y}=\mathbf{A}^{\top}\mathbf{x}^{\prime}$ (2)

We then consider the impact of adding a reflector at the distal facet of the
fibre. In general, the reflector is considered to be spatially heterogeneous
in terms of its localised Jones reflection matrices: there may be uncorrelated
Jones matrices describing reflections at each spatial point. Further, if the
reflector is offset from the fibre, light may couple between spatial positions
due to diffraction. This behaviour is linear and so is represented by a
partial reflector matrix (PM), $\mathbf{R}\in\mathbb{C}^{2M\times 2M}$, that
relates $\mathbf{y}^{\prime}$ and $\mathbf{x}^{\prime}$ (see Figure 1):

$\mathbf{x}^{\prime}=\mathbf{R}\mathbf{y}^{\prime}$ (3)

Next, combining Equations 1, 2 and 3 we can determine the optical field
exiting the proximal facet in the reverse direction:

$\mathbf{y}=\mathbf{A}^{\top}\mathbf{R}\mathbf{A}\mathbf{x}=\mathbf{C}\mathbf{x}$
(4)

where $\mathbf{C}=\mathbf{A}^{\top}\mathbf{R}\mathbf{A}\in\mathbb{C}^{2M\times
2M}$ is termed a ‘reflection matrix’ (RM). Physically speaking, $\mathbf{C}$
represents light taking a complete round-trip (or double-pass): forwards down
the fibre, off a given reflector and back up the fibre, as illustrated in
Figure 1. When $\mathbf{C}$ is measured in practice, the measurements may lead
to a non-square approximation to $\mathbf{C}$, denoted as
$\tilde{\mathbf{C}}$, hence $\tilde{\mathbf{C}}$ first needs to be downsampled
into a square matrix $\tilde{\tilde{\mathbf{C}}}$, which is then used as a
surrogate for $\mathbf{C}$ in what follows (see Appendix A.2).

As our first objective is to estimate the TM of the fibre $\mathbf{A}$ we
would ideally use pairs of measured vectors,
$(\mathbf{x},\mathbf{y}^{\prime})$ through standard a calibration procedure
Gordon  _et al._ (2018); Gataric  _et al._ (2019); Carpenter  _et al._ (2014).
Since we now only have access to the proximal facet of the fibre, we instead
need to estimate the RM, $\mathbf{C}=\mathbf{A}^{\top}\mathbf{R}\mathbf{A}$,
based on the measured vectors $(\mathbf{x},\mathbf{y})$. Given that we can
characterise $\mathbf{R}$ in advance (Appendix A.3.1) and it will remain fixed
throughout use, the goal is then to recover $\mathbf{A}$ based on measurements
of $\mathbf{C}$ and prior characterisations of $\mathbf{R}$, however, with
only this information it is not in general possible to recover $\mathbf{A}$
unless the fibre TM is unitary Gu  _et al._ (2015); Gu (2017). We therefore
propose a novel approach that uses several different reflectors, with PMs
$\mathbf{R}_{n},~{}n=1..Q$, and show that this enables unambiguous recovery of
$\mathbf{A}$ in the more general non-unitary case.

For experimental simplicity, it would be desirable to use as few different
reflectors, $\mathbf{R}_{n},~{}n=1..Q$, and associated RMs,
$\mathbf{C}_{n},~{}n=1..Q$, as possible for TM characterisation and to switch
between these without any physical moving parts, for example, by modulating
wavelength. The case with a single reflector reduces to the method presented
by Gu et al. Gu  _et al._ (2015) and does not produce unique solutions for
non-unitary TMs, $\mathbf{A}$. Therefore, we need $Q\geq 2$ reflectors.

We will consider in what follows the case of three reflectors, i.e. $Q=3$. For
TM characterisation (with reference to the appropriate optical path in Figure
1) there will be three associated PMs $\mathbf{R}_{n},~{}n=1..3$. The TMs
associated with each reflector during characterisation,
$\mathbf{A}_{n},n=1..3$, may in general be different, for example if
wavelength modulation is used to activate different reflectors. We can then
record 3 associated RMs, $\mathbf{C}_{n},~{}n=1..3$ measured for example using
the process of Appendix A.3.2. Based on Equation 4, we express these RMs in
terms of the PMs $\mathbf{R}_{n},~{}n=1..3$ and associated TMs
$\mathbf{A}_{n},~{}n=1..3$ as:

$\mathbf{C}_{1}=\mathbf{A}_{1}^{\top}\mathbf{R}_{1}\mathbf{A}_{1}$ (5)

$\mathbf{C}_{2}=\mathbf{A}_{2}^{\top}\mathbf{R}_{2}\mathbf{A}_{2}$ (6)

$\mathbf{C}_{3}=\mathbf{A}_{3}^{\top}\mathbf{R}_{3}\mathbf{A}_{3}$ (7)

In order to solve these equations for $\mathbf{A}_{n},n=1..3$ we must assume a
known relationship between these TMs, for example a wavelength-dependent phase
shift, and must also have measured $\mathbf{R}_{1}$, $\mathbf{R}_{1}$ and
$\mathbf{R}_{1}$ _a priori_ using a characterisation process such as that in
Appendix A.3.1.

Following TM recovery, the final objective is to performing imaging. This
involves using an extended optical path (Figure 1) in which light propagates
down the fibre via a possibly different fourth TM, $\mathbf{A}_{4}$, is
partially reflected off the reflector via a PM, $\mathbf{R}_{4}$, and is
partially transmitted through the reflector via a TM, $\mathbf{A}_{refl}$.
$\mathbf{R}_{4}$ and $\mathbf{A}_{r}efl$ can be measured _a priori_ using a
characterisation process such as that in Appendix A.3.1. Following possible
further propagation and reflection off the sample, the returned light at the
proximal facet is used for image reconstruction.

Again, a known physical relationship is required to recover $\mathbf{A}_{4}$
from $\mathbf{A}_{1}$, $\mathbf{A}_{2}$ or $\mathbf{A}_{3}$ and hence
reconstruct images. We therefore consider two simple models to approximate
this relationship between the calibration and imaging TMs and enable image
recovery: a ‘zeroth-order model’ and a ‘first-order model’.

###  II.2 Zeroth-order model

The zeroth-order model assumes that the TMs under the different reflectors are
the same, i.e.:

$\mathbf{A}=\mathbf{A}_{1}=\mathbf{A}_{2}=\mathbf{A}_{3}=\mathbf{A}_{4}$ (8)

This could correspond to physically switching reflectors without disturbing
the fibre TM. Alternatively, if modulating wavelength to switch between
reflectors, this corresponds to the assumption that the fibre TM does not
change with wavelength. Using this model, the TM can be recovered in a
computationally straightforward way relying largely on analytical steps
(Appendix B), providing a number of key insights. First, if the number of
reflectors is $n\geq 3$ then analytical solutions for the TMs can be
recovered. Numerical and analytical investigations for the case with $n=2$
reflectors were inconclusive. Second, the eigenvalues of each PM,
$R_{n},n=1..3$ must be distinct for unambiguous TM recovery, in agreement with
previous findings Gu  _et al._ (2015). This means that light must be coupled
between different modes and polarisations and so, for example, a conventional
partial mirror reflector would not work because it preserves polarisation
states leading to repeated matrix eigenvalues. This insight influences the
potential design of any reflector.

Though this model may be used if reflectors are physically switched, for the
wavelength modulation approach, Equation 8 is physically unrealistic because
it neglects wavelength-induced phase-shifts that are significant even over
small physical distances (see Appendix E). However, even if these phase-shifts
are neglected the recovered TMs still provide an approximate starting point
for more complex models relating TMs at different wavelengths.

###  II.3 First-order model

Considering again the idea of switching reflectors by modulating wavelength, a
more physically accurate relationship between TMs is constructed by
considering linear phase shift as wavelength is varied. In general the
relationship between fibre TMs at different wavelengths is complex Redding
_et al._ (2013) so this linear approximation is limited to be within the
‘spectral bandwidth’ of the fibre. We can then use a coupled-mode theory
treatment of optical fibres Yariv (1973); Huang (1994). This perturbational
approach models a length of MMF as a sequence of infinitesimal segments that
introduce field coupling. The differential change in field in each mode can be
modelled by a vector, $d\mathbf{m}\in\mathbb{C}^{2M\times 1}$, such that:

$d\mathbf{m}=d\mathbf{A}\mathbf{m}$ (9)

where $d\mathbf{A}\in\mathbb{C}^{2M\times 2M}$ represents an infinitesimal
coupling matrix, and $\mathbf{m}\in\mathbb{C}^{2M\times 1}$ represents the
field in each of $M$ spatial modes over 2 orthogonal polarisations. This
matrix equation represents a system of $2M$ linear differential equations. The
solution to these equations for any arbitrary input set of modes,
$\mathbf{x}$, after travelling a distance $\ell$ along the central axis of the
fibre is given by a matrix exponential Yariv (1973); Loterie  _et al._ (2017):

$\mathbf{y}=e^{d\mathbf{A}\ell}\mathbf{x}$ (10)

Equivalently, this model can be thought of as the Lie algebra $d\mathbf{A}$
constructing the Lie group of invertible matrices
$e^{d\mathbf{A}\ell}\in\mathbb{C}^{2M\times 2M}$ Dragt (1982). We can estimate
$d\mathbf{A}$ from the TM at characterisation wavelength $\lambda_{1}$,
$\mathbf{A}_{1}$ as:

$d\mathbf{A}=\frac{\log\mathbf{A}_{1}}{\ell_{1}}$ (11)

where $\log\mathbf{A}_{1}$ represents a matrix logarithm and $\ell_{1}$ is the
optical path length along the central fibre axis. In general, matrix
logarithms of complex matrices produce degenerate solutions, analogous to
logarithms of complex numbers. This means that in estimating $d\mathbf{A}$
each eigenvalue, $\gamma$, has an equivalence class of degenerate forms:
$\gamma+i2\pi n/\ell_{1}$, $n\in\mathbb{Z}$. However, under the bandwidth-
limited scenario presented here, the change in wavelength results in an
overall relative phase shift between any two matrix elements of $<2\pi$
(Appendix E) and so therefore $n=0$ is the appropriate choice.

If we assume constant dispersion within the relatively small ($<10$nm)
spectral bandwidth considered, a reasonable assumption for typical glasses
Payne and Gambling (1975), we can find the equivalent optical thicknesses at
the other characterisation wavelengths ($\lambda_{2}$, $\lambda_{3}$) and
imaging wavelength ($\lambda_{4}$) relative to the first characterisation
wavelength ($\lambda_{1}$):

$\ell_{2}=\frac{\lambda_{1}}{\lambda_{2}}\ell_{1}$ (12)

$\ell_{3}=\frac{\lambda_{1}}{\lambda_{3}}\ell_{1}$ (13)

$\ell_{4}=\frac{\lambda_{1}}{\lambda_{4}}\ell_{1}$ (14)

Rearranging and expressing in terms of $\mathbf{A}_{1}$ we can write

$\mathbf{A}_{2}=e^{d\mathbf{A}\ell_{2}}=e^{d\mathbf{A}\frac{\lambda_{1}}{\lambda_{2}}\ell_{1}}=e^{\frac{\lambda_{1}}{\lambda_{2}}\log\mathbf{A}_{1}}$
(15)

$\mathbf{A}_{3}=e^{\frac{\lambda_{1}}{\lambda_{3}}\log\mathbf{A}_{1}}$ (16)

$\mathbf{A}_{4}=e^{\frac{\lambda_{1}}{\lambda_{4}}\log\mathbf{A}_{1}}$ (17)

A key requirement for this model is that the wavelength sweep range is kept
within the spectral bandwidth of the fibre, i.e. the bandwidth within which
far-field speckle patterns from a fibre remain correlated Mosk  _et al._
(2012) or equivalently over which the eigenmodes of the fibre remain
approximately the same. We show in Appendix E that this implies the maximum
change in path length for any mode is $<2\pi$.

The eigendecomposition of $d\mathbf{A}$ is:

$d\mathbf{A}=\mathbf{Q}\left[\begin{array}[]{c c c}\Lambda_{1}&\cdots&0\\\
\vdots&\ddots&\vdots\\\
0&\cdots&\Lambda_{2M}\end{array}\right]\mathbf{Q}^{-1}$ (18)

where $\mathbf{Q}\in\mathbb{C}^{2M\times 2M}$ is a matrix whose columns,
$\mathbf{q}_{m}$ with $m=1..2M$, are the eigenvectors of $d\mathbf{A}$ and
$\Lambda_{m}$ with $m=1..2M$ are the corresponding eigenvalues. The
eigendecomposition of $\mathbf{A}_{n}$ with $n=2,3,4$ is therefore Hall
(2015):

$\mathbf{A}_{n}=e^{\frac{\lambda_{1}}{\lambda_{n}}\log\mathbf{A}_{1}}=\mathbf{Q}\left[\begin{array}[]{c
c c}e^{\frac{\lambda_{1}}{\lambda_{n}}\Lambda_{1}}&\cdots&0\\\
\vdots&\ddots&\vdots\\\
0&\cdots&e^{\frac{\lambda_{1}}{\lambda_{n}}\Lambda_{2M}}\end{array}\right]\mathbf{Q}^{-1}$
(19)

In general, matrices constructed using matrix exponentials (as in Equations 15
– 17) share the same eigenvectors but the corresponding eigenvalues are
related by an exponential Hall (2015). Our model therefore assumes that the
eigenvectors of the fibre do not change appreciably over the selected
bandwidth, which is consistent with previous experimental work Carpenter  _et
al._ (2015). This is in contrast to the zeroth-order model that requires the
additional assumption that eigenvalues stay constant within the spectral
bandwidth.

It is not straightforward to solve Equation 5 to Equation 7 analytically to
obtain $\mathbf{A}_{1}$, $\mathbf{A}_{2}$, $\mathbf{A}_{3}$ or
$\mathbf{A}_{4}$ when they are related by matrix exponentials as in Equations
15 – 17. Therefore, we present an optimisation-based approach. We begin by
using the zeroth-order approximation to obtain an approximate solution for
$\mathbf{A}_{1}$, termed $\mathbf{\hat{A}}_{1}$. This is achieved by
performing an optimisation based on Equation 5:

$\mathbf{\hat{A}}_{1}^{\prime}=\arg\min_{\mathbf{A}}\left\lVert\mathbf{A}^{\top}\mathbf{R}_{1}\mathbf{A}-\mathbf{C}_{1}\right\rVert_{F}$
(20)

where $\left\lVert..\right\rVert_{F}$ indicates the Frobenius norm. The
solution to this optimization problem, $\mathbf{\hat{A}}^{\prime}_{1}$, may be
found using an iterative gradient descent solver initialised by the zeroth-
order approximation. We know that because the true $\mathbf{A}_{1}$ may be
non-unitary, this equation may not have a unique solution Gu  _et al._ (2015).
We therefore use the result of Equation 20, $\mathbf{\hat{A}}_{1}^{\prime}$,
to generate an estimate of $\mathbf{A}_{2}$, as
$\mathbf{\hat{A}}_{2}=e^{\frac{\lambda_{1}}{\lambda_{2}}\log(\mathbf{\hat{A}}_{1}^{\prime})}$.
This is then used as the starting point for a second optimisation:

$\mathbf{\hat{A}}_{2}^{\prime}=\arg\min_{\mathbf{A}}\left\lVert\mathbf{A}^{\top}\mathbf{R}_{2}\mathbf{A}-\mathbf{C}_{2}\right\rVert_{F}$
(21)

Finally,
$\mathbf{\hat{A}}_{3}=e^{\frac{\lambda_{2}}{\lambda_{3}}\log(\mathbf{\hat{A}}_{2}^{\prime})}$
is used as the starting point for a third optimisation:

$\mathbf{\hat{A}}_{3}^{\prime}=\arg\min_{\mathbf{A}}\left\lVert\mathbf{A}^{\top}\mathbf{R}_{3}\mathbf{A}-\mathbf{C}_{3}\right\rVert_{F}$
(22)

To complete the iteration,
$\mathbf{\hat{A}}_{1}=e^{\frac{\lambda_{3}}{\lambda_{1}}\log(\mathbf{\hat{A}}_{3}^{\prime})}$
is used as a starting point for Equation 20, and the process is repeated. For
numerous simulated and experimentally measured TMs we observe that
$\mathbf{\hat{A}}_{1}$ converges to the true value of $\mathbf{A}_{1}$,
thereby enabling calculation of $\mathbf{A}_{4}$, required for image
reconstruction. This process is summarised in the flow-chart of Figure 2. In
this work, for simulated TMs the convergence criteria is a fixed number of
iterations (both 300 and 500 are examined) while for experimentally measured
TMs the convergence criteria is an error gradient (terminate when change in
error is $<$1% of value at first iteration).

Figure 2: Flowchart detailing the iterative optimisation routine to estimate
the TM required for imaging, $\mathbf{A}_{4}$, using the first-order model.

Each of Equation 20–22 is in itself degenerate due to the quadratic matrix
form $\mathbf{A}^{\top}\mathbf{R}\mathbf{A}$ (see Section II.1). However, this
iterative process may avoid local minima because such degeneracies are not
typically the same for all three equations. The exception to this is a global
sign error: $\pm\mathbf{A}_{n}$ with $n=1..3$, is optimal for each respective
equation. However, this degeneracy corresponds to a global phase-shift (i.e. a
constant factor) so can be neglected for almost all imaging modalities.

Local optima can still arise when the optimisation of Equation 22 exactly
reverses the optimisations of Equations 20 and 21. It is observed empirically
that using the zeroth-order approximation as a starting point tends to avoid
this and produces reliable convergence to the global minimum.

The algorithm presented here is observed to converge reliably to the correct
solution and has the computational advantage that the relatively expensive
matrix logarithm operation is performed only between minimisations, as opposed
to in each objective function. It is therefore used for the remainder of this
paper.

###  II.4 Imaging

Considering the case in which reflectors are switched by modulating wavelength
and the first-order model is used to reconstruct the fibre TM at the imaging
wavelength ($\lambda_{4}$), $\mathbf{A}_{4}$, the next step is to reconstruct
an image of a sample using data recorded at the imaging wavelength. With
reference to the physical model (Figure 1) this requires prior measurement of
the TM and PM of the reflector, $\mathbf{A}_{refl}$ and $\mathbf{R}_{4}$
respectively, at the imaging wavelength, $\lambda_{4}$, via a calibration
process (e.g. Appendix A.3.1). These quantities are assumed to stay constant
throughout operation. The imaging process is summarised in the flow chart of
Figure 3. An experimental system in which such a workflow could be realised is
described in Appendix A.1, where a spatial light modulator (SLM) is used to
project optical fields onto the proximal fibre facet.

Figure 3: Flowchart detailing how an image of a sample is recorded.

There are two approaches to illuminating the sample via the imaging fibre. The
first is to use the measured TM at the imaging wavelength, $\mathbf{A}_{4}$,
to produce a known distal illumination profile. This produces controllable
illumination profiles that could enable, for example, various structured
illumination techniques but assumes that the TM reconstruction process is
complete before imaging. In cases where computing the TM introduces an
undesirable time delay, a known proximal illumination profile could used to
produce a semi-random distal illumination profile that can later be corrected
by averaging over several random illuminations, a technique known as speckle
averaging Choi  _et al._ (2012), or by using knowledge of $\mathbf{A}_{4}$
reconstructed offline.

Considering our physical model (Figure 1), we project a random or pre-
calculated illumination vector, $\mathbf{x}_{illum}$, onto the proximal facet
of the fibre, giving at the distal facet:

$\mathbf{y}_{illum}^{\prime}=\mathbf{A}_{4}\mathbf{x}_{illum}$ (23)

Next, the light passes through the full reflector structure, which in the
imaging configuration must be partially transmissive. In general
$\mathbf{A}_{refl}\in\mathbb{C}^{2N\times 2M}$ with $N\geq M$, recalling that
$\mathbf{A}_{4}\in\mathbb{C}^{2M\times 2M}$. This allows for oversampled
illumination fields. The field exiting the reflector stack is then:

$\mathbf{y}^{\prime\prime}_{illum}=\mathbf{A}_{refl}\mathbf{A}_{4}\mathbf{x}_{illum}$
(24)

with $\mathbf{y}^{\prime\prime}_{illum}\in\mathbb{C}^{2N}$. This can be
propagated by some linear operator, $\mathbf{G}\in\mathbb{C}^{2N\times 2N}$,
through free-space (e.g. Fresnel or Fraunhofer propagation) before reaching
the target where the field will be:

$\mathbf{y}_{illum}^{\prime\prime\prime}=\mathbf{G}\mathbf{A}_{refl}\mathbf{A}_{4}\mathbf{x}_{illum}$
(25)

In general $\mathbf{G}$ is parameterised by the distance, $d$, between the
subject and the distal surface of the reflector structure (Figure 1. In
general, $d$ is not known _a priori_ and so must be estimated during
operation, for example using an automatic focus detection algorithm such as
the Brenner algorithm Brenner  _et al._ (1976); Gordon  _et al._ (2018).

For simplicity we will assume the target is a reflective surface, though a
volume scattering target could be characterised using several known
illumination profiles in sequence as in spatial frequency domain imaging
Cuccia  _et al._ (2009). We therefore let the target be represented by vector,
$\mathbf{r}_{target}\in\mathbb{C}^{2N}$. An element-wise (Hadamard) product
with the illumination gives:

$\mathbf{x}_{target}^{\prime\prime\prime}=\mathbf{y}_{illum}^{\prime\prime\prime}\circ\mathbf{r}_{target}$
(26)

This in turn, produces a field at the distal fibre facet of:

$\mathbf{x}_{target}^{\prime}=\mathbf{A}_{refl}^{\top}\mathbf{G}^{\top}\mathbf{x}_{target}^{\prime\prime\prime}$
(27)

Next, we consider the illumination light that is reflected back from the
reflector stack:

$\mathbf{x}_{refl}^{\prime}=\mathbf{R}_{4}\mathbf{A}_{4}\mathbf{x}_{illum}$
(28)

We sum the two fields and propagate back through the fibre giving:

$\mathbf{y}_{total}=\mathbf{A}_{4}^{\top}(\mathbf{x}_{target}^{\prime}+\mathbf{x}_{refl}^{\prime})$
(29)

The goal is to recover $\mathbf{r}_{target}$ from $\mathbf{y}_{total}$, the
raw measured data. To do this, we first recover
$\mathbf{x}_{target}^{\prime\prime\prime}$ by rearrangement of Equation 29 and
substitution of Equations 27 and 28:

$\left(\mathbf{A}_{4}^{\top}\right)^{-1}\mathbf{y}_{total}=(\mathbf{x}_{target}^{\prime}+\mathbf{x}_{refl}^{\prime})$

$\mathbf{x}_{target}^{\prime}=\left(\mathbf{A}_{4}^{\top}\right)^{-1}\mathbf{y}_{total}-\mathbf{R}_{4}\mathbf{A}_{4}\mathbf{x}_{illum}$

$\mathbf{x}_{target}^{\prime\prime\prime}=\left(\mathbf{A}_{refl}^{\top}\mathbf{G}^{\top}\right)^{+}\left(\left(\mathbf{A}_{4}^{\top}\right)^{-1}\mathbf{y}_{total}-\mathbf{R}_{4}\mathbf{A}_{4}\mathbf{x}_{illum}\right)$
(30)

where $(..)^{+}$ represents a generalised inverse because $\mathbf{A}_{refl}$
may be non-square. Quantities on the RHS of Equation 30 are known prior to
imaging, except for $\mathbf{y}_{total}$ which is recorded during imaging and
possibly $\mathbf{G}$, which may need to be estimated as described above.
Therefore, we could naively recover an estimate of $\mathbf{r}_{target}$,
which we term $\mathbf{\tilde{r}}_{target}$, from
$\mathbf{x}_{target}^{\prime\prime\prime}$ through an element-wise (Hadamard)
division:

$\mathbf{\tilde{r}}_{target}=\mathbf{x}_{target}^{\prime\prime\prime}\oslash\mathbf{y}_{illum}^{\prime\prime\prime}$

$\mathbf{\tilde{r}}_{target}=\mathbf{x}_{target}^{\prime\prime\prime}\oslash\left(\mathbf{G}\mathbf{A}_{refl}\mathbf{A}_{4}\mathbf{x}_{illum}\right)$

However, this approach may erroneously amplify noise in areas where
illumination power, $|\mathbf{y}_{illum}^{\prime\prime\prime}|^{2}$, is low
producing speckle-like artefacts in the recovered
$\mathbf{\tilde{r}}_{target}$. To compensate for these artefacts we apply a
probabilistic approach. Using the model derived in Section III.6, the noise at
pixel $(x,y)$ is considered circularly-symmetric complex Gaussian distributed
with standard deviation $\sigma$. The noise power at that pixel, termed
$\mathrm{n}(x,y)$, is therefore drawn from a $\chi^{2}$ distribution with 1
degree of freedom:

$\frac{\mathrm{n}(x,y)}{\sigma^{2}}\sim\chi^{2}(1)$

Using an empirically measured value for $\sigma$ we use the $\chi^{2}$
cumulative distribution function, $\chi^{2}_{CDF}$, to estimate for each pixel
a probability that it is generated by noise, $p_{n}(x,y)$:

$p_{n}(x,y)=\chi^{2}_{CDF}\left(\frac{||\mathbf{x}_{target}^{\prime\prime\prime}(x,y)||^{2}}{\sigma^{2}}\right)$

We can order $p_{n}(x,y)$ across all $(x,y)$ into a vector denoted as
$\mathbf{p}_{n}$. Separating the illumination into amplitude and phase parts
as:

$\mathbf{a}_{illum}=|\mathbf{y}_{illum}^{\prime\prime\prime}|$

$\mathbf{\theta}_{illum}=\arg(\mathbf{y}_{illum}^{\prime\prime\prime})$

we then use element-wise division by a probabilistically weighted term as
follows:

$\mathbf{\tilde{r}}_{target}=\mathbf{x}_{illum}^{\prime\prime\prime}\oslash\left(\left(\left(1-\mathbf{p}_{n}\right)\circ\mathbf{a}_{illum}+\mathbf{p}_{n}\right)\circ
e^{i\mathbf{\theta}_{target}}\right)$ (31)

With this approach, pixels with smaller illumination amplitudes tend only to
modify the phase of $\mathbf{x}_{illum}^{\prime\prime\prime}$ when correcting
the speckle-like artefacts during recovery of $\mathbf{\tilde{r}}_{target}$.
Pixels with larger illumination amplitudes, more likely to contain useful
information, experience both amplitude and phase correction. In this way, the
speckle-like artefacts in $\mathbf{\tilde{r}}_{target}$ are reduced. Areas of
low illumination will still have poor signal-to-noise ratio but this could be
corrected by using a more uniform pre-determined illumination profile or,
again, by implementing speckle averaging Choi  _et al._ (2012).

##  III Methods

###  III.1 Implementation of reflector

The switchable reflector structure described in Section II.1 could be
implemented in a number of ways. Reflectors could be switched using a physical
mechanism, e.g by translating or rotating a large distal plate relative to the
facet. Though this would enable use of the computationally simpler zeroth-
order reconstruction approach, it would require actuators and a large plate at
the distal facet, adding undesired bulk and compromising the ultra-thin form
factor. Here, we propose instead that switching be achieved through modulating
wavelength using the concept of a reflector stack (Figure 4a) with different
reflectance behaviour at several closely spaced wavelengths.

The envisaged stack comprises layers of long-pass optical filters so that
light penetrates to different depths at different wavelengths (Figure 4b).
Though the first-order reconstruction approach must be used on account of the
wavelength modulation, it is still advantageous to satisfy the requirements
for the zeroth-order model so that it can be used as a starting point for
optimisation and because it is a special case of the first-order model.
Therefore, we require that the stack should couple light between different
modes of the fibre, including polarisation modes, in a pseudo-random manner to
ensure with high probability that the eigenvalues of the TM are distinct as
required (Section II.2). This could be realised by placing atop each filter a
thin (20–50nm) layer comprising spatially heterogenous wire-grid polarisers
(Figure 4c) fabricated, for example, using elctron beam lithography Williams
_et al._ (2016). To act as optical polarisers the fabricated wires must have
width and pitch less than the wavelength of light ($<\lambda/2$) and so these
layers are classed ‘metasurfaces’. Spatial heterogeneity would be achieved
using arrays of partial polariser cells that have varying linear diattenuation
(or dichroism) and diattenuation axis orientations (Figure 4c).

Figure 4: Design of a wavelength-dependent reflector stack: a) Proposed stack
structure, attached to distal facet of optical fibre, comprising layers of
glass and long-pass optical filters with wire-grid polariser metasurfaces
sandwiched between. Light paths are shown for the 4 different wavelenths used:
the characterisation wavelengths $\lambda_{1}$, $\lambda_{2}$ and
$\lambda_{3}$ each producing 3 reflector matrices, $\mathbf{R}_{n}~{}n=1..3$;
and the imaging wavelength $\lambda_{4}$ that passes the whole way through the
stack. The stack behaviour can be understood by considering optical fields
exiting the fibre, propagating through each layer via diffraction and being
partially reflected at the metasurfaces generating backward propagating
wavefronts. The sum of all backward propagating components represents the
field coupled back into the fibre. b) Spectral transmission profiles of the 3
long-pass optical filters with example imaging and characterisation
wavelengths labelled. c) Example design of wire-grid polariser metasurface
comprising an array of diattenuating cells of size 1–2$\mu$m containing
nanowire gratings with varying pitch (50–200nm), duty cycle (25%–75%) and
orientation (0–2$\pi$). The pitch and duty cycle determine the linear
diattenuation of each cell, $D$, and the orientation determines the
diattenuation axis, $\theta_{D}$. These are used to form a Jones matrix for
each cell that determines its optical transmission and reflection.

Because of the optical filters, the reflector stack provides a different PM at
each wavelength, $\mathbf{R}_{n},~{}n=1..4$. This is seen by considering the
different sum of wavefronts that occurs due to reflection and absorption from
different layers in Figure 4a. Sweeping a tunable laser effectively switches
between reflectors without any actuation at the distal facet. The solid nature
of these structures means they can be characterised before use, for example
using the process in Appendix A.3.1, and are expected to remain stable
throughout operation. The linear nature of this reflector stack design ensures
that each PM, $\mathbf{R}_{n},~{}n=1..4$, is symmetric due to optical
reciprocity Gu  _et al._ (2015).

###  III.2 Selection of wavelengths

For a given reflector stack the optimal filter centre wavelengths, RM
characterisation wavelengths and imaging wavelength must be determined. We
assume here that all filters are long-pass but a system using short-pass
filters would also be feasible. First, we approximate the long-pass filter
transmission function with a sigmoid:

$f(\lambda,\lambda_{filt})=\frac{\tau_{max}}{1+e^{-\frac{\alpha}{\tau_{max}}\left(\lambda-\lambda_{filt}\right)}}$
(32)

where $\tau_{max}$ is the maximum power transmission (fixed here as 0.8 based
on real component data), $\alpha$ determines the steepness of the filter
(measured to be $\sim$0.73 based on numerous commercially available filters,
e.g. ThorLabs FELH0850, which represents a $\sim$5nm wavelength change for 10%
to 90% normalised transmission), and $\lambda_{filt}$ is the centre wavelength
of the filter. For the case of 3 RM characterisation wavelengths
($\lambda_{1}<\lambda_{2}<\lambda_{3}$), this leads to the following matrix:

$\mathbf{F}=\left(\begin{array}[]{c c
c}f(\lambda_{1},\lambda_{a})&f(\lambda_{1},\lambda_{b})&f(\lambda_{1},\lambda_{c})\\\
f(\lambda_{2},\lambda_{a})&f(\lambda_{2},\lambda_{b})&f(\lambda_{2},\lambda_{c})\\\
f(\lambda_{3},\lambda_{a})&f(\lambda_{3},\lambda_{b})&f(\lambda_{3},\lambda_{c})\\\
\end{array}\right)$ (33)

where $\lambda_{a}<\lambda_{b}<\lambda_{c}$ are the centre wavelengths of the
3 filters with transmission curves described by Equation 32. To maximise the
difference between the filters’ behaviour at the interrogation wavelengths, we
minimise the condition number of $\mathbf{F}$, which gives our objective
function, $g$:

$g(\lambda_{1},\lambda_{2},\lambda_{3},\lambda_{a},\lambda_{b},\lambda_{c})=\kappa(\mathbf{F})=\lVert\mathbf{F}^{-1}\rVert_{2}\lVert\mathbf{F}\rVert_{2}$
(34)

where $\kappa(..)$ represents condition number and $\lVert..\rVert_{2}$ is the
$\ell_{2}$ norm. We also want to ensure that at the imaging wavelength
$\lambda_{4}$ ($>\lambda_{3}$) the reflector stack transmits significantly
more light than at the longest characterisation wavelength ($\lambda_{3}$).
This is essential for efficient illumination and imaging. To achieve this, we
define a ratio, $\beta$:

$\beta=\frac{f(\lambda_{4},\lambda_{c})}{f(\lambda_{3},\lambda_{c})}$ (35)

The full optimisation problem is then:

$\begin{array}[]{c}\min\limits_{\lambda_{1},\lambda_{2},\lambda_{3},\lambda_{a},\lambda_{b},\lambda_{c}}g(\lambda_{1},\lambda_{2},\lambda_{3},\lambda_{a},\lambda_{b},\lambda_{c})\\\
\mathrm{s.t.~{}~{}~{}~{}}\lambda_{min}<\lambda_{1}<\lambda_{2}<\lambda_{3}<\lambda_{4}<\lambda_{max}\\\
\lambda_{2}-\lambda_{1},\lambda_{3}-\lambda_{2},\lambda_{4}-\lambda_{3}>\Delta\lambda_{las},~{}~{}~{}~{}\lambda_{4}-\lambda_{1}<\Delta\lambda_{fib}\\\
\lambda_{min}<\lambda_{a}<\lambda_{b}<\lambda_{c}<\lambda_{max},~{}~{}~{}~{}\beta>\beta_{min}\end{array}$
(36)

where $\lambda_{min}$ and $\lambda_{max}$ are the minimum and maximum range of
the tunable laser, $\Delta\lambda_{las}$ is the minimum tuning step of the
laser, and $\Delta\lambda_{fib}$ is the spectral bandwidth of the fibre.

Here, we set $\lambda_{min}$ and $\lambda_{max}$ to near-infrared wavelengths
of 840nm and 860nm respectively Gordon  _et al._ (2018). We set
$\Delta\lambda_{las}$ to 0.2nm to reflect the tuning step of commercial
tunable lasers. $\Delta\lambda_{fib}$ is set to be 7nm representing a
reasonable spectral bandwidth for a fibre of length $\sim$2m Carpenter  _et
al._ (2015); Thompson (2012). Finally, we choose $\beta_{min}=10$ as providing
acceptable isolation between characterisation and imaging. With the objective
and constraints defined, optimisation is performed using a genetic algorithm,
implemented in MATLAB 2016b.

###  III.3 Simulated fibre matrices

For a first proof-of-principle, fibre TMs of size 32$\times$32 are simulated.
This corresponds to 16 pixels of spatial resolution in two orthogonal
polarisations. The matrices are based on a model of a MCF with large core-to-
core coupling. The 16 spatial pixels are considered to represent light guiding
cores organised in a rectangular $4\times 4$ grid with unit spacing (arbitrary
units) between adjacent elements. The power coupling between a given core at
the input and a given core at the output is modelled as decreasing
exponentially with the squared lateral distance between them, i.e.:

$p_{\alpha,\beta}=e^{\left(-\frac{\left(x_{\alpha}-x_{\beta}\right)^{2}+\left(y_{\alpha}-y_{\beta}\right)^{2}}{\sigma_{fib}^{2}}\right)}$
(37)

where $p_{\alpha,\beta}$ is the power coupled between core $\alpha$ at the
input (with coordinates $(x_{\alpha},y_{\alpha})$) and core $\beta$ at the
output (with coordinates $(x_{\beta},y_{\beta})$). $\sigma_{fib}^{2}$ is a
coupling parameter which we set to be 3. This is a Gaussian power coupling
model, which reflects empirical observations of real MCFs Gordon  _et al._
(2019). The power coupling profile is duplicated for the second polarisation
via a Kronecker product with a $2\times 2$ matrix of ones. The complex phase
of each element is drawn randomly from a uniform distribution,
$\theta_{\alpha,\beta,p,s}\sim U(-\pi,\pi)$ where $p=1,2,~{}~{}s=1,2$ indicate
the input and output polarisations respectively. This gives a combined
expression for each matrix element:

$a_{\alpha,\beta,p,s}=e^{\left(\frac{\left(x_{\alpha}-x_{\beta}\right)^{2}+\left(y_{\alpha}-y_{\beta}\right)^{2}}{\sigma_{fib}^{2}}\right)}\cdot
e^{i\theta_{\alpha,\beta,p,s}}$ (38)

The simulated TM, $\mathbf{A}$, is formed of all elements
$a_{\alpha,\beta,p,s}$ in some ordering such that each pair ($\beta$, $s$)
defines a unique row index, and each pair ($\alpha$, $p$) defines a unique
column index. $\mathbf{A}$ is then decomposed using singular value
decomposition:

$\mathbf{A}=\mathbf{U}_{A}\mathbf{S}_{A}\mathbf{V}_{A}^{H}$ (39)

where $\mathbf{U}_{A}$ and $\mathbf{V}_{A}$ are unitary, $\mathbf{S}_{A}$ is a
diagonal matrix containing singular values of $\mathbf{A}$, and $(..)^{H}$
represents a Hermitian transpose. To ensure $\mathbf{A}$ is non-unitary, new
singular values, $s_{n},~{}n=1..32$, are generated such that
$s_{n}=n/32,~{}n=1..32$. This ensures an invertible non-unitary matrix.

If using the zeroth-order model the simulated matrix is used as the TM at all
wavelengths. If using the first-order model the simulated matrix is used as
the TM at wavelength $\lambda_{1}$, $\mathbf{A}_{1}$, and the TMs at the other
wavelengths are generated from the matrix logarithm of $\mathbf{A}_{1}$
($d\mathbf{A}$ from Equation 11) and then applying Equations 15–17 to obtain
$\mathbf{A}_{2}$, $\mathbf{A}_{3}$, $\mathbf{A}_{4}$.

Although this method of generating matrices simulates a MCF, it could just as
easily represent a MMF through an appropriate change of mode basis, e.g. by
redefining the spatial pixels to be Laguerre–Gauss modes Gataric  _et al._
(2019); Carpenter  _et al._ (2014). Simulated TMs are generated using custom
code written in MATLAB 2016b.

###  III.4 Measured fibre matrices

The proposed methodology is further validated by testing the first-order model
and associated recovery using TMs measured from two real fibres. The first is
a MCF with TM of size $1648\times 1648$ measured at a single wavelength of
850nm in Gordon  _et al._ (2018). This TM is used as $\mathbf{A}_{1}$ and TMs
at other wavelengths ($\mathbf{A}_{2}$, $\mathbf{A}_{3}$, $\mathbf{A}_{4}$)
are generated from this by finding the matrix logarithm of $\mathbf{A}_{1}$
($d\mathbf{A}$ from Equation 11) and then applying Equations 15 – 17 to obtain
$\mathbf{A}_{2}$, $\mathbf{A}_{3}$, $\mathbf{A}_{4}$. These are subsequently
used to generate RMs from which the TMs are recovered. This process of
adjusting TMs and generating RMs is implemented using custom code written in
MATLAB 2016b. This case validates that the recovery method is feasible for
large TMs.

The second is a 2m piece of step-index MMF with $420\times 420$ TMs measured
across a wavelength range 1525 – 1567nm in steps of 0.08nm (taken from
Carpenter (2016)). Because of the step-index nature of the fibre, the higher-
order modes have relatively narrow spectral bandwidth (see Appendix F).
Therefore, to improve bandwidth performance only the 110 lowest-order modes
are used, giving submatrices of size $110\times 110$. Here, we use a
wavelength spacing of $\Delta\lambda=0.9$nm giving TMs measured at
characterisation wavelengths $\lambda_{1}=1525.6$nm, $\lambda_{2}=1526.5$nm
and $\lambda_{4}=1527.4$nm, and imaging wavelength $\lambda_{4}=1528.3$nm.
This wavelength separation is sufficient to enable the use of multiple off-
the-shelf optical filters in the reflector stack.

###  III.5 Simulating reflectors

When recovering simulated TMs of size $2M\times 2M$ ($2M=32$ here, as
described in Section III.2), we simulate reflectors by first creating a
diagonal matrix with integers 1 to $2M$ along the main diagonal in some
permutation. Distinct permutations are used for the main diagonals of the
three separate reflectors. Random permutations of the integers 1 to $2M-1$ are
then used to create the subdiagonals and superdiagonals of these reflector
matrices, mimicking the mode-coupling that would be expected in real
reflectors. The singular value decomposition of this matrix is computed, and
the matrix is then reconstructed from the left- and right-singular vector
matrices but with the singular values replaced by a new set, $s_{m}=m/2M$ with
$m=1..2M$. This ensures that the requirement for distinct eigenvalues (see
Section II.2) is satisfied because of the square-root relationship between
singular values and eigenvalues.

Next, we present a method of simulating physically realistic reflectors like
that shown in Figure 4a. This is achieved by extending the propagation matrix
method used for simulating multilayer optical materials (Bragg reflectors or
stacks) Orfanidis (2016). The reflector is modelled as a series of layers, as
shown in Figure 4a, and the propagation operator used is 2D Fresnel
propagation, applied separately to each polarisation. This is in contrast to
conventional propagation matrix method in which the propagation operator is
the complex exponential propagation operator $e^{-ikz}$ where $z$ is distance
and $k$ is wavenumber. The Fresnel propagation operator enables computation of
output fields $\mathbf{p}_{H}^{\prime}$ and $\mathbf{p}_{V}^{\prime}$ (where
subscripts $H$ and $V$ denote horizontal and vertical polarisations
respectively), resulting from input fields $\mathbf{p}_{H}$ and
$\mathbf{p}_{V}$ propagating a distance $l$ through a layer as:

$\mathbf{p}_{H}^{\prime}=\mathcal{F}^{-1}\left(\mathcal{F}\left(\mathbf{p}_{H}\right)e^{-i\frac{2\pi^{2}}{k}\left(\xi_{x}^{2}+\xi_{y}^{2}\right)l}\right)$

where $\mathcal{F}$ is the discrete Fourier transform, $n$ is the refractive
index of the layer, $\lambda$ is the wavelength, $k=\frac{2\pi n}{\lambda}$,
$\xi_{x}$ and $\xi_{y}$ are coordinates in the Fourier plane, and constant
factors are neglected Goodman (1996). This is repeated for $\mathbf{p}_{V}$,
giving a 2D complex vector at each point.

Using this modified method, the reflector stack is simulated using three types
of layer: pure glass, glass with a wire-grid polariser metasurface on the top,
and optical filters. Here, we simulate absorptive optical filters but
reflective filters would work equally well because their behaviour in the
passband is the same. The filters have transmission curves typical of
commercially available components and the centre wavelengths are determined
using the process in Section III.2.

The stack simulated here comprises a first 1mm thick layer of glass, followed
by the 3 filters in succession, each 3mm thick. The metasurfaces, negligibly
thin for propagation purposes, are placed between the glass and the first
filter, between the first and second filters, and between the second and third
filters. The glass and filters have matched refractive indices of 1.52 and for
simplicity we neglect the imaginary part of refractive index, i.e. distance-
dependent loss. Therefore, the only internal reflections occur at the
metasurfaces, reflective filters (if used), and the final glass-air interface.
Each metasurface is simulated by randomly generating diattenuating partial-
polariser Jones matrices for each sampled point of the field. Each point is
assigned a diattenuation angle drawn from a uniform distribution,
$\theta_{D}\sim U(-\pi,\pi)$, and a diattenuation drawn from a uniform
distribution, $D\sim U(0,1)$. The transmitted light at each pixel is computed
by multiplying the sampled field vector by the relevant Jones matrix:

$\left(\begin{array}[]{c}p_{tH}\\\
p_{tV}\end{array}\right)=\sqrt{\alpha}\left(\begin{array}[]{c
c}\cos\theta_{D}&-\sin\theta_{D}\\\
\sin\theta_{D}&\cos\theta_{D}\end{array}\right)\left(\begin{array}[]{c
c}1&0\\\ 0&\sqrt{\frac{1-D}{1+D}}\end{array}\right)$
$\cdot\left(\begin{array}[]{c c}\cos\theta_{D}&\sin\theta_{D}\\\
-\sin\theta_{D}&\cos\theta_{D}\end{array}\right)\left(\begin{array}[]{c}p_{H}^{\prime}\\\
p_{V}^{\prime}\end{array}\right)$ (40)

The light reflected is then given by:

$\left(\begin{array}[]{c}p_{rH}\\\
p_{rV}\end{array}\right)=\sqrt{\alpha}\left(\begin{array}[]{c
c}\cos\theta_{D}&-\sin\theta_{D}\\\
\sin\theta_{D}&\cos\theta_{D}\end{array}\right)\left(\begin{array}[]{c
c}0&0\\\ 0&1-\sqrt{\frac{1-D}{1+D}}\end{array}\right)$
$\cdot\left(\begin{array}[]{c c}\cos\theta_{D}&\sin\theta_{D}\\\
-\sin\theta_{D}&\cos\theta_{D}\end{array}\right)\left(\begin{array}[]{c}p_{H}^{\prime}\\\
p_{V}^{\prime}\end{array}\right)$ (41)

where $\alpha$ is a parameter representing overall power loss in the
metasurface due to plasmonic losses, set to a typical value of $~{}0.8$ here
Williams  _et al._ (2016).

The next step is to determine the overall PMs at the imaging and
characterisation wavelengths, $\mathbf{R}_{n},n=1...4$, and the TM of the
reflector at the imaging wavelength, $\mathbf{A}_{refl}$. This is achieved by
iteratively propagating different input fields forwards and backwards through
the stack until a steady-state solution is reached using custom code written
in MATLAB. This process is repeated at each of the desired wavelengths. The
coupling between spatial modes introduced by diffractive propagation between
layers and the coupling between polarisations introduced by the metasurfaces
produces PMs that with high probability have distinct eigenvalues, which is
necessary for approximate zeroth-order TM recovery that precedes first-order
TM recovery (Section II.2).

The input fields chosen may be derived from the columns of the fibre TM (the
modes of the fibre), because any field exiting the fibre must be a linear
combination thereof. However, the modes of the fibre may be expressed in some
particular basis, e.g. Laguerre–Gauss modes for MMF, and so must first be
converted to Cartesian basis so that Fresnel propagation can be applied. The
light coupled back into the fibre can only be coupled into a linear
combination of the reverse-propagating (i.e. complex conjugated) fibre modes.
Therefore, the steady-state reflected field is re-expressed in the fibre mode
basis. This means the reflector matrices have the same size as the fibre TM
($\in\mathbb{C}^{2M\times 2M}$), but $\mathbf{A}_{refl}$ will in general be
$\in\mathbb{C}^{2N\times 2M}$. Because the $2N$-dimensional vector output of
$\mathbf{A}_{refl}$ is in a Cartesian pixel basis, it is straightforward to
compute sample illumination as per Section II.4.

Both methods of simulating reflector PMs are implemented in custom written
code in MATLAB 2016b.

###  III.6 Noise model

We now consider the effect of noise on the presented TM characterisation and
imaging system. Noise is inherently present in any realistic implementation of
this system, such as that proposed in Appendix A.1, and so must be examined to
demonstrate feasibility and robustness.

The main source of noise in such measurements (for example data acquired using
the system in Figure 10) is assumed to be electronic noise introduced by the
image sensor. However, there are many processing steps (mostly matrix
operations) that transform these raw images into the RMs used for
reconstruction Gordon  _et al._ (2018). The noise in each image pixel is
assumed to be independent, which is consistent with previous fibre imaging
work Gordon  _et al._ (2018, 2019). By the central limit theorem, the repeated
multiplication and addition of noisy variables via matrix operations will
produce circularly symmetric complex zero-mean Gaussian noise. The relatively
strong optical confinement (i.e. wave-guiding) during RM characterisation will
result in large detected optical power and therefore approximately Gaussian
shot noise. Laser power fluctuation (relative intensity noise) is not
inherently Gaussian, but after time averaging and repeated matrix operations
creates a Gaussian contribution to overall noise Yamamoto (1983).

Therefore, we model the noise in the system by adding random matrices with
independent circularly symmetric complex Gaussian elements to each PM
($\mathbf{R}_{1}$, $\mathbf{R}_{2}$, $\mathbf{R}_{3}$) and each RM
($\mathbf{C}_{1}$, $\mathbf{C}_{2}$, $\mathbf{C}_{3}$). For example, a noisy
version of $\mathbf{R}_{3}$ is given by:

$\mathbf{\hat{R}}_{3}=\mathbf{R}_{3}+\mathbf{N}_{3}$ (42)

where each element of $\mathbf{N}_{3}$ is a complex random variable, $n=a+ib$,
drawn from the distributions:

$a\sim\mathcal{N}(0,\sigma_{n})~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}b\sim\mathcal{N}(0,\sigma_{n})$
(43)

The same distributions, parameterised by $\sigma_{n}$, are used to generate
noise matrices $\mathbf{\hat{R}}_{2}$, $\mathbf{\hat{R}}_{3}$,
$\mathbf{\hat{C}}_{1}$, $\mathbf{\hat{C}}_{2}$ and $\mathbf{\hat{C}}_{3}$. We
vary $\sigma_{n}$ over a wide range, including values reflecting empirically
observed experimental noise (see Figure 6) Gordon  _et al._ (2018).

To model noise in image reconstruction we add circularly symmetric complex
Gaussian noise to $\mathbf{A}_{refl}$, $\mathbf{R}_{4}$ and the measured field
at the proximal fibre facet, $\mathbf{y}_{total}$, to get noisy version
$\mathbf{\hat{A}}_{refl}$, $\mathbf{\hat{R}}_{4}$ and
$\mathbf{\hat{y}}_{total}$ respectively. Combined with the noisy estimate of
the TM at the imaging wavelength, $\mathbf{\hat{A}}_{4}$, we can then get a
noisy estimate of the illuminated target, $\mathbf{\hat{r}}_{target}$. A noisy
estimate of the illumination can then be generated using
$\mathbf{\hat{A}}_{4}$ and $\mathbf{\hat{A}}_{refl}$ using Equation 25, and we
can recover a noisy estimate of $\mathbf{\hat{r}}_{target}$ using Equation 31.
Simulation of noise is implemented using custom code in MATLAB 2016b.

###  III.7 Performance metrics

To validate the performance of the presented method, we define several error
metrics. The first is the error between the recovered TM at the imaging
wavelength, $\mathbf{\tilde{A}}_{4}$, and the ideal TM, $\mathbf{A}_{4}$,
given as:

$\epsilon_{mat}=\frac{||\mathbf{\tilde{A}}_{4}-\mathbf{A}_{4}||_{F}}{||\mathbf{A}_{4}||_{F}}=\frac{\left\lVert\left(\mathbf{\tilde{A}}_{4}^{-1}\mathbf{A}_{4}-\mathbf{I}\right)\mathbf{A}_{4}\right\rVert_{F}}{||\mathbf{A}_{4}||_{F}}$
(44)

where $||..||_{F}$ represents the Frobenius norm. Dividing by the norm of the
ideal matrix, $\mathbf{A}_{4}$ ensures the metric presents normalised (or
proportional) error. We then define a second error metric that represents the
error introduced when using a recovered matrix to reconstruct an image:

$\epsilon_{image}=\frac{||\mathbf{\tilde{A}}_{4}^{-1}\mathbf{A}_{4}\mathbf{r}_{target}-\mathbf{r}_{target}||_{2}}{||\mathbf{r}_{target}||_{2}}=\frac{\left\lVert\left(\mathbf{\tilde{A}}_{4}^{-1}\mathbf{A}_{4}-\mathbf{I}\right)\mathbf{r}_{target}\right\rVert_{2}}{||\mathbf{r}_{target}||_{2}}$
(45)

where $\mathbf{r}_{target}$ is the ideal image. A further metric measures
error due to noise in image recovery:

$\epsilon_{image,n}=\frac{||\mathbf{\tilde{r}}_{target}-\mathbf{r}_{target}||_{2}}{||\mathbf{r}_{target}||_{2}}$
(46)

where $\mathbf{\tilde{r}}_{target}$ is the estimated recovered image as
defined in Equation 30.

###  III.8 Recovering TMs

Having produced the necessary instantaneous RMs ($\mathbf{C}_{n},~{}n=1..3$),
based on real or simulated TMs, the next step is to solve for the TMs using _a
priori_ known reflector PMs, $\mathbf{R}_{n}$ with $n=1..3$. For TMs and RMs
generated using the zeroth-order model, the equations of Section II.2 are
solved using code written in MATLAB 2016b, including solving the systems of
linear equations generated by Equation 60.

In the first-order case, the recovery algorithm depicted in Figure 2 is
implemented in the TensorFlow package via Python because of its efficient GPU
utilisation and hence superior speed. Initially, we use the Adam stochastic
gradient descent optimiser Kingma and Ba (2014), changing to a conventional
gradient descent optimiser for fine adjustment in the last few iterations.

The zeroth-order recovery algorithm typically finishes in less than 1s.
However, first-order iterative algorithm is significantly slower and can take
several hours for 1500 iterations on a $1648\times 1648$ matrix on a single
Tesla K40 GPU. For smaller $32\times 32$ matrices requiring about 500
iterations, the running time is typically $<5$ seconds.

##  IV Results

###  IV.1 Simulated TMs

We first randomly generate a single $32\times 32$ non-unitary matrix,
$\mathbf{A}_{1}$, and then $\mathbf{A}_{2}$, $\mathbf{A}_{3}$,
$\mathbf{A}_{4}$ using the first-order model. We then randomly generate
reflector PMs $\mathbf{R}_{n}$ with $n=1...4$ and use $\mathbf{A}_{1}$ to
simulate RMs $\mathbf{C}_{n}$ with $n=1...3$ and perform TM recovery using the
first-order method (Figure 5a) with element-wise error $<0.05$ (Figure 5b).
The singular values of this matrix vary from 1 down to 0.03 (Figure 5c),
verifying that this recovery can be performed on very non-unitary matrices.

Figure 5: Recovery of a simulated TM using the first-order recovery algorithm:
a) The original simulated TM is compared with the recovered TM (noise level
$\sigma_{n}=10^{-3}$) demonstrating successful recovery. b) Proportional
element-wise error in TM reconstruction (maximum $<0.05$) after 300 iterations
of first-order model. c) Singular values of TM demonstrating non-unitarity.

Next, we examine the impact of noise by repeating this recovery process for 64
different noise realisations with increasing noise power (see Section III.6).
Matrices are generated using either the zeroth- or first-order model and then
reconstructed using the respective recovery algorithm. The resultant error in
image reconstruction using the recovered TM (Figure 6a) for the zeroth-order
model, in which the recovery is largely analytical, deteriorates rapidly in
performance beyond a noise threshold (visualised as a ‘hump’ in Figure 6a). In
the first-order model, in which recovery relies on an optimisation process,
such instability is not observed, thus it is suspected to arise from numerical
error accumulation when computing analytical results. The first-order model
exhibits improved stability at higher noise levels (such as those encountered
experimentally Gordon  _et al._ (2018)) because the optimisation process, by
definition, does not terminate until error (including numerical error) is
minimised. At lower noise levels, first-order performance is limited by the
number of iterations of the algorithm, which is evident from the impact of
increasing the number of iterations from 300 to 500.

Figure 6: Impact of measurement noise and different perturbation conditions on
recovery. a) Error in reconstructed images for different levels of noise in
measured quantities ($\epsilon_{image}$ and $\sigma_{N}$ are defined in
Section III.7). At low noise levels the error of the first-order model is seen
to be limited by the number of iterations, but at higher noise levels
converges to a trend consistent with the zeroth-order model. The first-order
approach is more robust at higher noise levels due to the use of iterative
optimisation. A typical upper bound for noise encountered experimentally
Gordon  _et al._ (2018) is indicated to show that recovery from real
experimental data is feasible. b) Histogram showing effect of 200 different
random TM realisations, i.e. different fibre perturbations, under the first-
order model. Error after the initial zeroth-order approximation (see Figure 2)
can be large ($>1$) but after 500 iterations, 50% of TMs have errors $<0.1$
and 95% have errors $<0.4$.

Furthermore, testing the effect of 200 different random TM realisations that
could result from perturbations of the fibre or manufacturing variations
(Figure 6b) shows that after 500 iterations using the first-order recovery
algorithm, the error converges towards zero with $>$50% of matrices having
error $<0.1$ and 90% having error $<0.4$. Negligibly small noise power is used
so as to only observe the effect of perturbations.

###  IV.2 Real MCF TMs

Having established the potential of the first-order model using simulated
matrices, we then applied it to an experimentally measured MCF TM Gordon  _et
al._ (2018). Using the optimisation procedure of III.2, we determine the
wavelengths to be used as follows: filters centred at 840.0nm, 842.7nm and
849.3nm, characterisation wavelengths of $\lambda_{1}$ = 843.4nm,
$\lambda_{2}$ = 844.9nm and $\lambda_{3}$ = 846.5nm, and imaging wavelength
$\lambda_{4}$ = 850.4nm.

The non-square TM is first ‘downsampled’ to a square form of size $1648\times
1648$ (Figure 7a). This downsampling process can be achieved via many methods,
for example that detailed in Appendix A.2, and can be reversed at the end of
recovery to obtain an estimate for the original non-square TM.

The measured TM is used as $\mathbf{A}_{1}$ and we apply the first-order model
to simulate $\mathbf{A}_{2}$, $\mathbf{A}_{3}$ and $\mathbf{A}_{4}$. Next,
realistic reflector PMs, $\mathbf{R}_{n}$ with $n=1...4$, are simulated (see
Section III.5) and RMs $\mathbf{C}_{1}$, $\mathbf{C}_{2}$ and $\mathbf{C}_{3}$
computed using Equations 5–7. Using the first-order algorithm we then recover
the TM at the imaging wavelength, $\mathbf{A}_{4}$ (Figure 7a). After 1500
iterations, the maximum element-wise error of the recovered matrix is
$<10^{-3}$ (Figure 7b). Because of the dense core spacing ($<4\mu$m) there is
significant coupling between cores at the wavelength used, 852nm Gordon  _et
al._ (2018), meaning the MCF begins to exhibit mode-coupling properties more
like MMF. The singular values (Figure 7d) show that the matrix is non-unitary
but has a low condition number ($\sim$4), again similar to experimental
measurements of MMF Carpenter  _et al._ (2014).

Figure 7: First-order recovery of an experimentally-measured MCF TM. a)
Original TM using data from Gordon  _et al._ (2018) and the first-order
recovered TM, showing high visual similarity indicating successful recovery.
b) Proportional element-wise error in TM reconstruction (maximum $<10^{-3}$).
c) Singular values of TM demonstrating non-unitarity.

Having recovered the TM, we next test the imaging procedure using random
illumination. A $28\times 28$ image is simulated and then successfully
recovered (Figure 8). The image, $\mathbf{r}_{target}$ with reference to
Figure 1, comprises 28 $x$ positions $\times$ 28 $y$ positions $\times$ 2
polarisations $=$ 1568 degrees of freedom (i.e. $2N=1568$ as defined in
Section II.4), which ensures that full reconstruction through a TM of
dimension 1648$\times$1648 (i.e. $2M=1648$ as defined in Section II.4) is
possible. The total image reconstruction error, $\epsilon_{n,image}$ (Equation
46), is $<0.1$ and arises primarily from illumination correction. Lower
illumination levels result in higher noise (see Section II.4) and could thus
be improved by averaging over multiple random illumination conditions i.e.
speckle averaging Choi  _et al._ (2012). Under conditions of perfectly uniform
illumination, we calculate that $\epsilon_{n,image}$ becomes $<0.002$, giving
a lower bound on error dominated by other facts e.g. computational
limitations.

Figure 8: Successful reconstruction of an amplitude, phase and polarisation
image of a target at the distal facet using a TM recovered using the first-
order model. The recovered image noise is largely driven by illumination
correction (see Section II.4).

###  IV.3 Real multi-wavelength MMF TMs

Using a multi-wavelength fibre TM measured from a 2m piece of step-index MMF
Carpenter (2016) to generate RMs, $\mathbf{C}_{n},~{}n=1..3$, also led to
successful TM, $\mathbf{A}_{4}$, recovery (Figure 9a). In this case, the
characterisation wavelengths chosen are $\lambda_{1}=1525.6$nm,
$\lambda_{2}=1526.5$nm and $\lambda_{3}=1527.4$nm, and the imaging wavelength
is $\lambda_{4}=1528.3$nm. A $110\times 110$ subset of the full TM data set is
used in order to limit spectral bandwidth appropriately (see Appendix F).
After 475 iterations taking 14.5 minutes, the algorithm converges with
significant visual similarity (Figure 9a). The matrix error, $\epsilon_{mat}$
(Equation 44) reaches a value of 0.26. Looking at the element-wise error
(Figure 9b) we see that error is typically $<0.05$ except in the 3 lowest-
order modes, in which the value reaches 0.2. This error is observed to
manifest as phase errors, with the amplitude values being largely accurate.
Importantly, the singular value distribution agrees with that expected of MMF
(Figure 9c).

Figure 9: First-order recovery of an experimentally-measured step-index MMF
TM, using TMs recorded at multiple wavelengths (1525.6nm, 1526.5nm, 1527.4nm
and 1528.3nm) Carpenter (2016). a) Original TM using data from Carpenter
(2016) and the first-order recovered TM, showing high visual similarity
indicating successful recovery. b) Proportional element-wise error in TM
reconstruction values typically $<$0.05 except for the very lowest order modes
that exhibit error $\sim$0.2 due to phase error. c) Singular values of
measured TM demonstrating non-unitarity.

##  V Discussion

The ability to form images through hair-thin optical fibres is currently
limited due to the need for _in situ_ TM calibration. To overcome this
limitation, we propose a method for recovering the instantaneous TM of a fibre
based only on instantaneous reflection-mode measurements and _a priori_
characterisation reflection- and transmission-mode measurements. Calibration
is achieved using three reflectors placed at the distal tip of the fibre, each
of which produces a different reflection matrix based on a small modulation in
wavelength, from which the imaging TM can be determined. As the fibre TM will
differ slightly at each of these wavelengths, we developed a first-order
method of modelling and compensating this change using matrix exponentials,
which is valid within the spectral bandwidth of the fibre. Under these
assumptions, we demonstrate that it is possible to recover TMs and perform
imaging using realistic optical components with acceptable noise tolerance for
both simulated and measured TMs, including an experimentally measured MCF TM
and a multi-wavelength MMF TM.

Nonetheless, several challenges must be overcome for experimental deployment
of the proposed method. Given the demonstrated validity of the current model
within the spectral bandwidth of the fibre, at 850nm the wavelength modulation
would need to be within $\sim 5$nm for a typical 1m length of MCF Thompson
(2012) or at 1550nm within $\sim 10$nm for a typical 2m length of MMF
Carpenter  _et al._ (2015). In biomedical endoscopy, fibre lengths are
typically 1–2m, however, for applications in industrial inspection, greater
lengths may be needed. To extend the model to work for smaller bandwidth
fibres, e.g. longer fibres or MMFs with very large numbers of modes, either
the characterisation bandwidth could be reduced or a more complex (‘high
order’) fibre model used. The characterisation bandwidth is typically limited
by the ‘sharpness’ of available optical filters, which enables significant
modulations in reflectance/absorbance over very small wavelength ranges.
Custom-fabricated filters may offer improved performance over off-the-shelf
products. Alternatively, more complex models of the fibre could be used.
Complete propagation models of graded-index MMF have been used to accurately
compensate bending Plöschner  _et al._ (2015) but require precise _a priori_
knowledge of the refractive index profile. A major advantage of the approach
presented here is that no such prior knowledge of the fibre is required,
enabling the use of more complex refractive index profiles such as MCF. It may
be possible to develop a more complex model that uses the differential changes
in the TM with respect to wavelength to model the fibre over a larger
bandwidth for example using machine learning techniques Rahmani  _et al._
(2018); Yang  _et al._ (2018). As they do not require a reference model,
machine learning techniques may also enable the use of a smaller reflector
stack, with only 2 reflectors, though our preliminary investigations have not
conclusively shown that this produces a unique solution when using an
optimisation process to recover the TM.

Another challenge compared to existing transmission-mode fibre
characterisation systems is the introduction of additional characterisation
steps, both before and during imaging for instantaneous TM recovery (see
Appendix A.3). For example, to recover a single TM will require 3 times as
many measurements because of the 3 distinct wavelengths. Speed can be improved
by experimental design, for example using high-speed digital micromirror
devices instead of liquid-crystal SLMs Mitchell  _et al._ (2016), or by
exploiting _a priori_ known sparsity structures of fibres to parallelise
spatial measurements Gordon  _et al._ (2018). Further, because the
characterisation measurements use laser light at distinct wavelengths, they
could be captured simultaneously by placing appropriate filters in front of 3
separate image sensors, one for each wavelength. Computational recovery of TMs
and images also presents a speed limitation that could be addressed by using
additional GPUs or specialised hardware like an FPGA or ASIC, fine-tuning the
algorithm and solver settings, or only reconstructing a single polarisation (a
4-fold speed increase). If an analytical solution for the first-order case
could be found, a dramatic speed increase would result.

A final challenge is fabricating and installing an appropriate reflector stack
at the distal tip of the fibre. While optical filters can be purchased off-the
shelf, wire-grid polariser metasurfaces most likely require custom fabrication
using nanofabrication techniques such as electron-beam lithography. Single-
layer devices of this nature have already been demonstrated Williams  _et al._
(2016, 2017), but extending these to make reflector stacks will involve a
number of additional deposition, lithography and processing steps. Though the
metasurfaces must have specific polarisation transmission and reflection
properties, they need not be fabricated to a specific design with high
precision because they can be characterised _in situ_.

Because the method does not place any requirements on the singular value
distribution of the TM other than that the TM is invertible (i.e. unitarity is
not required), it may also be applicable to more general scattering media. We
provide further validation of this by demonstrating recovery of highly non-
unitary simulated TMs with a very uneven distribution of singular values,
similar to the quarter-circle distribution found in random scattering media
Popoff  _et al._ (2011). Nonetheless, it would require a precisely
characterised reflector to be placed inside the medium, which would require an
invasive step in tissue imaging applications.

##  VI Conclusion

In summary, we have developed a new approach to determine the transmission
matrix of a multi-mode optical fibre without requiring access to the distal
facet. We have demonstrated successful transmission matrix recovery and
imaging using realistic optical components with acceptable noise tolerance for
simulated transmission matrices, and those experimentally measured from a
multi-core fibre and a multi-wavelength multi-mode fibre. The proposed method
paves the way for experimental realisation of lensless imaging through hair-
thin optical fibres.

###### Acknowledgements.

GSDG acknowledges funding from Cancer Research UK (C47594/A21102,
C55962/A24669); and a pump-priming award from the Cancer Research UK Cambridge
Centre Early Detection Programme (A20976). CW acknowledges funding from a
Cancer Research UK Pioneer Award (C55962/A24669). MG acknowledges funding from
Engineering and Physical Sciences Research Council for the Centre for
Mathematical and Statistical Analysis of Multimodal Clinical Imaging
(EP/N014588/1). RPM acknowledges funding from the Engineering and Physical
Sciences Research Council (EP/L015889/1). SEB acknowledges funding from Cancer
Research UK (C47594/A16267, C14303/A17197, C47594/A21102); and the EU FP7
agreement (FP7-PEOPLE-2013-CIG-630729). GSDG would like to thank Peter
Christopher from the University of Cambridge for useful discussions.

##  Appendix A Conceptual experimental implementation

###  A.1 Experimental set-up

A conceptual experimental set-up for reflection-mode TM recovery is shown in
Figure 10. There are two key sections: the illumination arm and the detection
arm. These perform similar functions to previous work Gordon  _et al._ (2018)
but are both located at the proximal fibre facet. There is also an additional
input that allows characterisation in transmission mode prior to use. Though
our previous fibre TM characterisation set-ups have used non-interferometric
phase retrieval Gordon  _et al._ (2018), an optional reference arm for
interferometric imaging is included as it is widely used method in fibre
characterisation Cižmár and Dholakia (2012). The tunable laser is required to
provide a high-coherence source at the different characterisation and imaging
wavelengths. These can typically be tuned in steps of 0.1–0.2nm, sufficient to
cover a range of fibre TMs with different spectral bandwidths.

Figure 10: Conceptual experimental set-up for fibre TM characterisation in
reflection mode configuration using distal reflector stack. SLM: spatial light
modulator, PBS: polarising beam splitter, NPBS: non-polarising beam splitter,
$\lambda/2$: half waveplate.

Proposed designs of a reflector stack are shown in Figure 4. The plasmonic
metasurface reflectors can be fabricated on a glass substrate using electron
beam lithography patterning followed by metal deposition to create high-
resolution wire-grid polarisers Williams  _et al._ (2016). Fabrication could
be scaled-up for commercial applications using, for example, UV lithography or
interference lithography Neisser and Wurm (2015). The absorptive or reflective
long-pass filters used for the stack can be purchased off-the-shelf or made to
order with sufficiently small gradations in centre wavelength. Stability of
the stack over time is important to ensure the reflector matrices remain
relatively constant. This can be achieved by using stable metals such as gold
and depositing protective layers, e.g. MgF, on top of metasurfaces.

###  A.2 Downsampling Reflection Matrices

In practical implementations we may measure RMs that are non-square, denoted
$\tilde{\mathbf{C}}_{n}$ with $n=1...4$, due to different sampling schemes for
$\mathbf{x}$ and $\mathbf{y}$ of Figure 1. For example, the number of pixels
on the camera may greatly exceed the number of different input calibration
fields used to characterise the RM – in our previous work characterising TMs
there were $\sim 10^{3}$ times as many pixels as calibration inputs Gordon
_et al._ (2018).

In this scenario, consider $2M$ calibration samples
$(\tilde{\mathbf{x}}_{s},\tilde{\mathbf{y}}_{s})_{s=1}^{2M}$ where
$\tilde{\mathbf{x}}_{s}\in\mathbb{C}^{P\times 1}$ and
$\tilde{\mathbf{y}}_{s}\in\mathbb{C}^{Q\times 1}$. Denoting
$\tilde{\mathbf{C}}_{n}\in\mathbb{C}^{Q\times P}$ as a mapping of the fibre
$\tilde{\mathbf{y}}=\tilde{\mathbf{C}}_{n}\tilde{\mathbf{x}}$ (with reference
to Figure 1), we wish to determine a downsampling process to obtain a square
RM, $\mathbf{C}\in\mathbb{C}^{2M\times 2M}$ that may be used in our recovery
algorithms. One way of doing this is to estimate the largest singular values
and corresponding left singular vectors of
$\tilde{\mathbf{C}}_{n}\in\mathbb{C}^{Q\times P}$, which should match those of
$\mathbf{C}_{n}\in\mathbb{C}^{2M\times 2M}$. This is done in the following
way. Let

$\tilde{\mathbf{X}}_{\mathrm{cal}}:=[\tilde{\mathbf{x}}_{1},\ldots,\tilde{\mathbf{x}}_{2M}]\in\mathbb{C}^{P\times
2M},\quad\tilde{\mathbf{Y}}_{\mathrm{cal}}:=[\tilde{\mathbf{y}}_{1},\ldots,\tilde{\mathbf{y}}_{2M}]\in\mathbb{C}^{Q\times
2M}$ (47)

so that
$\tilde{\mathbf{Y}}_{\mathrm{cal}}=\tilde{\mathbf{C}}\tilde{\mathbf{X}}_{\mathrm{cal}}$.
If the coefficients of an input signal
$\tilde{\mathbf{x}}\in\mathbb{C}^{P\times 1}$ with respect to the
$2M$-dimensional calibration basis are denoted by
$\tilde{\mathbf{x}}_{\mathrm{cal}}\in\mathbb{C}^{2M\times 1}$, i.e.,
$\tilde{\mathbf{x}}:=\tilde{\mathbf{X}}_{\mathrm{cal}}\tilde{\mathbf{x}}_{\mathrm{cal}}$,
then its output signal is given as
$\tilde{\mathbf{y}}=\tilde{\mathbf{C}}\tilde{\mathbf{x}}=\tilde{\mathbf{C}}\tilde{\mathbf{X}}_{\mathrm{cal}}\tilde{\mathbf{x}}_{\mathrm{cal}}=\tilde{\mathbf{Y}}_{\mathrm{cal}}\tilde{\mathbf{x}}_{\mathrm{cal}}$.
In other words, $\tilde{\mathbf{Y}}_{\mathrm{cal}}$ is a mapping of the fibre
from the $2M$-dimensional calibration basis to $\mathbb{C}^{Q\times 1}$, and
it is hoped that the largest $2M$ singular values and corresponding left
singular vectors of $\tilde{\mathbf{Y}}_{\mathrm{cal}}$ approximate those of
$\tilde{\mathbf{C}}_{n}$ and $\mathbf{C}_{n}$. We therefore first truncate the
singular values and vectors of $\tilde{\mathbf{C}}_{n}$ to produce
$\tilde{\tilde{\mathbf{C}}}_{n}\in\mathbb{C}^{Q\times 2M}$. Next, we further
compress the matrix $\tilde{\tilde{\mathbf{C}}}_{n}$ into a square matrix
$\mathbf{C}_{n}\in\mathbb{C}^{2M\times 2M}$ which preserves the same singular
values and left singular vectors. This can be done via a process such as that
described in Gordon  _et al._ (2019). Following this process, $\mathbf{C}_{n}$
is square and so can be used for TM recovery as described in Sections II.2 and
II.3.

This enables a factorisation of $\tilde{\mathbf{C}}_{n}$ into the product of
$\mathbf{C}_{n}\in\mathbb{C}^{2M\times 2M}$ and ‘downsampling matrix’,
$\mathbf{T}_{n}\in\mathbb{C}^{P\times 2M}$ such that
$\tilde{\mathbf{C}}_{n}=\mathbf{T}_{n}\mathbf{C}_{n}$. Once the TM has been
determined, $\mathbf{T}_{n}$ can be used to reconstruct
$\tilde{\mathbf{C}}_{n}$ if desired. The process of downsampling may then be
repeated for the other wavelengths, $n=2..4$.

###  A.3 Operating process

Figure 11 shows a flow chart for the full operating process, using the set-up
of Figure 10, for recovering TMs and performing imaging.

Figure 11: Flowchart detailing the full process of TM recovery and imaging
using the reflection mode endoscope set-up. Details of specific sub-tasks are
given in figures and sections as indicated.

####  A.3.1 Characterisation before first use

Before using the system for the first time, we must take accurate measurements
of the reflector stack RMs at all wavelengths and its TM at the imaging
wavelength (Figure 12). These are considered to be stable over long periods so
this is considered a ‘one-off’ process. A flowchart detailing this process is
given in Figure 12. First, the fibre-only TMs ($\mathbf{A}_{n},~{}n=1..4$) are
characterised using the transmission mode input. Then, the reflector stack is
carefully fixed to the distal fibre facet without perturbing the fibre.
Finally, the characterisation process is repeated to measure the RMs of the
fibre with the reflector attached ($\mathbf{C}_{n},~{}n=1..4$). Using the
fibre-only TMs, the PMs ($\mathbf{R}_{n}$ with $n=1..4$) can be determined
from the measured RMs. The TM of the reflector stack at the imaging wavelength
($\mathbf{A}_{s}tack$) is also determined at this stage.

Figure 12: Flowchart detailing the series of measurements required to
accurately characterise the reflector matrices of the reflector stack at
different wavelengths prior to first use of the system. This is achieved using
the experimental set-up of Figure 10.

Non-square measured TMs $\tilde{\mathbf{A}}_{n}$ and RMs
$\tilde{\mathbf{C}}_{n}$ with $n=1..4$ may, if needed, be downsampled to form
square approximations as described in Appendix A.2.

####  A.3.2 Recovering instantaneous transmission matrix

The sequence of experimental measurements required to record instantaneous RMs
at each respective wavelength ($\mathbf{C}_{1}$, $\mathbf{C}_{2}$ and
$\mathbf{C}_{3}$ of Equation 5 – Equation 7) is illustrated in Figure 13.
Fundamentally, the process involves sending known calibration patterns into
the proximal facet using SLM1 of Figure 10, and recording the resultant field
after a round-trip (down the fibre, off the distal reflector, back through the
fibre) on the camera in the detection arm of Figure 10.

Figure 13: Flowchart detailing the series of measurements required to recover
instantaneous RMs using the set-up of Figure 10.

Although the RMs at different wavelengths are here measured sequentially, to
speed up operation it may be possible to measure them in parallel, owing to
linear nature of all optical elements involved. This could, for example, be
achieved by using multiple laser sources multiplexed together and using
several image sensors, each with its own optical filter, such that light from
any given laser source (having a specific wavelength) is detected on only one
sensor.

##  Appendix B Derivation of zeroth-order model

Using the zeroth-order assumption implied by Equation 8, we can rewrite
Equations 5 to 7 as:

$\mathbf{C}_{1}=\mathbf{A}^{\top}\mathbf{R}_{1}\mathbf{A}$ (48)

$\mathbf{C}_{2}=\mathbf{A}^{\top}\mathbf{R}_{2}\mathbf{A}$ (49)

$\mathbf{C}_{3}=\mathbf{A}^{\top}\mathbf{R}_{3}\mathbf{A}$ (50)

Starting from Equation 48 and Equation 49 we derive:

$\mathbf{C}_{2}^{-1}\mathbf{C}_{1}=\mathbf{A}^{-1}\mathbf{R}_{2}^{-1}\mathbf{R}_{1}\mathbf{A}$

which we compactly write as:

$\mathbf{C}_{\alpha}=\mathbf{A}^{-1}\mathbf{R}_{\alpha}\mathbf{A}$ (51)

where
$\mathbf{C_{\alpha}}=\mathbf{C}_{2}^{-1}\mathbf{C}_{1}\in\mathbb{C}^{2M\times
2M}$,
$\mathbf{R}_{\alpha}=\mathbf{R}_{2}^{-1}\mathbf{R}_{1}\in\mathbb{C}^{2M\times
2M}$. Equation 51 can then be transformed into:

$\mathbf{R}_{\alpha}\mathbf{A}-\mathbf{A}\mathbf{C}_{\alpha}=\mathbf{0}$ (52)

In order to avoid the trivial solution $\mathbf{A}-\mathbf{A}=\mathbf{0}$ to
Equation 52, $\mathbf{R}_{\alpha}$ and $\mathbf{C}_{\alpha}$ must not be equal
to the identity matrix, which is ensured by choosing different reflectors
$\mathbf{R}_{1}$ and $\mathbf{R}_{2}$. It is important to note from Equation
51 that $\mathbf{C}_{\alpha}$ and $\mathbf{R}_{\alpha}$ are similar matrices
and so have the same eigenvalues. In a similar fashion we can derive from
Equations 49 and 50:

$\mathbf{R}_{\beta}\mathbf{A}-\mathbf{A}\mathbf{C}_{\beta}=\mathbf{0}$ (53)

where $\mathbf{C_{\beta}}=\mathbf{C}_{3}^{-1}\mathbf{C}_{2}$ and
$\mathbf{R}_{\beta}=\mathbf{R}_{3}^{-1}\mathbf{R}_{2}$. Finally, from
Equations 48 and 50:

$\mathbf{R}_{\gamma}\mathbf{A}-\mathbf{A}\mathbf{C}_{\gamma}=\mathbf{0}$ (54)

where $\mathbf{C_{\gamma}}=\mathbf{C}_{3}^{-1}\mathbf{C}_{1}$ and
$\mathbf{R}_{\gamma}=\mathbf{R}_{3}^{-1}\mathbf{R}_{1}$. Equations 52–54 are
examples of Sylvester equations so we can exploit existing methods to solve
them. Because Equations 52–54 have RHS $=\mathbf{0}$, the trivial solution
$\mathbf{A}=\mathbf{0}$ satisifies all three. However, we know that there must
exist at least one non-trivial solution for $\mathbf{A}$ because if
$\mathbf{A}=\mathbf{0}$ then
$\mathbf{C}_{1},\mathbf{C}_{2},\mathbf{C}_{3}=\mathbf{0}$ by Equations 5–7.
This would imply a system with 100% power loss which we know physically is not
the case. Therefore, we aim to determine the space of non-trivial solutions
(i.e. the nullspace of Equation 52, 53 or 54) by implementing the
Bartels–Stewart algorithm for solving Sylvester equations of the form
$\mathbf{D}\mathbf{X}-\mathbf{X}\mathbf{E}=\mathbf{F}$ for the special case
where $\mathbf{F}=\mathbf{0}$ Bartels and Stewart (1972).

We begin by using the Schur matrix decomposition Horn and Johnson (2013). This
states that any square matrix, $\mathbf{H}$ can be decomposed as:

$\mathbf{H}=\mathbf{Q}\mathbf{U}\mathbf{Q}^{-1}$ (55)

where $\mathbf{Q}$ is unitary and $\mathbf{U}$ is an upper triangular matrix.
Importantly, the diagonal of $\mathbf{U}$ comprises the eigenvalues of
$\mathbf{H}$ sorted in descending order of magnitude. With reference to
Equation 52 we can then write:

$\mathbf{R}_{\alpha}=\mathbf{Q}_{R_{\alpha}}\mathbf{U}_{R_{\alpha}}\mathbf{Q}_{R_{\alpha}}^{-1}$
(56)

$\mathbf{C}_{\alpha}=\mathbf{Q}_{C_{\alpha}}\mathbf{U}_{C_{\alpha}}\mathbf{Q}_{C_{\alpha}}^{-1}$
(57)

Substituting Equation 56 and Equation 57 back into Equation 52 gives:

$\mathbf{Q}_{R_{\alpha}}\mathbf{U}_{R_{\alpha}}\mathbf{Q}_{R_{\alpha}}^{-1}\mathbf{A}-\mathbf{A}\mathbf{Q}_{C_{\alpha}}\mathbf{U}_{C_{\alpha}}\mathbf{Q}_{C_{\alpha}}^{-1}=\mathbf{0}$

$\mathbf{U}_{R_{\alpha}}\mathbf{Q}_{R_{\alpha}}^{-1}\mathbf{A}\mathbf{Q}_{C_{\alpha}}-\mathbf{Q}_{R_{\alpha}}^{-1}\mathbf{A}\mathbf{Q}_{C_{\alpha}}\mathbf{U}_{C_{\alpha}}=\mathbf{0}$

$\mathbf{U}_{R_{\alpha}}\mathbf{A}^{\prime}-\mathbf{A}^{\prime}\mathbf{U}_{C_{\alpha}}=\mathbf{0}$
(58)

where

$\mathbf{A}^{\prime}=\mathbf{Q}_{R_{\alpha}}^{-1}\mathbf{A}\mathbf{Q}_{C_{\alpha}}$
(59)

and $\mathbf{A}^{\prime}\in\mathbb{C}^{2M\times 2M}$. Because of the
degenerate nature of Equation 58 (i.e. RHS $=\mathbf{0}$) there are multiple
solutions to this equation – the trivial solution plus at least one non-
trivial. We reduce the space of possible solutions by requiring that the
eigenvalues of $\mathbf{R}_{\alpha}$ are distinct, with the result that only
the diagonal elements of $\mathbf{A}^{\prime}$ are indeterminate (see Appendix
B.1 for full derivation). This is arranged by suitable design of the
reflectors $\mathbf{R}_{1}$ and $\mathbf{R}_{2}$, which is not in general
difficult because even randomly generated matrices produce distinct
eigenvalues with high probability Marčenko and Pastur (1967); Popoff  _et al._
(2010).

Equation 58 is another Sylvester equation but, crucially,
$\mathbf{U}_{R_{\alpha}}$ and $\mathbf{U}_{C_{\alpha}}$ are upper triangular
matrices so we can solve for $\mathbf{A}^{\prime}$ element by element as it is
typically done in the Bartels–Stewart algorithm. It follows that elements of
$\mathbf{A}^{\prime}$ are related to elements of $\mathbf{U}_{C_{\alpha}}$ and
$\mathbf{U}_{R_{\alpha}}$ as:

$\begin{array}[]{c}\left[ur_{P,Q}a_{Q,Q}-\sum_{s=P}^{Q-1}uc_{s,Q}a_{P,s}\right]+(ur_{P,P}-uc_{Q,Q})a_{P,Q}+\\\
\sum_{q=P+1}^{Q-1}ur_{P,q}a_{q,Q}=0\end{array}$ (60)

where $uc_{P,Q}$ represents an element of $\mathbf{U}_{C_{\alpha}}$ at row $P$
and column $Q$, $ur_{P,Q}$ represents an element of $\mathbf{U}_{R_{\alpha}}$
at row $P$ column $Q$, and $a_{P,Q}$ represents an element of
$\mathbf{A}^{\prime}$ at row $P$ and column $Q$ (see Appendix B.1 for full
derivation).

Varying $P$ from $1$ to $Q-1$ in Equation 60 produces $Q-1$ equations. If we
arbitrarily set the diagonal elements, $a_{m,m}$, of $\mathbf{A^{\prime}}$, we
are left with $Q-1$ unknowns and can solve this system of linear equations for
all elements of the $Q^{th}$ column of $\mathbf{A^{\prime}}$. By varying $Q$
from $2$ to $2M$ and solving each resultant system of equations, we can fully
determine all non-diagonal elements of $\mathbf{A}^{\prime}$ given a set of
arbitrary diagonal elements, $a_{m,m}$ ($m=1..2M$). This repeated solution of
systems of linear equations could be expressed as a single matrix equation,
but the matrix would be of size $M(2M-1)\times M(2M-1)$ creating a memory
requirement $\sim M^{4}$. The equation systems are therefore solved
sequentially for computational tractability.

When an arbitrary set of $2M$ complex numbers (i.e. a vector
$\in\mathbb{C}^{2M}$) is selected for the diagonal elements of
$\mathbf{A}^{\prime}$ in Equation 58 (or equivalently elements $a_{m,m}$ where
$m=1..2M$ of Equation 60) the remaining elements of $\mathbf{A}^{\prime}$ can
be solved. This implies that the ‘true’ solution for $\mathbf{A}^{\prime}$ can
be expressed as a linear combination of at most $2M$ matrices,
$\mathbf{A}_{m}^{\prime}$ where $m=1..2M$, each of which is a solution of
Equation 58. Together, these $2M$ matrices comprise a basis for
$\mathbf{A}^{\prime}$ (equivalently, they comprise a nullspace of Equation
58). Using Equation 59 this basis can be converted to a basis for $\mathbf{A}$
(equivalently, a nullspace for Equation 52), denoted $\mathbf{\tilde{A}}_{m}$
where $m=1..2M$. To improve robustness to noise and numerical error an
iterative Gram–Schmidt orthogonalisation procedure, modified to account for
the conversion from $\mathbf{A^{\prime}}_{m}$ to $\mathbf{\tilde{A}}_{m}$, can
be applied to create an orthogonal set of $\mathbf{\tilde{A}}_{m}$ where
$m=1..2M$ (see Appendix C). The process of determining a basis for
$\mathbf{A}$ is summarised in Figure 14.

Figure 14: Flowchart detailing the process of finding a basis for the matrix
$\mathbf{A}$ using the zeroth-order solution method.

Having determined a suitable basis we express the matrix we wish to recover,
$\mathbf{A}$, as:

$\mathbf{A}=w_{1}\mathbf{\tilde{A}}_{1}+..+w_{2M}\mathbf{\tilde{A}}_{2M}$ (61)

where $w_{m}$ represents the complex-valued weight of the $m^{th}$ basis
element, the matrix $\mathbf{\tilde{A}}_{m}$. The solution space now has only
$2M$ degrees of freedom (represented by $w_{m}$) so computational complexity
can be further reduced by considering a subset of elements of $\mathbf{A}$. We
can select $B$ ($\geq 2M$) matrix elements of $\mathbf{A}$, either arbitrarily
or according to some prior information e.g. the elements with the largest mean
based on a statistical model. The corresponding $B$ elements of
$\mathbf{\tilde{A}}_{m}$ are then ordered into a column vector
$\mathbf{b}_{m}\in\mathbb{C}^{B}$ for every $m=1..2M$. These column vectors
form a matrix $\mathbf{B}_{\alpha}\in\mathbb{C}^{B\times 2M}$:

$\mathbf{B}_{\alpha}=\left[\begin{array}[]{c c c
c}\mathbf{b}_{1}&\mathbf{b}_{2}&\cdots&\mathbf{b}_{2M}\end{array}\right]$ (62)

so that we can write:

$\mathbf{B}_{\alpha}\mathbf{w}^{\top}=\mathbf{b}_{est}$ (63)

where $\mathbf{b}_{est}$ is an estimate of the $B$ selected elements of the
true TM, $\mathbf{A}$, and $\mathbf{w}=[w_{1}~{}~{}\cdots~{}~{}w_{N}]$ is a
vector containing the complex weights of Equation 61. Since
$\mathbf{B}_{\alpha}$ is either a square or tall matrix, we pre-multiply by
its Moore-Penrose pseudo inverse, $\mathbf{B}_{\alpha}^{\dagger}$:

$\mathbf{B}_{\alpha}^{\dagger}\mathbf{B}_{\alpha}\mathbf{w}^{\top}=\mathbf{B}_{\alpha}^{\dagger}\mathbf{b}_{est}$

$\mathbf{w}^{\top}=\mathbf{B}_{\alpha}^{\dagger}\mathbf{b}_{est}$

We then multiply both sides by $\mathbf{B}_{\alpha}$ to get:

$\mathbf{B}_{\alpha}\mathbf{w}^{\top}=\mathbf{B}_{\alpha}\mathbf{B}_{\alpha}^{\dagger}\mathbf{b}_{est}$

and substitute in Equation 63 to obtain a recursive expression:

$\mathbf{b}_{est}=\mathbf{B}_{\alpha}\mathbf{B}_{\alpha}^{\dagger}\mathbf{b}_{est}$
(64)

Following the same steps to derive Equation 64 from Equation 52 but starting
from Equation 53 gives:

$\mathbf{b}_{est}=\mathbf{B}_{\beta}\mathbf{B}_{\beta}^{\dagger}\mathbf{b}_{est}$
(65)

Similarly, starting from Equation 54 gives:

$\mathbf{b}_{est}=\mathbf{B}_{\gamma}\mathbf{B}_{\gamma}^{\dagger}\mathbf{b}_{est}$
(66)

Clearly, $\mathbf{b}_{est}$ is an eigenvector of
$\mathbf{B}_{\alpha}\mathbf{B}_{\alpha}^{\dagger}$,
$\mathbf{B}_{\beta}\mathbf{B}_{\beta}^{\dagger}$ and
$\mathbf{B}_{\gamma}\mathbf{B}_{\gamma}^{\dagger}$ with an eigenvalue of 1.
Physical considerations of power conservation suggest the existence of at
least one non-trivial solution for $\mathbf{b}_{est}$. Our empirical
investigations with both simulated and real TMs (see Section IV) suggest that
when using at least two out of Equations 64, 65 and 66 we can identify a
single non-trivial solution, the ‘true’ solution: we arbitrarily choose
Equations 64 and 65. Next we apply a variant of the ‘power method’ for finding
dominant eigenvalues of a matrix Bodewig (1959) and recursively substitute
Equation 64 into Equation 65, substitute the result into Equation 64 etc. to
give:

$\mathbf{b}_{est}^{t}=\mathbf{B}_{\alpha}\mathbf{B}_{\alpha}^{\dagger}\mathbf{B}_{\beta}\mathbf{B}_{\beta}^{\dagger}\dots\mathbf{B}_{\alpha}\mathbf{B}_{\alpha}^{\dagger}\mathbf{b}_{est}^{0}$
(67)

where $t$ denotes the number of matrix multiplications. For large $t$, the
output $\mathbf{b}_{est}^{t}$ is expected to converge to the desired solution
$\mathbf{b}_{est}$ for any input $\mathbf{b}_{est}^{0}$. Alternatively, for
large $t$, the dominant eigenvector of
$\mathbf{B}_{\alpha}\mathbf{B}_{\alpha}^{\dagger}\mathbf{B}_{\beta}\mathbf{B}_{\beta}^{\dagger}\dots\mathbf{B}_{\alpha}\mathbf{B}_{\alpha}^{\dagger}$
will approximate the solution well. In this work, we find empirically that 2–3
iterations, i.e. $t\in\\{4,6\\}$, is sufficient for faithful TM recovery.

To obtain the desired full TM, $\mathbf{A}$, we determine the weights,
$\mathbf{w}$, from approximation $\mathbf{b}_{est}^{t}$ using Equation 63 and
then compute the appropriate weighted sum of $\mathbf{A}_{m}$ with $m=1..2M$
using Equation 61.

An alternative to such a method could be to directly optimise the weights of
Equation 61 to best satisfy Equations 48–50, or to minimise a prior
probability distribution placed on the matrix (e.g. an empirically determined
sparsity structure such as that presented in Gordon  _et al._ (2019)).

###  B.1 Sylvester equations element by element using Bartels-Stewart
algorithm

The Bartels-Stewart algorithm is a method of solving matrix equations of the
form $\mathbf{D}\mathbf{X}-\mathbf{X}\mathbf{E}=\mathbf{F}$ (where
$\mathbf{D}$, $\mathbf{E}$, $\mathbf{F}$ $\in\mathbb{C}^{N\times N}$), termed
_Sylvester equations_ , to find the matrix $\mathbf{X}\in\mathbb{C}^{N\times
N}$ Bartels and Stewart (1972). Here, we wish to solve the special case where
$\mathbf{F}=\mathbf{0}$. We begin with the Schur decompositions of the
matrices $\mathbf{R}_{\alpha}$ and $\mathbf{C}_{\alpha}$ (adapted from
Equations 56 and 57):

$\mathbf{R}=\mathbf{Q}_{R}\mathbf{U}_{R}\mathbf{Q}_{R}^{-1}$ (68)

$\mathbf{C}=\mathbf{Q}_{C}\mathbf{U}_{C}\mathbf{Q}_{C}^{-1}$ (69)

where $\mathbf{Q}_{R}$ and $\mathbf{Q}_{R}$ is unitary and $\mathbf{U}_{R}$
and $\mathbf{U}_{C}$ are upper triangular matrices. Importantly, the diagonal
of $\mathbf{U}_{C}$ comprises the eigenvalues of $\mathbf{C}$ sorted in
descending order of magnitude, and likewise for $\mathbf{U}_{R}$. The Schur
decomposition of a matrix with entries sorted in this way can be computed, for
example, in MATLAB. This enables us to write the following Sylvester equation:

$\mathbf{U}_{R}\mathbf{A}^{\prime}-\mathbf{A}^{\prime}\mathbf{U}_{C}=\mathbf{0}$
(70)

where $\mathbf{A}^{\prime}$ is linked to the original transmission matrix,
$\mathbf{A}$ by:

$\mathbf{A}^{\prime}=\mathbf{Q}_{R}^{-1}\mathbf{A}\mathbf{Q}_{C}$ (71)

Recalling that all matrices here are $\in\mathbb{C}^{N\times N}$, we first
consider the element in row $N$, column $1$, i.e. $(N,1)$, of the zero matrix
on the RHS of Equation 70. We see that:

$ur_{N,N}a_{N,1}-a_{N,1}uc_{1,1}=0$ (72)

where $ur_{m,n}$ is the element in row $m$, column $n$ of $\mathbf{U}_{R}$,
$uc_{m,n}$ is the element in row $m$, column $n$ of $\mathbf{U}_{C}$, and
$a_{m,n}$ is the element in row $m$, column $n$ of $\mathbf{A}^{\prime}$.

Next, we require that the eigenvalues of $\mathbf{R}$ are distinct, which is
ensured by suitable design of reflectors (see Sections II.1 and III.1). We
also know that the eigenvalues of $\mathbf{R}$ are equal to those of
$\mathbf{C}_{\alpha}$ because the two are similar matrices (see, for example,
Equation 51). Therefore, the values on the diagonal of $\mathbf{U}_{R}$ are
distinct from one another, and are equal to the values on the diagonal of
$\mathbf{U}_{C}$. The only way that Equation 72 can hold is if $a_{N,1}=0$,
providing the solution for that element.

We then consider element $(N-1,1)$, giving:

$ur_{N-1,N-1}a_{N-1,1}+ur_{N-1,N}a_{N,1}-a_{N-1,1}uc_{1,1}=0$

Since we know that $a_{N,1}=0$, we apply again the above reasoning and
conclude that $a_{N-1,1}=0$. Continuing on up this column, we find that all
elements must be zero until we reach the first row, where:

$ur_{1,1}a_{1,1}-a_{1,1}uc_{1,1}=0$

Now, we know that $ur_{m,m}=uc_{m,m}$ for every $m$ as discussed above.
Therefore, every possible $a_{1,1}\in\mathbb{C}$ satisfies this equation,
meaning the element is indeterminate. We now consider the second column of the
zero matrix on the RHS of Equation 70, starting with element $(N,2)$:

$ur_{N,N}a_{N,2}-(a_{N,1}uc_{1,2}+a_{N,2}uc_{2,2})=0$

We know $a_{N,1}=0$ so can write:

$ur_{N,N}a_{N,2}-a_{N,2}uc_{2,2}=0$

and therefore $a_{N,2}=0$. Next, we consider element $(N-1,2)$:

$ur_{N-1,N-1}a_{N-1,2}+ur_{N-1,N}a_{N,2}-(a_{N-1,1}uc_{1,2}+a_{N-1,2}uc_{2,2})=0$

Using previously known elements, we find that:

$ur_{N-1,N-1}a_{N-1,2}-a_{N-1,2}uc_{2,2}=0$

Therefore $a_{N-1,2}=0$. This continues up column 2 of the zero matrix until
we get to:

$ur_{2,2}a_{2,2}-a_{2,2}uc_{2,2}=0$

$a_{2,2}$ can be anything since $ur_{2,2}=uc_{2,2}$. Now, we consider element
$(1,2)$:

$ur_{1,1}a_{1,2}+ur_{1,2}a_{2,2}-(a_{1,1}uc_{1,2}+a_{1,2}uc_{2,2})=0$

$(ur_{1,1}-uc_{2,2})a_{1,2}+ur_{1,2}a_{2,2}-a_{1,1}uc_{1,2}=0$

By repeatedly applying this logic to each column of the zero matrix on the RHS
of Equation 70, we conclude that the matrix $\mathbf{A}^{\prime}$ must have an
upper triangular form. It then follows that element $(3,3)$ gives:

$ur_{3,3}a_{3,3}-a_{3,3}uc_{3,3}=0$

and element $(2,3)$ gives:

$ur_{2,1}a_{1,3}+ur_{2,2}a_{2,3}+ur_{2,3}a_{3,3}-(a_{2,1}uc_{1,3}+a_{2,2}uc_{2,3}+a_{2,3}uc_{3,3})=0$

$ur_{2,1}a_{1,3}+(ur_{2,2}-uc_{3,3})a_{2,3}+ur_{2,3}a_{3,3}-a_{2,2}uc_{2,3}=0$

Considering element $(1,3)$ gives:

$ur_{1,1}a_{1,3}+ur_{1,2}a_{2,3}+ur_{1,3}a_{3,3}-(a_{1,1}uc_{1,3}+a_{1,2}uc_{2,3}+a_{1,3}uc_{3,3})=0$

$(ur_{1,1}-uc_{3,3})a_{1,3}+ur_{1,2}a_{2,3}+ur_{1,3}a_{3,3}-a_{1,1}uc_{1,3}-a_{1,2}uc_{2,3}=0$

Therefore, these elements are all functions of elements that have already been
computed when solving previously columns of the zero matrix on the RHS of
Equation 70, and the indeterminate diagonal elements of $\mathbf{A}^{\prime}$.
By repeated application of this process, we can generate an equation for each
upper triangular element, $(P,Q)$ with $P<Q$, of the zero matrix on the RHS of
Equation 70 as:

$(ur_{P,P}-uc_{Q,Q})a_{P,Q}+\sum_{q=P+1}^{Q}ur_{P,q}a_{q,Q}-\sum_{s=P}^{Q-1}uc_{s,Q}a_{P,s}=0$

Rearranging, we can write:

$\left[-\sum_{s=P}^{Q-1}uc_{s,Q}a_{P,s}\right]+(ur_{P,P}-uc_{Q,Q})a_{P,Q}+\sum_{q=P+1}^{Q}ur_{P,q}a_{q,Q}=0$

$\left[ur_{P,Q}a_{Q,Q}-\sum_{s=P}^{Q-1}uc_{s,Q}a_{P,s}\right]+(ur_{P,P}-uc_{Q,Q})a_{P,Q}+\sum_{q=P+1}^{Q-1}ur_{P,q}a_{q,Q}=0$
(73)

In summary, we find that the diagonal elements, $a_{m,m}$ with $m=1..N$, are
indeterminate but that all other elements in the matrix $\mathbf{A}^{\prime}$
can be computed from these diagonals by solving a series of linear equations.

##  Appendix C Orthogonal basis generation

To achieve maximum signal-to-noise ratio during optimisation or any other
solution method, the basis generated for reconstructing the TM, $\mathbf{A}$
in Section II.2, should ideally be orthogonal. This basis is constructed
starting from an arbitrary vector $\in\mathbb{C}^{N}$ used for the diagonal
elements (Figure 14). In general, this will produce a complete but non-
orthogonal basis. Elements of this basis may be very similar resulting in
numerical instabilities during the zeroth-order recovery process.

To produce an orthogonal basis, we first define a function that gives the
solution of the system of equations defined in Equation 60. This function,
which could be implemented by any linear equation solver, takes as input a
vector representing diagonal elements of the solution matrix, $\mathbf{u}$,
and returns the full matrix, $\mathbf{B}$:

$\mathbf{B}=f(\mathbf{u})$ (74)

We can then start with an arbitrary set of complex diagonal elements,
$\mathbf{u}_{1}\in\mathbb{C}^{N}$, to get a basis element for
$\mathbf{A}^{\prime}$ (defined in Equation 71), which we term
$\mathbf{B}_{1}^{\prime}=f(\mathbf{u}_{1})$. Next, we invert Equation 71 to
find the associated basis element for $\mathbf{A}$, given by
$\mathbf{B}_{1}=\mathbf{Q}_{R}\mathbf{B}_{1}^{\prime}\mathbf{Q}_{C}^{-1}$. To
get the next basis element we then generate a random complex matrix,
$\mathbf{B}_{2}\in\mathbb{C}^{N\times N}$ and perform a Gram-Schmidt
orthogonalisation step:

$\mathrm{vec}(\mathbf{\hat{B}}_{2})=\mathrm{vec}(\mathbf{B}_{2})-\mathrm{proj}_{\mathrm{vec}(\mathbf{B}_{1})}\mathrm{vec}(\mathbf{B}_{2})$
(75)

where $\mathrm{vec}(..)$ represents ordering of matrix elements into a vector.
We then use Equation 71 to convert back to an estimated basis element for
$\mathbf{A}^{\prime}$:

$\mathbf{\hat{B}}_{2}^{\prime}=\mathbf{Q}_{R}^{-1}\mathbf{\hat{B}}_{2}\mathbf{Q}_{C}$
(76)

and find the diagonal elements of this basis element:

$\mathbf{\hat{u}}_{2}=\mathrm{diag}\left(\mathbf{\hat{B}}_{2}^{\prime}\right)$
(77)

We then use these diagonal elements to generate a new matrix that is a valid
solution to 60 (i.e. a basis element for $\mathbf{A}^{\prime}$):

$\mathbf{B}_{2}^{\prime}=f(\mathbf{\hat{u}}_{2})$ (78)

and convert back to a basis element for $\mathbf{A}$:

$\mathbf{B}_{2}=\mathbf{Q}_{R}\mathbf{B}_{2}^{\prime}\mathbf{Q}_{C}^{-1}$ (79)

We then substitute $\mathbf{B}_{2}$ back into Equation 75 and iterate several
times. The result of this will be a basis matrix, $\mathbf{B}_{2}$, that is
orthogonal to $\mathbf{B}_{1}$. This process is then repeated for subsequent
basis elements using all previously generated basis elements. For example, the
Gram-Schmidt orthogonalisation procedure step for the $m^{th}$ basis element
becomes:

$\mathrm{vec}(\mathbf{\hat{B}}_{m})=\mathrm{vec}(\mathbf{B}_{m})-\sum_{j=1}^{m}\mathrm{proj}_{\mathrm{vec}(\mathbf{B}_{j})}\mathrm{vec}(\mathbf{B}_{m})$
(80)

In this way, an orthogonal basis for $\mathbf{A}$ can be generated.

##  Appendix D Alternative method of finding the null space of Sylvester
equations

In some cases, finding a basis for the TM can suffer from numerical error due
to the large number of computation steps involved (see Section II.2). To
derive an alternative approach, we begin with Equation 52:

$\mathbf{R}_{\alpha}\mathbf{A}-\mathbf{A}\mathbf{C}_{\alpha}=\mathbf{0}$ (81)

We can rewrite this in the following form Horn and Johnson (2013):

$\mathbf{M}~{}~{}\mathrm{vec}(\mathbf{A})=0$ (82)

where $\mathrm{vec}(..)$ represents ordering of matrix elements into a vector
and

$\mathbf{M}=\mathbf{I}_{N\times
N}\otimes\mathbf{R}_{\alpha}+\mathbf{R}_{\alpha}\otimes\mathbf{I}_{N\times N}$
(83)

where $\otimes$ is the Kronecker product. By finding the nullspace of
$\mathbf{M}$ we find an orthogonal basis for $\mathbf{A}$. If $\mathbf{A}$ is
of size $N\times N$, the resultant matrix $\mathbf{M}$ is of size $N^{2}\times
N^{2}$ and so for larger problems (e.g. the 1648$\times$1648 MCF TMs of
Section IV.2), the memory usage becomes impractically large. Therefore, in
this work this method is used for small simulated matrices.

##  Appendix E Approximate bandwidth calculation

Let us call the centre wavelength of the fibre $\lambda_{0}$ and the bandwidth
of interest $\Delta\lambda$. The bandwidth in frequency units is:

$f_{1}-f_{2}=\frac{c}{n_{eff}(\lambda_{0}-\Delta\lambda/2)}-\frac{c}{n_{eff}(\lambda_{0}+\Delta\lambda/2)}$

$\Delta
f=\frac{c}{n_{eff}}\left(\frac{1}{\lambda_{0}-\Delta\lambda/2}-\frac{1}{\lambda_{0}+\Delta\lambda/2}\right)$
(84)

where $n_{eff}$ is the effective refractive index of the fibe (e.g. $\sim$1.5
for glass) and $c$ is the speed of light in a vacuum. This bandwidth is the
range over which the eigenmodes of the fibre remain constant, which is
equivalent to the spectral correlation bandwidth Mosk  _et al._ (2012). This
is in turn related to the time difference between longest path lengths (time
of flight) by:

$\Delta t=\frac{1}{\Delta
f}=\frac{n_{eff}}{c}\frac{1}{\left(\frac{1}{\lambda_{0}-\Delta\lambda/2}-\frac{1}{\lambda_{0}+\Delta\lambda/2}\right)}$
(85)

If we assume the fibre has a physical length $\ell$, we find that time of
flight down the fibre for the fastest mode will be of the order of:

$t_{1}=\frac{\ell}{\frac{c}{n_{eff}}}$ (86)

The longest time of flight will then be:

$t_{2}=t_{1}+\Delta t$ (87)

This results in an effective length for the longest path of:

$\ell_{2}=\frac{c}{n_{eff}}t_{2}=\frac{c}{n_{eff}}\left(t_{1}+\Delta t\right)$
(88)

$\ell_{2}=\frac{c}{n_{eff}}\left(\frac{\ell}{\frac{c}{n_{eff}}}+\Delta
t\right)=\ell+\frac{c}{n_{eff}}\Delta t$ (89)

$\ell_{2}=\ell+\frac{1}{\left(\frac{1}{\lambda_{0}-\Delta\lambda/2}-\frac{1}{\lambda_{0}+\Delta\lambda/2}\right)}$
(90)

At the shortest wavelength, $\lambda_{0}-\Delta\lambda/2$, the phase shift
introduced by the shortest path is:

$\phi_{1,short}=\frac{2\pi\ell}{\lambda_{0}-\Delta\lambda/2}$ (91)

and by the longest path:

$\phi_{1,long}=\frac{2\pi\ell_{2}}{\lambda_{0}-\Delta\lambda/2}$ (92)

$\phi_{1,long}=\frac{2\pi\left(\ell+\frac{1}{\left(\frac{1}{\lambda_{0}-\Delta\lambda/2}-\frac{1}{\lambda_{0}+\Delta\lambda/2}\right)}\right)}{\lambda_{0}-\Delta\lambda/2}$
(93)

giving a difference of:

$\Delta\phi_{1}=\phi_{1,long}-\phi_{1,short}=2\pi\frac{\frac{1}{\left(\frac{1}{\lambda_{0}-\Delta\lambda/2}-\frac{1}{\lambda_{0}+\Delta\lambda/2}\right)}}{\lambda_{0}-\Delta\lambda/2}$
(94)

$\Delta\phi_{1}=2\pi\frac{1}{\left(1-\frac{\lambda_{0}-\Delta\lambda/2}{\lambda_{0}+\Delta\lambda/2}\right)}$
(95)

At the longest wavelength, $\lambda_{0}+\Delta\lambda/2$, the phase shift
introduced by the shortest path is:

$\phi_{2,short}=\frac{2\pi\ell}{\lambda_{0}+\Delta\lambda/2}$ (96)

and by the longest path:

$\phi_{2,long}=\frac{2\pi\ell_{2}}{\lambda_{0}+\Delta\lambda/2}$ (97)

$\phi_{2,long}=\frac{2\pi\left(\ell+\frac{1}{\left(\frac{1}{\lambda_{0}-\Delta\lambda/2}-\frac{1}{\lambda_{0}+\Delta\lambda/2}\right)}\right)}{\lambda_{0}+\Delta\lambda/2}$
(98)

giving a difference of:

$\Delta\phi_{2}=\phi_{2,long}-\phi_{2,short}=2\pi\frac{\frac{1}{\left(\frac{1}{\lambda_{0}-\Delta\lambda/2}-\frac{1}{\lambda_{0}+\Delta\lambda/2}\right)}}{\lambda_{0}+\Delta\lambda/2}$
(99)

$\Delta\phi_{2}=2\pi\frac{1}{\left(\frac{\lambda_{0}+\Delta\lambda/2}{\lambda_{0}-\Delta\lambda/2}-1\right)}$
(100)

Now we consider the difference in phase shift between the longest and shortest
path at the two wavelengths, which corresponds to the phase difference between
respective elements of matrices at the two different wavelengths. This gives:

$\Delta\phi_{1}-\Delta\phi_{2}=2\pi\left(\frac{1}{\left(1-\frac{\lambda_{0}-\Delta\lambda/2}{\lambda_{0}+\Delta\lambda/2}\right)}-\frac{1}{\left(\frac{\lambda_{0}+\Delta\lambda/2}{\lambda_{0}-\Delta\lambda/2}-1\right)}\right)$
(101)

$\Delta\phi_{1}-\Delta\phi_{2}=2\pi\left(\frac{1}{\left(1-\frac{\lambda_{0}-\Delta\lambda/2}{\lambda_{0}+\Delta\lambda/2}\right)}+\frac{1}{\left(1-\frac{\lambda_{0}+\Delta\lambda/2}{\lambda_{0}-\Delta\lambda/2}\right)}\right)$
(102)

$\Delta\phi_{1}-\Delta\phi_{2}=2\pi\left(\frac{\lambda_{0}+\Delta\lambda/2}{\left(\lambda_{0}+\Delta\lambda/2-\lambda_{0}+\Delta\lambda/2\right)}+\frac{\lambda_{0}-\Delta\lambda/2}{\left(\lambda_{0}-\Delta\lambda/2-\lambda_{0}-\Delta\lambda/2\right)}\right)$
(103)

$\Delta\phi_{1}-\Delta\phi_{2}=2\pi\left(\frac{\lambda_{0}+\Delta\lambda/2}{\Delta\lambda}+\frac{\lambda_{0}-\Delta\lambda/2}{-\Delta\lambda}\right)$
(104)

$\Delta\phi_{1}-\Delta\phi_{2}=2\pi\left(\frac{\lambda_{0}+\Delta\lambda/2}{\Delta\lambda}-\frac{\lambda_{0}-\Delta\lambda/2}{\Delta\lambda}\right)$
(105)

$\Delta\phi_{1}-\Delta\phi_{2}=2\pi$ (106)

This then shows that if the wavelength is kept within the spectral bandwidth,
the maximum additional phase shift introduced between two paths will be
$2\pi$, validating the assumption made when taking the matrix logarithm in
Section II.3.

Intuitively, if path lengths vary in phase by $\ll 2\pi$ when the wavelength
is varied, the diffracted far-field pattern will remain similar and so the
‘speckle correlation’ will be high. When the path length differs by $>2\pi$,
variations in wavelengths produce quite different diffraction patterns (as the
phases will be significantly altered) and thus the speckle patterns will
become decorrelated. This is described in more detail in Mosk  _et al._
(2012).

##  Appendix F Bandwidth of step-index MMF

To validate the first-order recovery algorithm, we use a dataset consisting of
TMs taken at multiple wavelengths (1525–1567nm in steps of 0.08nm) Carpenter
(2016). The TMs are taken from a step-index refractive index profile MMF with
420 modes (including polarisation modes). We first examine the bandwidth over
which the first-order model is approximately accurate. To do this, we consider
the TM as a function of wavelength, $\mathbf{A}_{\mathrm{MMF}}(\lambda)$,
where $\lambda$ is wavelength in nm. Using a measured TM at a reference
wavelength, $\lambda_{ref}$ (1526nm is used here), we estimate TM at some test
wavelength, $\lambda_{test}>\lambda_{ref}$ as:

$\mathbf{A}_{est}(\lambda_{test},\lambda_{ref})=e^{\left(\frac{\lambda_{ref}}{\lambda_{test}}\log\mathbf{A}_{\mathrm{MMF}}(\lambda_{ref})\right)}$
(107)

We then compute an error metric to test the predictive accuracy of the first-
order model:

$\epsilon_{est}(\lambda_{test},\lambda_{ref})=\sqrt{\frac{\left\lVert\mathbf{A}_{est}(\lambda_{test},\lambda_{ref})-\mathbf{A}_{\mathrm{MMF}}(\lambda_{test})\right\rVert_{F}}{\left\lVert\mathbf{A}_{\mathrm{MMF}}(\lambda_{test})\right\rVert_{F}}}$
(108)

where $\mathbf{A}_{\mathrm{MMF}}(\lambda_{test})$ is the measured MMF TM at
wavelength $\lambda_{test}$ and $\left\lVert..\right\rVert_{F}$ is the
Frobenius norm. Defining $\Delta\lambda=\lambda_{test}-\lambda_{ref}$ we can
rewrite Equation 108 as:

$\epsilon_{est}(\Delta\lambda)=\sqrt{\frac{\left\lVert\mathbf{A}_{est}(\Delta\lambda,\lambda_{ref})-\mathbf{A}_{\mathrm{MMF}}(\lambda_{ref}+\Delta\lambda)\right\rVert_{2}}{\left\lVert\mathbf{A}_{\mathrm{MMF}}(\lambda_{ref}+\Delta\lambda)\right\rVert_{2}}}$
(109)

Setting $\lambda_{ref}$ to 1526nm, we can then plot $\epsilon_{est}$ as a
function of $\Delta\lambda$ (Figure 15). We find empirically that changing
$\lambda_{ref}$ does not significantly change this curve.

The full $420\times 420$ matrix provides $\epsilon_{est}=0.5$ for a
$\Delta\lambda$ of $\sim$0.5nm. This is insufficient to to enable the design
of a reflector stack like that in Section III.1 using off-the-shelf optical
filters because the spacings between filter centre wavelengths would need to
be $<0.1$nm. To improve bandwidth, we restrict analysis to the lowest order
modes of the fibre and form an $N\times N$ sub-matrix comprising the first $N$
columns of the first $N$ rows of each $\mathbf{A}_{\mathrm{MMF}}$. The
resulting $\epsilon_{est}$ curves for $N=30,72,110,132,156,210,420$ are
plotted in Figure 15. We find that using the 110 lowest-order modes gives a
$\Delta\lambda$ of 5nm for $\epsilon_{est}=0.5$. This is ample to allow design
of a realistic filter stack in line with the principles of Section III.1 and
III.2 and so is used here.

Figure 15: Validation of the first-order model of Section A.1 using measured
multi-wavelength TMs from a step-index MMF Carpenter (2016). The normalised
error of the TM predicted by the first-order model compared with the actual
measured TM ($\epsilon_{est}$) is plotted as the difference in reference and
test wavelengths ($\Delta\lambda$) is increased. Different sub-matrices of
only the lowest-order modes are considered: that is, sub-matrices of size
$N\times N$ with $N=30,72,110,132,156,210,420$. For lower-order sub-matrices
the first-order model is seen to be valid over a larger wavelength range.

In reality dispersion will further change the effective optical thicknesses at
the different wavelengths. From the dataset here, we compute that over a 5nm
range the error introduced is $<$3% and so is therefore neglected in this
work.

## References

  * Cižmár and Dholakia (2012) T. Cižmár and K. Dholakia, Exploiting multimode waveguides for pure fibre-based imaging., Nature Communications 3, 1027 (2012).
  * Gataric  _et al._ (2019) M. Gataric, G. S. D. Gordon, F. Renna, A. G. C. P. Ramos, M. P. Alcolea, and S. E. Bohndiek, Reconstruction of Optical Vector-Fields With Applications in Endoscopic Imaging, IEEE Transactions on Medical Imaging 38, 955 (2019), arXiv:1804.10636 .
  * Gordon  _et al._ (2018) G. S. D. Gordon, J. Joseph, M. P. Alcolea, T. Sawyer, A. J. Macfaden, C. Williams, C. R. M. Fitzpatrick, P. H. Jones, M. di Pietro, R. C. Fitzgerald, T. D. Wilkinson, and S. E. Bohndiek, Quantitative phase and polarisation endoscopy applied to detection of early oesophageal tumourigenesis,  (2018), arXiv:1811.03977 .
  * Loterie  _et al._ (2015) D. Loterie, S. Farahi, I. Papadopoulos, A. Goy, D. Psaltis, and C. Moser, Digital confocal microscopy through a multimode fiber, Optics Express 23, 23845 (2015), arXiv:1502.04172 .
  * Rotter and Gigan (2017) S. Rotter and S. Gigan, Light fields in complex media: Mesoscopic scattering meets wave control, Reviews of Modern Physics 89, (2017), arXiv:1702.05395 .
  * Dunning and Lind (1982) G. J. Dunning and R. C. Lind, Demonstration of image transmission through fibers by optical phase conjugation., Optics letters 7, 558 (1982).
  * Di Leonardo and Bianchi (2011) R. Di Leonardo and S. Bianchi, Hologram transmission through multi-mode optical fibers., Optics express 19, 247 (2011).
  * Vasquez-Lopez  _et al._ (2018) S. A. Vasquez-Lopez, R. Turcotte, V. Koren, M. Plöschner, Z. Padamsey, M. J. Booth, T. Čižmár, and N. J. Emptage, Subcellular spatial resolution achieved for deep-brain imaging in vivo using a minimally invasive multimode fiber, Light: Science & Applications 7, 110 (2018).
  * Farahi  _et al._ (2013) S. Farahi, D. Ziegler, I. N. Papadopoulos, D. Psaltis, and C. Moser, Dynamic bending compensation while focusing through a multimode fiber, Optics Express 21, 22504 (2013).
  * Stasio  _et al._ (2015) N. Stasio, D. B. Conkey, C. Moser, and D. Psaltis, Light control in a multicore fiber using the memory effect, Optics Express 23, 30532 (2015).
  * Plöschner  _et al._ (2015) M. Plöschner, T. Tyc, and T. Čižmár, Seeing through chaos in multimode fibres, Nature Photonics 9, 529 (2015).
  * Spillman and Udd (2014) W. B. Spillman and E. Udd,  _Field Guide to Fiber Optic Sensors_ , Smart Sensors, Measurement and Instrumentation, Vol. 21 (SPIE, Bellingham, WA, 2014).
  * Udd (2011) E. Udd, Fiber Shape Sensing Systems and Methods (2011).
  * Chen  _et al._ (2018) H. Chen, N. K. Fontaine, R. Ryf, D. T. Neilson, P. Winzer, S. Gross, N. Riesen, M. Withford, J. Li, J. C. Alvarado, S. Wittek, J. Du, Z. He, J. E. Antonio Lopez, and R. Amezcua-Correa, Remote Mode-Forming Over Multimode Fiber Employing Single-Ended Channel Estimation, in  _2018 European Conference on Optical Communication (ECOC)_, Vol. 1 (IEEE, 2018) pp. 1–3.
  * Gu  _et al._ (2015) R. Y. Gu, R. N. Mahalati, and J. M. Kahn, Design of flexible multi-mode fiber endoscope, Optics Express 23, 26905 (2015).
  * Lee  _et al._ (2018) S. Lee, B. E. Bouma, and M. Villiger, Reciprocity in Measuring Multimode Fiber Transmission, in  _2018 IEEE Photonics Conference (IPC)_ (IEEE, 2018) pp. 1–2.
  * Warren  _et al._ (2016) S. C. Warren, Y. Kim, J. M. Stone, C. Mitchell, J. C. Knight, M. A. A. Neil, C. Paterson, P. M. W. French, and C. Dunsby, Adaptive multiphoton endomicroscopy through a dynamically deformed multicore optical fiber using proximal detection, Opt. Express 24, 21474 (2016).
  * Gu (2017) R. Y. Gu,  _Rigid and Flexible Multimode Fiber Endoscopes_ , Ph.D. thesis, Stanford University (2017).
  * Carpenter  _et al._ (2014) J. Carpenter, B. Eggleton, and J. Schröder, 110x110 optical mode transfer matrix inversion, Optics Express 22, 96 (2014).
  * Gordon  _et al._ (2014) G. Gordon, M. Crisp, R. Penty, T. Wilkinson, and I. White, Feasibility Demonstration of a Mode-Division Multiplexed MIMO-enabled Radio-over-Fiber Distributed Antenna System, Journal of Lightwave Technology 8724, 1 (2014).
  * Fan and Kahn (2005) S. Fan and J. M. Kahn, Principal modes in multimode waveguides., Optics letters 30, 135 (2005).
  * Goldstein (2010) D. H. Goldstein,  _Polarized Light_ , 3rd ed. (Boca Ranton, 2010).
  * Redding  _et al._ (2013) B. Redding, S. M. Popoff, and H. Cao, All-fiber spectrometer based on speckle pattern reconstruction, Optics Express 21, 6584 (2013).
  * Yariv (1973) A. Yariv, Coupled-mode theory for guided-wave optics, IEEE Journal of Quantum Electronics 9, 919 (1973).
  * Huang (1994) W. Huang, Coupled-mode theory for optical waveguides: an overview, Journal of the Optical Society of America A 11, 963 (1994), arXiv:1502.06504v1 .
  * Loterie  _et al._ (2017) D. Loterie, D. Psaltis, and C. Moser, Bend translation in multimode fiber imaging, Optics Express 25, 6263 (2017).
  * Dragt (1982) A. J. Dragt, Lie algebraic theory of geometrical optics and optical aberrations, Journal of the Optical Society of America 72, 372 (1982).
  * Payne and Gambling (1975) D. N. Payne and W. A. Gambling, Zero material dispersion in optical fibre, Electronics Letters 11, 176 (1975).
  * Mosk  _et al._ (2012) A. P. Mosk, A. Lagendijk, G. Lerosey, and M. Fink, Controlling waves in space and time for imaging and focusing in complex media, Nature Photonics 6, 283 (2012).
  * Hall (2015) B. C. Hall,  _Lie groups, Lie algebras, and representations: An elementary introduction_ , 2nd ed. (Springer, 2015).
  * Carpenter  _et al._ (2015) J. Carpenter, B. J. Eggleton, and J. Schröder, Observation of Eisenbud–Wigner–Smith states as principal modes in multimode fibre, Nature Photonics 9, 751 (2015).
  * Choi  _et al._ (2012) 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, Physical Review Letters 109, 203901 (2012).
  * Brenner  _et al._ (1976) J. F. Brenner, B. S. Dew, J. B. Horton, T. King, P. W. Neurath, and W. D. Selles, An automated microscope for cytologic research a preliminary evaluation., Journal of Histochemistry & Cytochemistry 24, 100 (1976).
  * Cuccia  _et al._ (2009) D. J. Cuccia, F. Bevilacqua, A. J. Durkin, F. R. Ayers, and B. J. Tromberg, Quantitation and mapping of tissue optical properties using modulated imaging, Journal of Biomedical Optics 14, 024012 (2009).
  * Williams  _et al._ (2016) C. Williams, R. Bartholomew, G. Rughoobur, G. S. D. Gordon, A. J. Flewitt, and T. D. Wilkinson, Fabrication of nanostructured transmissive optical devices on ITO-glass with UV1116 photoresist using high-energy electron beam lithography, Nanotechnology 27, 485301 (2016).
  * Thompson (2012) A. Thompson,  _Developing endoscopic instrumentation and techniques for in vivo fluorescence lifetime imaging and spectroscopy_ , Ph.D. thesis, Imperial College London (2012).
  * Gordon  _et al._ (2019) G. Gordon, R. Mouthaan, T. Wilkinson, and S. Bohndiek, Coherent Imaging Through Multicore Fibres with Applications in Endoscopy, Under review (2019).
  * Carpenter (2016) J. Carpenter, A Matlab file containing the optical transfer function of a multimode optical fiber (2016).
  * Orfanidis (2016) S. J. Orfanidis,  _Electromagnetic Waves and Antennas_ (Rutgers University, Piscataway, 2016) pp. 216–223.
  * Goodman (1996) J. W. Goodman,  _Introduction To Fourier Optics_ , 2nd ed., McGraw-Hill physical and quantum electronics series (McGraw-Hill, New York, New York, USA, 1996).
  * Yamamoto (1983) Y. Yamamoto, AM and FM quantum noise in semiconductor lasers–Part I: Theoretical analysis, IEEE Journal of Quantum Electronics 19, 34 (1983).
  * Kingma and Ba (2014) D. P. Kingma and J. Ba, Adam: A Method for Stochastic Optimization,  (2014), arXiv:1412.6980 .
  * Rahmani  _et al._ (2018) B. Rahmani, D. Loterie, G. Konstantinou, D. Psaltis, and C. Moser, Multimode optical fiber transmission with a deep learning network, Light: Science & Applications 7, 69 (2018), arXiv:1805.05134 .
  * Yang  _et al._ (2018) M. Yang, Z.-h. Liu, Z.-d. Cheng, J.-s. Xu, C.-f. Li, and G.-c. Guo, Deep Hybrid Scattering Image Learning,  (2018), arXiv:1809.07706 .
  * Mitchell  _et al._ (2016) K. J. Mitchell, S. Turtaev, M. J. Padgett, T. Čižmár, and D. B. Phillips, High-speed spatial control of the intensity, phase and polarisation of vector beams using a digital micro-mirror device, Optics Express 24, 29269 (2016), arXiv:1607.08236 .
  * Williams  _et al._ (2017) C. Williams, G. Rughoobur, A. J. Flewitt, and T. D. Wilkinson, Nanostructured plasmonic metapixels, Scientific Reports 7, 7745 (2017).
  * Popoff  _et al._ (2011) S. M. Popoff, G. Lerosey, M. Fink, A. C. Boccara, and S. Gigan, Controlling light through optical disordered media: Transmission matrix approach, New Journal of Physics 13, 123021 (2011), arXiv:1107.5285 .
  * Neisser and Wurm (2015) M. Neisser and S. Wurm, ITRS lithography roadmap: 2015 challenges, Advanced Optical Technologies 4, (2015).
  * Bartels and Stewart (1972) R. Bartels and G. W. Stewart, Algorithm 432: Solution of the Matrix Equation {A X + X B = C}, Comm. ACM 15, 820 (1972).
  * Horn and Johnson (2013) R. A. Horn and C. R. Johnson,  _Book_ (Cambridge University Press, New York, 2013) p. 607.
  * Marčenko and Pastur (1967) V. A. Marčenko and L. A. Pastur, Distribution of Eigenvalues for Some Sets of Random Matrices, Mathematics of the USSR-Sbornik 1, 457 (1967).
  * Popoff  _et al._ (2010) 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, Physical Review Letters 104, 100601 (2010).
  * Bodewig (1959) E. Bodewig,  _Matrix Calculus_ (North Holland Publishing Company, Amsterdam, 1959).