\providecommand{\main}{..}

\documentclass[main.tex]{subfiles}
\graphicspath{ {../images/} }
\begin{document}
% Background section
\section{Background}

% Subs 1
\subsection{Microscopy}

Microscopy has been an area of study for centuries. Robert Hooke, in 1665 \cite{backg_micrographia}, presented illustrations of a new world unseen before, including detailed drawings of insects and plants. This book is also the first to coin the term cell, which is the foundation of biology. Robert Hooke’s work came a few decades after the Dutch lens-maker Hans Janssen and his son Zacharias are attributed to have invented the compound microscope \cite{backg_early_history_microscope}. In this contraption, light which passes through the sample forms a magnified image of the sample by using two lenses (ocular/eyepiece and objective), the individual magnification of the two lenses being multiplied.

The most common form of microscopy is widefield microscopy, in which the whole sample is illuminated. The image can then be formed either using light transmitted through or reflected from the sample. However, other modes of illumination have been introduced over time. One of the drawbacks of widefield microscopy is that parts of the sample which are out of focus are still imaged and observed in a blurred state, leading to decreased quality of images acquired. In order to reject out-of-focus light, confocal microscopy was introduced \cite{backg_confocal_microscopy}. The first step in reducing out-of-focus light is to restrict laser illumination to a single spot. This method allows for the light captured from the sample to then be refocused. Then, a pinhole acts as a spatial filter, as the centre of the refocus contains only light which originates from the in-focus position of the sample. Light originating from out-of-focus locations, focuses away from the centre and is blocked. While this method rejects out-of-focus light, only a small area is illuminated at any time. Thus, in order to image the entire sample, the illumination focus needs to be relocated to different positions and the different observations of the filtered light form together the entire image. This method is named scanning microscopy.

The earliest microscopy methods relied on observation of transmitted or reflected light. However, other methods such as fluorescence have become widely used. Fluorescence is a phenomenon in which a photon (or multiple) is absorbed, then after a few nanoseconds another photon is emitted, which has lower energy as some is lost in the process \cite{backg_fluorescence_microscopy}. The main advantage of this process is that it is dependent on molecules, named fluorophores. Thus, the process can be localised by having these fluorophores only in specific structures of the sample, imaging only specific areas of interest, such as a cell membrane for example. In addition, fluorophores can also be used for revealing the function of the cell and observe cellular activity \cite{backg_genetical_neuronal}.  The process can be implemented in microscopic systems such as the confocal microscope described above. The emitted light, which is at a different wavelength to the illumination, can then be filtered through the use of a dichroic mirror, separating the optical paths.

The emission of light from fluorescence is dependent on the density of photons in the illuminated area. Processes such as two-photon fluorescence are non-linearly dependent on the illumination intensity due to being a second-order process, as using light with double the wavelength compared with 1-photon fluorescence requires the use of two photons to excite the same fluorophore. The advantage of two-photon fluorescence is reduced out-of-focus emission and larger penetration depth as a longer wavelength is less susceptible to scattering.

%Subs 2
\subsection{\textit{In vivo} imaging}

In the past 100 years, neuroscience has strived to understand the processes taking place inside the brain which give rise to cognition. While these processes can be studied at behavioural and computational level, they are the output of a complex neural network \cite{backg_understanding_brain}. \textit{In vivo} cellular imaging has provided the means for observing individual neurons while keeping intact the surrounding tissue and connections. Advancements in imaging technologies have extended the limits of \textit{in vivo} observation to what was previously only achievable \textit{in vitro} \cite{backg_imaging_in_vivo}.

One of the limitations of \textit{in vivo} imaging is the strong light scattering of different tissue components. Multi-photon fluorescence differs from single-photon, by utilising multiple photons of longer wavelength to achieve the same amount of energy required to excite the fluorophore. The advantage is that longer wavelength light is less susceptible to scattering, hence can penetrate further into tissue than shorter wavelength light. With the use of two-photon fluorescence microscopy, imaging depths of up to 1.6 mm have been achieved in cortex vasculature \cite{backg_in_vivo_2Ph}. In the case of neuronal imaging, the focus power decreases drastically after $700 \; \mu m$ due to scattering in the white matter \cite{backg_1300_nm_3Ph}. Three-photon fluorescence provides a better alternative. The drawback is that for larger depths, in order to keep the focal energy high, as required for multi-photon fluorescence, cortex heating can pose a problem \cite{mmf_neuron}. A longer excitation wavelength also leads to a lower imaging resolution due to the diffraction limit. However, if scattering is the main limitation in imaging, using a waveguide to transmit the excitation light deep inside the brain eliminates the problem. The drawback is that this method is invasive. 

An alternative approach is to use an imaging probe that is inserted into the tissue which minimises invasiveness. The best choice for minimising the footprint of the probe is to use an optical fibre \cite{mmf_neuron}. The minimum invasiveness is of paramount importance in order for the biological system observed to perform as close as possible to homeostasis.

In order for a multimode optical fibre to be used for imaging, the transmission through it needs to be characterised, due to the light being scrambled by the fibre. This can be achieved by empirically measuring the fibre properties prior to implanting inside the brain \cite{mmf_neuron}. One of the limitations of optical fibres is the sensitivity of light transmission to changes in the fibre’s conformation. Thus, if the fibre bends, the prior calibration can become useless \cite{seeing_through_chaos}. Compensating for this limitation is a main requirement in the quest for deep brain \textit{in vivo} imaging in awake, behaving animals. Creating the technology necessary to enable the acquisition of \textit{in vivo} images in a behaving animal is the main objective of this doctorate project. Imaging neuronal structures and monitoring activity within them will allow for a better understanding of how the brain functions. While it is possible to produce images from animals where the head is fixed, or while the animal is under anaesthesia, as both approaches limit movement, neither strategy is as valuable as permitting an awake animal to move and behave freely.

% Subs 3
\subsection{Optics and fibre optics}

% Subsubs 3.1
\subsubsection{Light propagation}

This subsection follows the work presented in Chapter 1 of \cite{fourier_optics}.

Light is electromagnetic radiation. Thus, its evolution can be described with the use of Maxwell’s equations. In the form shown below, an assumption is made that the media contains no free charge, which is generally the case in microscopy. Here $\vec{E}$ and $\vec{H}$ are the electric and magnetic field vector, respectively, $\epsilon$ is the permittivity of the medium and $\mu$ is the permeability of the medium.

\begin{gather}
    \nabla \times \vec{E} = - \mu \frac{\partial \vec{H}}{\partial t} \nonumber \\
    \nabla \times \vec{H} = \mu \frac{\partial \vec{E}}{\partial t} \nonumber \\
    \nabla \cdot \epsilon \vec{E} = 0 \nonumber \\
    \nabla \cdot \mu \vec{H} = 0
\end{gather}

The relations above can then be transformed into the wave equations below, where $n$ is the refractive index of the media, which could vary spatially, and $c$ is the speed of light in vacuum. 

\begin{gather}
    \nabla^2 \vec{E} - \frac{n^2}{c^2} \frac{\partial^2 \vec{E}}{\partial t^2} = 0 \nonumber \\
    \nabla^2 \vec{H} - \frac{n^2}{c^2} \frac{\partial^2 \vec{H}}{\partial t^2} = 0
\end{gather}

In the case of a non-varying refractive index, further simplifying approximations can be made. More specifically, that there is no coupling between vector components of the fields. Thus, the wave equation above can be applied individually to the field's components. For example, below is shown the formulation for the x-component of the electric field.

\begin{equation}
    \nabla^2 E_x - \frac{n^2}{c^2} \frac{\partial^2 E_x}{\partial t^2} = 0
\end{equation}

Thus, in order to avoid deriving similar solutions for every component, a scalar field $u$ can be used instead, as shown below. It is important to note that the assumption is not valid if the refractive index varies. Scalar approximation also can induce errors at the interface between two media, where different boundary conditions might apply to different components. Later in the report, when describing light propagation through a fibre, the analysis will consider the full vectorial approach. However, for beams propagating in free space, these assumptions are practical.

\begin{equation}
    \nabla^2 u - \frac{n^2}{c^2} \frac{\partial^2 u}{\partial t^2} = 0
\end{equation}

Through the use of a Fourier transform in the time domain, the wave equation above can be modified into the Helmholtz equation form shown below, where $k^2=-\frac{\omega^2 n^2}{c^2}$ is the wave number, $\omega=2\pi \nu$ and $\nu$ is the temporal frequency of the light. As in the case of monochromatic light the periodic temporal component is constant across the optical field during propagation, it is usually omitted in derivations.

\begin{equation}
    \nabla^2 u + k^2 u = 0
\end{equation}

% subsubs 3.2
\subsubsection{Diffraction and Fourier optics}

This section summarises the work presented in Chapters 3 and 4 of \cite{fourier_optics}.

A general plane wave can be described as $u(x,y,z) = U_0 \exp[i (k_x x + k_y y + k_z z)]$, which is a solution to the Helmholtz equation, where $(k_x,k_y,k_z)$ are the wavenumber components, which satisfy $|\vec{k}|^2 = k_x^2 + k_y^2 + k_z^2 = \left( \frac{2 \pi}{\lambda} \right)$, and $U_0$ is a complex number describing the amplitude and phase offset of the plane wave. The reason for introducing plane waves is the possibility of decomposing a complex field into a superposition of such plane waves.

In imaging, it is of interest to describe the change in the optical field from one plane to another. Mathematically, we are interested in the transformation $u(x,y,z=0) \rightarrow u(x,y,z=z)$, where the starting position can be considered at $z=0$ without loss of generality. 

The method in which a field can be decomposed into planar waves is through the use of Fourier transform. Applying the Fourier transform to the field at $z=0$, result in the equation below, where $A(k_x,k_y,z=0)$ is the transverse angular frequency spectrum of the field.

\begin{equation}
    A(k_x, k_y, z=0) = \iint_{-\infty}^{\infty} u(x,y,z=0) e^{-i (x k_x + y k_y)} dx dy
\end{equation}

If the inverse Fourier transform is applied, as shown below, the equation shows that the field $u(x,y,z=0)$ is a superposition of plane waves.

Due to the restriction $|\vec{k}|^2 = k_x^2 + k_y^2 + k_z^2 = \left( \frac{2 \pi}{\lambda} \right)$ being applied, any plane wave defined by $(k_x,k_y)$ has a $k_z$ associated. Thus, each plane wave will have a direction of propagation associated. 

\begin{equation}
    u(x, y, z=0) = \iint_{-\infty}^{\infty} A(k_x,k_y,z=0) e^{i (x k_x + y k_y)} dk_x dk_y
\end{equation}

Repeating the Fourier transforms for the field at position $z$ and applying the Helmholtz equation leads to the result below. This is a multiplication with a transfer function in the angular domain.

\begin{equation}
    A(k_X, k_y, z=z) = A(k_x, k_y, z=0) \exp \left[ iz \sqrt{\frac{4 \pi^2}{\lambda^2} - k_x^2 - k_y^2}  \right]
\end{equation}

The term under the square root is not always positive. More specifically, it is negative for $k_x^2+k_y^2>\frac{4 \pi^2}{\lambda^2}$, which introduces another $i$ term into the exponential, thus transforming it from an oscillation to decaying function. In other words, plane wave components corresponding to high frequency values become an evanescent wave, which does not propagate. For this reason, there is a diffraction limit to the achievable resolution in an imaging system.

The Fresnel approximation can be made if the distance between planes is comparably larger than the wavelength, such that $\lambda^2 k_x^2,\; \lambda^2 k_y^2 << 4 \pi^2$, as the angles of interest of the component plane waves are small. Thus, a first order approximation using the Taylor expansion of the square root can be done: $\sqrt{1 - \frac{\lambda^2 k_x^2}{4 \pi^2} - \frac{\lambda^2 k_y^2}{4 \pi^2}} \approx 1 - \frac{\lambda^2}{8 \pi^2} (k_x^2 + k_y^2)$. 

Thus, the transfer function exponent becomes:

\begin{equation}
    H(f_x, f_y; z) = \exp \left[ i2 \pi z \frac{1}{\lambda} \left[ 1 - \frac{\lambda^2}{8 \pi^2} (k_x^2 + k_y^2) \right] \right] 
    = \exp{\frac{i2 \pi z}{\lambda} \exp \left[ -i \frac{z \lambda}{4 \pi} (k_X^2 + k_y^2) \right]}
\end{equation}

The first term is a constant global phase which can be ignored. Applying the inverse Fourier transform to the transform function and convolving the result with the $u(x,y,z=0)$ field leads to Fresnel’s diffraction integral.

One application of this theory is for lenses. The change in field between entering and exiting a lens can be generally described though the use of a multiplication, as shown below, where $t_{lens} (x,y)$ is the modulation of the field done by the lens.

\begin{equation}
    u_{exit} (x,y) = t_{lens} (x,y) u_{entry} (x,y)
\end{equation}

In the paraxial approximation, the modulation is given by the form below, where $f$ is the focal length of the lens. This modulation can also include a pupil function which blocks light passing outside the lens facet.

\begin{equation}
    t_{lens} (x,y) = \exp \left[ - \frac{ik}{2f} \frac{x^2 + y^2}{2} \right]
\end{equation}

The application of interest is the transformation of the field from a focal plane to the other. The form of such transformation in the angular spectrum domain would be in the form shown below, where the H transform is applied twice for the diffraction to and from the focal planes and a convolution at the lens.

\begin{equation}
    A(k_x, k_y, z=+f) = H(k_x, k_y; f) \left[ T_{lens}(k_x, k_y) \otimes \left[ H(k_x, k_y; f)A(k_x, k_y, z=-f) \right] \right]
\end{equation}

The same relationship can be written in the spatial domain by using the inverse Fourier transforms versions of the functions and swapping multiplication with convolution.

The final result in the spatial domain is shown below, where (x,y) are the transversal spatial coordinates of the input plane and (v,w) are the transversal spatial coordinates at the output plane.

\begin{equation}
    u(v,w,z=+f) = \frac{1}{i \lambda f} \iint_{-\infty}^{\infty} u(x,y,z=-f) \exp \left[ -i \frac{2\pi}{\lambda f} (xv+yw) \right] dx dy
\end{equation}

This result is a Fourier Transform, where the output spatial coordinates are proportional to the spatial frequencies of the input field. More specifically: $k_x=\frac{2 \pi v}{\lambda f}$, $k_y = \frac{2 \pi w}{\lambda f}$. This allows for access to the Fourier plane of an optical field by passing it through a lens.

% subsubs 3.3
\subsubsection{Optical fibres} \label{sss:backg_optical_fibres}

When describing the transmission through an optical fibre, due to the circular symmetry of an ideal fibre, a cylindrical coordinate system is best suited. Thus, the parameters used are $r$ (radius coordinate), $\phi$ (azimuthal angle coordinate), $z$ (axial coordinate) and $t$ (time coordinate). Limiting the analysis to only one specific value of the optical field in z-component of the spacial frequency $k_0$ and one value in the temporal frequency $\omega$ (eliminating the Fourier transform integrals), the electric field can be described as below, where the propagation constant $\beta=n k_0$ is used instead. The same relationship can be written for the magnetic field $\vec{H}$.

\begin{equation}
    \vec{E}(r, \phi, z, t) = \vec{E}(r, \phi) \exp [ i ( \omega t - \beta z) ]
\end{equation}

As the unit vectors in the radial and azimuthal direction are not constant, solving the Maxwell equations is hard. However, this is not true for the $z$-component, where the usual form of the resultant Helmholtz equation $(\nabla^2 + k^2) E_z = 0$ still holds, where $\nabla^2$ is the Laplacian operator in cylindrical coordinates and $k$ is the wavenumber. The same relationship applies to the magnetic field component $H_z$. Solving the partial differential equation and applying the boundary conditions will lead to a set of propagating fibre modes. The details described above and the subsequent derivation for a step-index optical fibre is described in Chapter 3 of Amnon Yariv’s book \cite{backg_optical_elec}. However, these modes, while they describe the propagation in an optical fibre, are difficult to work with, especially in an experimental setting, because the polarisation of light is not uniform across the transverse plane $(r-\phi)$. 

Working under the assumption that the difference between the core and cladding refractive indices is small (weakly guiding fibres), the boundary continuity conditions can be simplified and the propagation modes can be approximated to a set of linearly polarised modes (LP modes), which have the property of uniform polarisation. For the x-polarisation state, the transverse electric field takes the form: $E_x = A J_l (hr) \exp{il\phi} \exp [ i (\omega t - kz)]$, where $A$ is a constant given by boundary conditions, $J_l$ is the Bessel equation of first kind of order $l$, $l$ is azimuthal index taking integer values $(…,-2,-1,0,1,2,…)$, $h^2=k^2-\beta^2$. The magnetic transverse field takes the form: $H_y=\frac{\beta}{\omega \mu} E_x$, where $\mu$ is the permeability of the fibre core. As the light is linearly polarised, $E_y=0$,$H_x \approx 0$. The y-polarisation state analysis is identical, as the azimuthal symmetry of the fibre can be used to approach the problem with a rotated frame of reference, leading to the same results.

One situation which arises is that the x-polarisation and y-polarisation modes are degenerate, in other words, their propagation constant $\beta$ is equal. This results in coupling between the two modes even in an ideal straight fibre. However, as the transversal field envelope is identical, coupling between the two modes only changes the direction of the linear polarisation. For a short length of the fibre, the coupling between polarisation states is small, thus, the light can be assumed to maintain its linear polarisation and the scalar approximation becomes a reasonable model, by only focusing on the component of the electrical field on which the polarisation is aligned with. Mathematically, instead of working with the vector $(E_x, E_y = 0)$, $E_x$ can be used alone.

Polarisation coupling can present an issue experimentally, for example when measuring the phase of a field using interferometry. However, while keeping the property of uniform polarisation, circularly polarised (CP) modes can be used instead. Plöschner et al. \cite{mmf_neuron} have shown that CP modes have the same profile as the LP modes, however, due to spin-orbit interaction between the polarisation and orbital angular momentum generated by the phase profile, the propagation constants of the modes with different polarisation are different. Thus, the modes are not degenerate anymore. This leads to circular polarisation being preserved. The specific propagation constants can be approximated to the first order from the LP propagation constants. The project has used CP modes.

The general case described above is of the multimode fibre (MMF). A specific case which arises, is the single-mode fibre (SMF), where there is only one propagating fibre mode. The number of modes in a fibre is determined by specifications such as core diameter, numerical aperture (NA), wavelength of the light propagated through the fibre and its refractive index profile. In terms of the latter, different versions can be constructed: step-index fibre, where the core and cladding have different but uniform refractive indices, or graded-index fibre where the refractive index of the core gradually decreases from a maximum value to the cladding value \cite{backg_review_optical_fibres}.

The advantage of SMFs, is that there is no modal scrambling, as there is only one propagation mode. This can be important in communications applications, as it allows for a larger bandwidth of data transmission. In imaging applications, where controlling the position of the output focus is required, using one SMF is not useful. Multiple SMFs can be used in endoscopic applications in the form of a bundle of fibres, where each fibre corresponds to a pixel in the final image. However, the large size of the fibre bundle defeats the purpose of minimising invasiveness in neurological applications. Thus, a MMF is a better choice. However, a SMF can still be used for example to carry the reference beam required for the interferometry used in calibrating the MMF.

In terms of refractive index profile, it was shown that graded-index fibres are more robust to changes in transmission due to bending \cite{backg_graded_index_fibres} and that they are a good candidate to use in imaging. However, for graded-index fibres to be as resistant to changes in propagation due to bending, a precise refractive index profile is required. Unfortunately, current manufacturing precision is still not at a desired level. Thus, commercially available fibres still require correction to bending to be used in imaging. This project has only focused on the step-index fibre correction. However, the correction procedure is applicable to graded-index fibres, the only requirement being a correct description of the graded-index fibre modes.

% subsubs 3.4
\subsubsection{Calibrating an optical fibre} \label{ssec:fibre_calibration}

Input light to a MMF couples into the different fibre modes and in the ideal case is transmitted without coupling to the output. However, due to the propagation constants of the modes being different, the phases of the output components are desynchronised, which leads to observation of a seemingly random speckle pattern. However, this scrambling is not stochastic, but deterministic. The transmission can thus be measured empirically in an experimental setting.

Such measurement is usually formulated as a matrix, coined the transmission matrix. In order for a matrix to be used, the optical fields at the input and output of the fibre need to be composed of a finite number of sampled locations. Then, the transmission matrix describes the transmission of light from every input sample location to every output sample location. Mathematically, it can be described by a matrix multiplication: $\mathbf{y}=\mathbf{Tx}$, where $\mathbf{x}$ and $\mathbf{y}$ are the input and output fields respectively, each represented as a vector that maps onto the a Cartesian space of pixels, and $\mathbf{T}$ is the transmission matrix. An illustration of the process undergone in the application of the TM to relate the input to the output of a fibre is shown in figure \ref{fig:apply_TM_diagram}.

\begin{figure}
    \centering
    \includegraphics[width=\textwidth]{images/TM_diagram.png}
    \caption{Diagram illustrating the procedure to apply a given TM to relate an input field to the output of the fibre.}
    \label{fig:apply_TM_diagram}
\end{figure}

The matrix approach is similar to describing the propagation of light through a lens. Between the input and the output focal planes there is a Fourier transform relationship, which can be described through the use of a discrete Fourier transform matrix. Thus, changing the objective lens of a microscope system for a MMF leads to a similar approach to describe light transmission, through the use of a matrix multiplication.

A description of the calibrating algorithm is given in \cite{backg_calibration}. The process works by focusing light at the input facet into one sample location at a time. At the output end of the fibre, a camera captures the pattern observed. In order for phase information to be extracted, the output facet is interfered with a reference beam. Then the global phase of the input light is changed, by adding a constant phase across the spatial light modulator (SLM), scanning in the range of $[0,2\pi)$. As the reference beam has a constant phase and the output pattern only has its global phase changed, the relative phase information between output pixels can be retrieved by observing the interference pattern. One iteration of this process relates one input to all outputs. Thus, a column of the TM is measured. The process is then repeated for all the other input pixels, until the entire TM is measured.

In order to generate a focus at the output facet, the values of the input pixels is chosen such that their complex summation of their individual contribution to the output pixel maximises intensity (absolute value). This can be achieved by setting the input phases such that their output contributions are in phase at the output pixel location.
One important note to mention is that the algorithm does not require the output plane to be at the fibre facet. The acquisition camera can be set such that it focuses away from the facet. Calibrating the fibre in such a manner allows for the generation of a TM which relates the input field to an output field at the focal plane of the acquisition camera.

% subsubs 3.5
\subsubsection{Matrix representation of the transmission through an optical fibre}

For a given fibre mode of index $k$, let its profile be described as a column vector $\mathbf{\psi}_k$, where each element represents the optical field value sampled at a given location. The locations are taken to be a Cartesian mesh-grid to better overlap with what a camera might observe or how an SLM might modulate the light field. Thus, for a specific mode coefficient $A_k$, the optical field in the Cartesian pixel basis can be described as: $\mathbf{x_C} = \mathbf{\psi_k} A_k$.

The light propagating through an optical fibre, is a superposition of propagation modes, each with an associated coefficient. Thus, the superposition equation $\mathbf{x_C} = \sum_{k=0}^M \mathbf{\psi_k} A_k$, where $M$ is the number of propagation modes, can be represented as a matrix multiplication, where $\mathbf{\Psi}$ is the mode matrix and $\mathbf{x}_M$ is the column vector of mode coefficients, or in other words, the optical field in the fibre mode basis:

\begin{equation} \label{eq:psi_mode_to_cart}
    \mathbf{x}_C = 
    \begin{pmatrix}
        \mathbf{\psi_1} & \mathbf{\psi_2} & \dots & \mathbf{\psi_M} 
    \end{pmatrix}
    \begin{pmatrix}
        A_1 \\ A_2 \\ \vdots \\ A_M
    \end{pmatrix}
    = \mathbf{\Psi} \mathbf{x}_M
\end{equation}

Thus, $\mathbf{\Psi}$ represents the matrix transformation from the modal basis to the Cartesian basis. As the fibre modes form an orthogonal basis, then a normalisation can be chosen for the fibre mode profiles such that $\mathbf{\psi}_k^{\dagger} \mathbf{\psi}_k = 1$ and $\mathbf{\psi}_k^{\dagger} \mathbf{\psi}_j = 0$ for $k \neq j$, where $\dagger$ represents the conjugate transpose. Thus leading to the relationship:

\begin{equation}
    \mathbf{\Psi}^\dagger \mathbf{\Psi} = \mathbf{I}_M
\end{equation}

Multiplying equation \ref{eq:psi_mode_to_cart} by $\mathbf{\Psi}^\dagger$ on the left side leads to:

\begin{equation}
    \mathbf{\Psi}^\dagger \mathbf{x}_C = \mathbf{\Psi}^\dagger \mathbf{\Psi} \mathbf{x}_M = \mathbf{x}_M 
\end{equation}

Thus, $\mathbf{\Psi}^\dagger$ represents the matrix transformation from the Cartesian basis to the modal basis. It is important to note that $\mathbf{\Psi} \mathbf{\Psi}^\dagger \neq \mathbf{I}_M$. More specifically, this matrix transformation represents the projection of an arbitrary optical field into the modal basis. Thus, components of the optical field which do not couple into the fibre modes are lost. Specifically, it is the transformation from Cartesian basis to modal basis which generates this projection ($\mathbf{\Psi}^\dagger$). However, if the Cartesian optical field is already entirely within the mode basis space, then the matrix transformation has the same effect as an identity matrix (no transformation).

An ideal fibre will generate no coupling between the modes as light propagates through it and the propagation is defined by the propagation constants of the modes only. This relationship can be described in terms of a transmission matrix $\mathbf{TM}_M$, where $\mathbf{x}_M$ and $\mathbf{y}_M$ represent the input and output optical fields in the mode basis, $\beta_k$ represents the propagation constant of mode $k$ and $L$ represents the length of the fibre:

\begin{equation}
    \mathbf{y}_M = 
    \begin{pmatrix}
        e^{i\beta_1L}  & 0  & \dots & 0 \\
        0 & e^{i\beta_2 L} & \dots & 0 \\
        \vdots & \vdots & \ddots & \vdots \\
        0 & 0 & \dots & e^{i\beta_M L}
    \end{pmatrix} \mathbf{x}_M = \mathbf{TM}_M \mathbf{x}_M
\end{equation}

Making use of the Cartesian to mode basis transformation and considering that the output optical field is assured to be entirely within the mode basis space, thus, the projection to mode basis is acting as an identity matrix, a relationship can be described between the Cartesian based TM and the modal based TM.

\begin{equation}
    \mathbf{\Psi}^\dagger \mathbf{y}_C = \mathbf{TM}_M \mathbf{\Psi}^\dagger \mathbf{x}_C \nonumber
\end{equation}

\begin{equation}
    \mathbf{y}_C = \mathbf{\Psi} \mathbf{TM}_M \mathbf{\Psi}^\dagger \mathbf{x}_C
\end{equation}

\begin{equation}
    \mathbf{TM}_C = \mathbf{\Psi} \mathbf{TM}_M \mathbf{\Psi}^\dagger
\end{equation}

Multiplying by $\mathbf{\Psi}^\dagger$ on the left side and with $\mathbf{\Psi}$ on the right side and using the orthogonality of the modes gives the mode base TM as a function of Cartesian base TM:

\begin{equation}
    \mathbf{TM}_M = \mathbf{\Psi}^\dagger \mathbf{TM}_C \mathbf{\Psi}
\end{equation}

% subsubs 3.6 
\subsubsection{Analytical description of the transmission through a bent multimode fibre}

Once the propagation through a MMF is formulated, a closer look at the bending effects can be taken. M. Plöschner et al. \cite{seeing_through_chaos} stated that there are two main effects which take place: geometric and compressive. The geometric bending effect occurs due to a change in the optical path length through the fibre, as light propagating on the inside of the bending travels a shorter path than on the outside. The compressive bending effect occurs due to stresses in the fibre which change its refractive index, thus affecting the optical path length as well. These effects were formulated theoretically for the case of a circular fibre bend of a given curvature radius. 

The specific analysis can be found in the methods section of the cited article, however, the main conclusion is that the bent fibre will have a different set of propagation modes, with a different set of propagation constants associated. While these could be computed, it would need to be done for all bending radii, individually. A more straightforward approach is to keep the straight fibre mode basis, but acknowledge that due to the change in actual propagation mode basis, the transmission is no longer without mode coupling.

The transmission through the circularly bent fibre in the x-z plane can be formulated though the use of a helper matrix $\mathbf{B}$ whose elements are given below, where $\beta_k$ is the straight fibre propagation constant of mode $k$, $\delta_{kj}=1$ for $j=k$ and $0$ otherwise, $n_core$ is the core refractive index, $k_0$ is the vacuum wavenumber, $\rho$ is the bending radius, $\psi_k$ and $\psi_j$ are the mode profiles of modes $k$ and $j$, respectively. $\xi = 1 - (1-2\sigma) \frac{n_{core}-1}{n_{core}}$ is a factor which arises due to the compressive bending effect, where $\sigma$ is the Poisson ratio of the core material.

\begin{equation}
    B_{kj} = \beta_k \delta_{kj} - \frac{n_{core} k_0}{\rho} \langle{\psi_k} | x \xi | \psi_j \rangle  
\end{equation}

The bra-ket term can be expanded into:

\begin{equation}
    \langle \psi_k | x \xi | \psi_j \rangle = \frac{\iint_{-\infty}^{\infty} \psi_k^* x \xi \psi_j dx dy}{N_k N_j}
\end{equation}

\noindent where:

\begin{equation}
    N_k = \sqrt{\iint_{-\infty}^{\infty} |\psi_k|^2 dx dy }
\end{equation}

Finally, the mode base transmission matrix on the bent fibre is given by the matrix exponential:

\begin{equation}
    \mathbf{TM}_{M,bent} = e^{i\mathbf{B}L}
\end{equation}

An arbitrary shape could then be approximated as a series of small circularly bent sections, thus allowing for the calculation of the transmission through a fibre with a specific shape.

In the case where the bending is not done in the x-z plane, but at an angle $\theta$ counter-clockwise with respect to x-axis, the same approach can be used, but substituting $x$ in the bra-ket term with $x \cos{\theta} + y \sin{\theta}$ (equivalent to an axis rotation).

% subsubs 3.7
\subsubsection{Correcting for bent optical fibres}

S. Farahi et al. \cite{backg_dynamic_bending_comp} have implemented a method for measuring the conformation of a fibre by delivering light at the distal end of the fibre, the end further from the laser source, through the use of micro-optics. The output at the proximal end, towards the laser source, was measured and compared against a large databank of images. The drawbacks of this method are the large dataset required before imaging and the need for access to the distal end. The latter would be impossible in an in vivo biological imaging context as it would be implanted inside the mouse brain.

Ruo Yu Gu et al. \cite{backg_design_flexible_endoscope} have proposed a method to characterise the transmission matrix in a bent fibre with the use of a checkerboard reflector at the distal end of the fibre, which eliminates the need for access at that end. Starting with an initial calibration, a new transmission matrix is calculated continuously to measure the change. The drawback is that the new calibration needs to be fast in order to account for the continuous change in fibre shape. This can prove difficult for a large number of modes, as the matrix itself is large. Noise can also prove to be an obstacle, as increasing the number of modes leads to a requirement for better signal to noise ratio.

G. Gordon et al. \cite{backg_metasurface_approx_bent_TM} have described a method which builds on the previous one. It relies on stacking multiple checkerboard reflectors, which provide further encoding based on polarisation and wavelength. This is a more general approach as it would a apply to non-unitary transmission matrices, which would occur in cases where energy is lost during transmission. The matrix is inferred by using an iterative algorithm based on a first order model approximation from recording done at different wavelengths. The main drawback is the time requirement, as in their case 14.5 minutes were needed for convergence. This is extremely slow if the technology target is to strive for correcting rapid changes in fibre shape, which might arise during imaging of more active specimens.

Given that the main limitation in the above-mentioned two methods is the time requirement for correction during imaging, a better approach could be through the use of machine learning. In this way, computational time is shortened at the expense of lengthy training, which can be done prior to imaging.

Y. Liu et al. \cite{backg_bending_deep_learning} have trained a convolutional neural network on speckle patterns recorded at the distal end of the fibre, which was deformed by translating the proximal end in the z-dimension (the central axis of the fibre). The data was acquired at 21 different proximal end positions, resulting in 21 classes for which a classifier network was trained. The data was split into a training dataset and a validation dataset. The training result was quantified through measuring the accuracy of inference on both the training dataset and validation dataset, the accuracy being defined as the percentage of correct classifications. While this method uses the distal end, it did show that the accuracy of a trained network can be over $90\%$ for a specific fibre. The authors trained two networks for two different fibres, but did not investigate the network accuracy on data obtained from a different fibre compared to the training dataset. 

P. Fan et al. \cite{backg_multimode_fibres_deep_learning} have trained a network on recorded images at the distal end, while inputting at the proximal end a set of binary images using digits or letters. There were three different bending states used, each with its own trained network. The classifier’s accuracy was over $90\%$, but lacked when applied to a state different to its trained one. This shows that there is a low correlation coefficient of the output between the states. This reduces the possibility of using a simple linear regression to determine the correct change in transmission, while at the same time it could increase the chances of a more complex network to differentiate between states, which is my goal.

In order to train a network to correctly calculate the change in transmission, a large amount of data is needed. This could prove difficult and time consuming to do experimentally. Thus, through the use of a simulation, training data can be generated. While the project will not focus on the use of machine learning, the simulation program can be used for a network approach.

% subsubs 3.8
\subsubsection{Experimental optical system} \label{sssec:optical_system}

The project, while heavily focused on theoretical and computational approaches, was meant to have a experimental application on an existing optical system capable of imaging through a MMF. A simplified schematic of the microscope system is illustrated in figure \ref{fig:optical_system_diagram}. The system can be split into three modules: the laser module, the MMF module and the calibration module. The purpose of the laser module was to split the light from the laser and send it to the other two stages. The laser (CrystalLaser DL488-030-SO) generated a continuous waveform beam of linearly polarised light with a wavelength of 488 nm. The halfwave plate (HWP1) changed the direction of the linear polarisation, such that the following polarising beam splitter (PBS) was able to send two orthogonal linearly polarised beams to the optical fibres for light transport.

Light to the MMF module was sent through a polarisation maintaining optical fibre (PMF), as the light reaching the SLM needed to be aligned with the polarisation axis of the SLM. An iris was placed in the Fourier plane of the SLM in order to filter only one diffraction order of the light reflected by the modulator. The filtered light then passed through a quarter-wave plate (QWP1) which changed the polarisation of light from linear to circular. The reason for using circularly polarised light through the MMF was due to circular polarisation being better preserved than linear polarisation 
(section \ref{sss:backg_optical_fibres}). The filtered light was then passed through a 4f system to properly couple the light into the proximal facet of the MMF.

\begin{figure}
    \centering
    \includegraphics[width=\textwidth]{images/simplified_optics_diagram.png}
    \caption{Simplified schematic of the optical system target for which the computational simulations were used. Abbreviations are defined in the main text.}
    \label{fig:optical_system_diagram}
\end{figure}

During calibration, a calibration module was placed at the output of the MMF. Light from the MMF was magnified by an objective lens (OL) after which it passed through another quarter-wave plate (QWP2) which changed the polarisation from circular to linear, which then was interfered with a reference beam onto a non-polarising beam splitter (NPBS). The reference beam was delivered from the laser module by a single-mode fibre (SMF). While the single-mode fibre maintained the linear polarisation, it did not maintain its direction. Thus, before the reference beam was interfered on the NPBS, a half-wave plate (HWP2) was used to change the polarisation direction to match the light received from the MMF. The resultant interference pattern was imaged onto a charged-coupled device (CCD) camera.
During imaging, the calibration module was replaced by the animal to be observed, by implanting the MMF in the brain tissue. At the output plane of the MMF, laser light was focused onto the tissue. Emission of fluorescent light was captured by the MMF and passed back through the 4f system. A dichroic mirror (DM), which only reflected light at 488 nm was used to filter the emitted light, which was at a different wavelength. The emitted light was captured and measured by a photo-multiplier tube (PMT). Each measurement of the PMT represented the intensity of one pixel in the final image. 

% subsub 3.9
\subsubsection{Simulation environment}
The thesis is based on a simulation environment created in Matlab, which allows for the numerical simulation of light propagating within a multimode fibre.

The first step in the program development has been to generate the linear polarisation modes in an optical fibre.  For this step the function present in the GitHub repository created by Szu-Yu Lee have been used \cite{LP_generation_code}. However, this is the only code in the simulation program which has been taken from an external source.

Then, using the mathematical approaches in article \cite{seeing_through_chaos}, functions were written to compute the changes in propagation constants and profile of the CP modes, allowing for future simulations to be done in a circular polarisation basis. Furthermore, using information from the same article, the functions computing the mode coupling for a circularly bent fibre were also written.

These functions represent the basis of the simulation program. Further additions are described in future chapters.

\end{document}