\providecommand{\main}{..}

\documentclass[main.tex]{subfiles}
\graphicspath{ {../images/} }
\begin{document}
% Bragg section
\section{Bragg gratings as point reflectors in multi-mode fibres}

%%%%
\subsection{Introduction}
One of the issues with using multiple reflectors is that the reflected light from each of them needs to be differentiated, in order to achieve a focus at a specific reflector. Using either small mirrors, point scatterers or fluorescent beads could lead to the same amount of reflected or emitted light whether the reflection or emission occurs at a specific location or from multiple locations at the same time. In the case of a reflector that maintains coherence, the phase information is preserved. Thus, phase information can not be used for differentiation as the phase pattern reflected is dependent on bending. As the bending is the scope of the correction, it is unknown and its effect on the change in phase cannot be known.

Bragg gratings can provide the solution to the problem at hand. The reflected light from a Bragg grating can be filtered based on its wavelength, as a function of the grating geometry \cite{bragg_book}. This would allow for differentiation by simply using a wavelength based filter to differentiate between light reflected from each grating.

However, two important points will need to be considered. Firstly, the memory effect approach requires focusing at a specific point, which requires that the grating act as a point reflector. Hence, the Bragg grating must be below a characteristic size determined by the physics of the light evolution in the fibre. This could have negative effects experimentally, as the amount of light reflected can also be small. Secondly, being able to achieve a focus at the grating allows for the generation of an assumed TM which is specific to the wavelength of light reflected. Thus, the change in TM due to wavelength needs to be considered before generating the input field for generating an excitation focus, if the focus is to be created using a different wavelength to that for which the Bragg grating was designed.

This chapter focuses on describing a more general approach in simulating a Bragg grating inside a multi-mode fibre. It starts from existing theoretical background and improves the approach to consider coupling between multiple modes, both co-propagating and counter-propagating. The aim is to reduce the problem to a system of de-coupled linear ordinary differential equations, finally, representing the Bragg grating effect on the optical field as a scattering matrix, relating the light exiting the grating section of the fibre to the light entering the grating section of the fibre.

Finally, using this approach, it was investigated whether A Bragg grating can be used as a reflector (localised reflection) and what is the bandwidth of such a grating.

%%%%
\subsection{Existing theoretical background} \label{ssec:bragg_theory_old}

The theory presented in this subsection follows the coupled-modes approach undertaken in chapter 4 of \cite{bragg_book}: ''Theory of Fiber Bragg Gratings" with the aim to describe the coupling between modes generated by the Bragg grating.

The approach starts from the Helmholtz equation for the electric field, where $\mathbf{P}$ is the polarisation field.
\begin{equation}
    \nabla^2 \mathbf{E} = \mu_0 \varepsilon_0 \frac{\partial^2 \mathbf{E}}{\partial t^2} + \mu_0 \frac{\partial^2 \mathbf{P}}{\partial t^2}
\end{equation}

The grating can be assumed as an introduced perturbation to the polarisation field, where `$up$` stands for unperturbed and `$p$` for perturbed.

\begin{equation} \label{eq:P_total}
    \mathbf{P} = \mathbf{P}_{up} + \mathbf{P}_p
\end{equation}

The unperturbed polarisation field represents the relative permittivity of the unperturbed material, which is directly related to the refractive index of the material.

\begin{equation}
    \mathbf{P}_{up} = \varepsilon_0 (\varepsilon_r-1) \mathbf{E} 
\end{equation}

This leads to the modified Helmholtz equation:
\begin{equation} \label{eq:mod_Helm}
    \nabla^2 \mathbf{E} = \mu_0 \varepsilon_0 \varepsilon_r \frac{\partial^2 \mathbf{E}}{\partial t^2} + \mu_0 \frac{\partial^2 \mathbf{P}_p}{\partial t^2}
\end{equation}

If the perturbation is an isotropic change in refractive index, thus, not introducing any coupling between orthogonal polarisations, then the vector notation can be dropped and the problem can be treated as a scalar field, by looking at only one specific light polarisation. 

The total electric field can be represented as a superposition of forward propagating modes (FP), noted with positive index, and backward propagating modes (BP), noted with negative index, where $A$ represents the mode coefficient, $\xi_h$ represents the profile of mode $h$ (mode $-h$ has the same profile, but flipped to overlap transversally with mode $h$), $\omega$ is the angular velocity of the temporal oscillation, $\beta$ is the absolute value of the propagation constant (identical for modes $h$ and $-h$), $N$ is the number of modes and $cc$ represents the complex conjugate terms (i.e. $\cos(x) = \frac{1}{2} ( e^{ix} + cc )$ ).

\begin{equation} \label{eq:total_E}
    E = \frac{1}{2} \left[ \sum_{h=1}^{N} A_h \xi_h e^{i (\omega t - \beta_h z)} + \sum_{-h=-1}^{-N} A_{-h} \xi_h e^{i (\omega t + \beta_h z)} \right] + cc
\end{equation}

A brief discussion on the modes $h$ and $-h$ is required. As mentioned above, in order for the modes $h$ and $-h$ to couple into each other, they must have the same profile transversally, however, as they propagate in opposite directions, their helicity is opposite. In other words, the modes have the same $m$ index but opposite $l$ index. The same observation can be made about polarisation as well. For the coupling to occur, the transversal polarisation must match, but as the modes propagate in opposite directions, the circular polarisation of the two modes must be orthogonal (one is left-hand while the other is right-hand and vice versa). However, as both the profile helicity and the polarisation is flipped between the modes, the propagation constants of the two modes have the same absolute value, but opposite signs.

Substituting equation \ref{eq:total_E} for $E$ in \ref{eq:mod_Helm}, expanding the Laplacian and assuming a slowly varying envelope approximation (equation \ref{eq:SVEA}) leads to equation \ref{eq:helm_expanded}.

\begin{equation} \label{eq:SVEA}
    \frac{\partial^2 A_h}{\partial z^2} \ll \beta_h \frac{\partial A_h}{\partial z} 
\end{equation}

\begin{equation} \label{eq:helm_expanded}
    \sum_{h=1}^{N} -i \beta_h \frac{\partial A_h}{\partial z} \xi_h e^{i (\omega t - \beta_h z)} + \sum_{-h=-1}^{-N} i \beta_h \frac{\partial A_{-h}}{\partial z} \xi_h e^{i (\omega t + \beta_h z)} + cc = \mu_0 \frac{\partial^2 P_p}{\partial t ^ 2}
\end{equation}

The next step is to use the orthogonality relationship of the modes (equation \ref{eq:mode_orth}) by multiplying both sides by $\xi_m^\ast$ and integrating over the transversal plane (over $x$ and $y$), which leads to equation \ref{eq:mode_coupling_base}.

\begin{equation} \label{eq:mode_orth}
    \iint_{-\infty}^{\infty} \xi_a \xi_b^\ast dx dy = \frac{2 \omega \mu_0}{\beta_a} \delta_{ab}
\end{equation}

\begin{equation} \label{eq:mode_coupling_base}
    \frac{\partial A_m}{\partial z} e^{i (\omega t - \beta_h z)} - \frac{\partial A_{-m}}{\partial z} e^{i (\omega t + \beta_h z)} + cc = \frac{i}{2 \omega} \iint_{-\infty}^{\infty} \frac{\partial^2 P_p}{\partial t ^ 2} \xi_m^\ast dxdy
\end{equation}

\noindent This new form has allowed for differential terms to be isolated rather than in a sum. 

The next step is to look into the form of the perturbed polarisation field. Writing the perturbed relative permittivity $\varepsilon_{r,total} = \varepsilon_r + \Delta \varepsilon$,  using the relationship $P = \varepsilon_0 (\varepsilon_{r,total}-1) E$ and equation \ref{eq:P_total} leads to:

\begin{equation} \label{eq:polarisation_perturbation}
    P_p = \varepsilon_0  \Delta \varepsilon E
\end{equation}

\noindent Using the general relationship $\varepsilon_r = n^2$ and writing the perturbation to the refractive index as $n_{total} = n_0 + \delta n$ leads to:
\begin{equation}
    \epsilon_r + \Delta \varepsilon = (n_0 + \delta n)^2 \approx n_0^2 + 2 n_0 \delta n \nonumber
\end{equation}
\begin{equation} \label{eq:var_epsilon}
    \Delta \varepsilon \approx 2 n_0 \delta n
\end{equation}

The grating is assumed to be a cosine variation in $z$ such that the change in refractive index is:
\begin{equation} \label{eq:delta_N_original}
    \delta n = \overline{\Delta n} + \frac{\Delta n}{2} \left[ e^{i \frac{2 \pi N}{\Lambda} z} + cc \right]
\end{equation}

\noindent where:
\begin{itemize}
    \item $\overline{\Delta n}$ is the mean refractive index change
    \item $\Delta n$ is the amplitude of the refractive index change
    \item $\Lambda$ is the grating period
    \item $N$ is the harmonic order
    \begin{enumerate}
        \item If the grating is sinusoidal, then $N$ can be assumed to be 1.
        \item However, for a non-sinusoidal, but periodic, grating, the refractive index change can be described in terms of a Fourier series expansion, which will have harmonics. Though, a later discussion will show why only one harmonic need be considered. 
    \end{enumerate}
\end{itemize}

Combining equations \ref{eq:mode_coupling_base}, \ref{eq:polarisation_perturbation}, \ref{eq:var_epsilon} and \ref{eq:delta_N_original} leads to the equation below:

\begin{equation} \label{eq:mode_coupling_base_expanded}
    \frac{\partial A_m}{\partial z} e^{i (\omega t - \beta_h z)} - \frac{\partial A_{-m}}{\partial z} e^{i (\omega t + \beta_h z)} + cc = \frac{i n_0 \varepsilon_0}{\omega} \iint_{-\infty}^{\infty} \delta n \frac{\partial^2 E}{\partial t^2} \xi_m^* dx dy
\end{equation}

The $\delta n$ term is composed of a sum between $\overline{\Delta n}$ (the DC term) and a term proportional to $\Delta n$ (the AC term). Thus, the right-hand side of equation \ref{eq:mode_coupling_base_expanded} can be written as a sum of two integrals, representing the DC and AC terms, respectively.

The book \cite{bragg_book} continues the derivation by looking into the coupling between a single FP mode and a single BP mode, deriving the more compact result below. Coupling between modes is heavily dependent on the propagation constants of the modes in question and the grating period. More specifically, in the case of co-propagating modes, the coupling is higher if the propagation constants are closer in value, while in the case of counter-propagating modes, the coupling is higher if the difference in propagation constants is close to a multiple of $\frac{2 \pi}{\Lambda}$. A more detailed discussion about detunement is done in sub-section \ref{ssec:3}. Thus, equation \ref{eq:mode_coupling_base_expanded} can be separated into two equations, one for the mode $m$ and another for the mode $-m$ as their propagation constants are very different in value, leading to their coefficient change in $z$ being dependent on different coupling terms. The main assumption made to arrive at the result below is that a specific mode only couples to itself in terms of co-propagation due to DC change in refractive index ($\overline{\Delta n}$) and to only a single counter-propagating mode due to the AC change in refractive index ($\Delta n$).

\begin{equation}
    \begin{cases}
        \frac{\partial A_f}{\partial z} &= -i \kappa_{dc,f} A_f - i \kappa_{ac}^\ast A_{-b} e^{i\Delta\beta z} \\
        \frac{\partial A_{-b}}{\partial z} &= i \kappa_{dc,-b} A_{-b} + i \kappa_{ac} A_f e^{-i\Delta\beta z}
    \end{cases}
\end{equation}

\noindent where:
\begin{equation}
    \kappa_{dc,-b} = n_0 \omega \varepsilon_0 \iint_{-\infty}^\infty  \overline{\Delta n} \xi_{b} \xi_{b}^{\ast} dx dy
\end{equation}
\begin{equation}
    \kappa_{dc,f} = n_0 \omega \varepsilon_0 \iint_{-\infty}^\infty  \overline{\Delta n} \xi_{f} \xi_{f}^{\ast} dx dy
\end{equation}
\begin{equation}
    \kappa_{ac} = n_0 \omega \varepsilon_0 \iint_{-\infty}^\infty  \frac{\Delta n}{2} \xi_{f} \xi_{b}^{\ast} dx dy
\end{equation}
\begin{equation}
    \Delta\beta = \beta_b + \beta_f - \frac{2 \pi N}{\Lambda}
\end{equation}

The terms above describe specific interactions generated by the grating with the two modes. $\kappa_{dc,-b}$ and $\kappa_{dc,f}$ describe how strongly the grating couples the two modes into themselves. $\kappa_{ac}$ describes how much the modes couple between themselves. The last term, $\Delta\beta$, represents a detunement factor, which also indicates how strongly the modes couple between themselves. 

In the case of $f=b$, this set of equations can be solved by eliminating the $e^i\Delta\beta z$ term, then by decoupling the differential equations. $z=0$ is set to be the beginning of the grating and $z=L$ the end of the grating. The resultant amplitude reflection coefficient becomes:
\begin{equation}
    \rho = \frac{A_f(0)}{A_{-f}(0)} =  \frac{-\kappa_{ac}\sinh(\alpha L)}{\delta \sinh(\alpha L) - i\alpha \cosh(\alpha L)}
\end{equation}

Where:
\begin{equation}
    \delta = \kappa_{dc,f} + \frac{\Delta\beta}{2}
\end{equation}
\begin{equation}
    \alpha = \sqrt{|\kappa_{ac}|^2 - \delta^2}
\end{equation}

This results allows for describing the reflection in a simple function which can then be analysed to infer the effect each parameter has on the result. While this is useful for a heuristic description of the coupling between two specific modes, it ignores the coupling between all the other modes.

%%%%
\subsection{Detunement effect on coupling} \label{ssec:3}
This subsection will describe the effect on coupling of detunement between two fibre modes. This will provide support for any assumption which can be made going forward about which pairs of modes can be ignored, as their detunement is large enough to assume that the modes do not couple into each other.

The first step in understanding the detunement effect is to understand the effect $\Delta\beta$ has on the coupling. Figure \ref{fig:dbeta_3d} plots $|\rho|$ as a function of $\Delta\beta$ and $L$. This was done using $\kappa_{dc} = \kappa_{ac} = 1$ as a base value.

\begin{sidewaysfigure}
    \centering
    \includegraphics[width=\textwidth]{images/dbeta_L_rho.png}
    \caption{Plot of $|\rho|$ as a function of $\Delta\beta$ and $L$}
    \label{fig:dbeta_3d}
\end{sidewaysfigure}

It is clear from figure \ref{fig:dbeta_3d}, that the reflection coefficient drastically decreases up to that point and limiting the search range for coupling modes up to a $\Delta\beta$ of $1000 \frac{\mathrm{rad}}{\mathrm{m}}$ would be a good assumption. It is observed that a longer grating results in larger reflectance, as expected. Additionally, a shorter grating will suffer more from detuning, which is to be expected, as then shorter coupling distances become comparable to the grating length. The plot will be symmetrical about the $\Delta \beta$ axis.

For the fibre parameters $\lambda = 488$ nm, $n_{core} = 1.45$, NA $= 0.22$ of the MMF to be used in practice, calculations lead to $n_{core} k_0 = 1.87 * 10^{7} \;\frac{\mathrm{rad}}{\mathrm{m}}$ and $n_{cladding} k_0 = 1.845 * 10^{7} \;\frac{\mathrm{rad}}{\mathrm{m}}$. This result gives a maximum possible $\Delta\beta$ of $216,840 \;\frac{\mathrm{rad}}{\mathrm{m}}$ within the set of forward or backward propagating modes. Thus, the range of propagation constants for either FP or BP modes is two orders of magnitude larger that the search range limit placed after observing the result in figure \ref{fig:dbeta_3d}. This result shows that an FP mode will only couple into a small subset of BP modes. The significance of this result is that the specificity of the coupling generated by the grating is very high, leading to a good wavelength isolation for the reflection.

%%%%
% Subsection detailing the full-form of the mode coupling equations
\subsection{Multi-modal approach} \label{ssec:multi_modal_approach}
The derivation of the coupling described before only describes the coupling between one forward propagating mode and one backward propagating mode. However, in an MMF, there is a large number of modes propagating in both directions. Thus, in order to analyse how a Bragg grating behaves in an MMF, the interactions between the modes needs to be taken into consideration simultaneously.

The aim of this subsection is to represent the mode coupling as a linear system of ordinary differential equations, where the first order derivative of the mode coefficients with respect to $z$ is equal to a linear summation of the mode coefficients at a specific axial location, which can then be decoupled to allow for solving the system in a straightforward manner.

Following the same approach as in section \ref{ssec:bragg_theory_old}, however, this time without any assumptions being made, by expanding both $\delta n$ and $E$ terms in equation \ref{eq:mode_coupling_base_expanded} and computing the second order partial time derivative of $E$, the coupled mode equation which arises is given in equation \ref{eq:fullform} 

\begin{align} \label{eq:fullform}
    \frac{\partial A_m}{\partial z} \exp \left[ - iz \beta_m \right] & - \frac{\partial A_{-m}}{\partial z} \exp \left[ iz \beta_m \right] + cc =  \nonumber \\
    & - i \sum_{h} A_h \kappa_{dc,(h,m)} \exp \left[ iz \left( - \beta_h \right) \right] \nonumber \\
    & - i \sum_{h} A_h \kappa_{ac,(h,m)} \exp \left[ iz \left( - \beta_h + \frac{2 \pi N}{\Lambda} \right) \right] \nonumber \\
    & - i \sum_{h} A_h \kappa_{ac,(h,m)} \exp \left[ iz \left( - \beta_h - \frac{2 \pi N}{\Lambda} \right) \right] \nonumber \\
    & - i \sum_{-h} A_{-h} \kappa_{dc,(-h,m)} \exp \left[ iz  \beta_h \right] \nonumber \\
    & - i \sum_{-h} A_{-h} \kappa_{ac,(-h,m)} \exp \left[ iz \left( \beta_h + \frac{2 \pi N}{\Lambda} \right) \right] \nonumber \\
    & - i \sum_{-h} A_{-h} \kappa_{ac,(-h,m)} \exp \left[ iz \left( \beta_h - \frac{2 \pi N}{\Lambda} \right) \right] + cc
\end{align}

Where:
\begin{equation}
    \kappa_{dc,(a,b)} = n_0 \omega \varepsilon_0 \iint_{-\infty}^\infty \overline{\Delta n} \xi_a \xi_b^{\ast} dx dy
\end{equation}
\begin{equation}
    \kappa_{ac,(a,b)}= n_0 \omega \varepsilon_0 \iint_{-\infty}^\infty \frac{\Delta n}{2} \xi_a \xi_b^{\ast} dx dy
\end{equation}

As the profile of the FP modes is identical to the profile of the BP modes, then:
\begin{align}
     \kappa_{dc,(a,b)} = \kappa_{dc,(-a,b)} &= \kappa_{dc,(a,-b)} = \kappa_{dc,(-a,-b)} \nonumber \\
     \kappa_{ac,(a,b)} = \kappa_{ac,(-a,b)} &= \kappa_{ac,(a,-b)} = \kappa_{ac,(-a,-b)}
\end{align}

Also, as the only complex parameters in the $\kappa$ terms are the mode profiles, then it is easy to deduce that:
\begin{align}
     \kappa_{dc,(a,b)} & = \kappa_{dc,(b,a)}^\ast \nonumber \\
     \kappa_{ac,(a,b)} & = \kappa_{ac,(b,a)}^\ast
\end{align}

% Subsection  discussing the terms in the full-form equations
\subsubsection{Discussion on the terms of the full-form equation}

The amount each term in the full-form equations contribute to coupling one of the modes present on the left side of equation \ref{eq:fullform} is dependent on how close the $\exp(iz)$ terms are. As on the left side of equation \ref{eq:fullform} the exponential terms contain a $-\beta_m$ and a $+\beta_m$ term, there is a very large difference in the propagation constants. Thus, if a specific term on the right hand side is tuned to one term in the left hand side, it will be detuned to the other and not contribute to it, leading to the possibility of treating the coupling into the forward $m$ mode separately from the coupling into the $-m$ mode.

First, an overview of the propagation constants is important. In a multi-mode fibre, the propagation constants of the FP modes are in the range of $-n_{core} k_0$ to $-n_{clad} k_0$ and in the range of $n_{clad} k_0$ to $n_{core} k_0$ for the BP modes. Figure \ref{fig:beta_FP_BP_diff} illustrates these intervals, while also providing values for a fibre with $n_{core} = 1.45$ and $NA = 0.22$ transmitting light of wavelength: $\lambda = 488\ \rm{nm}$. The illustration will provide the basis for visualising detunement.

\begin{figure} [!htbp]
    \centering
    \includegraphics[width=\textwidth]{images/beta_line_1.png}
    \caption{Illustration of the difference between the propagation constants of forward and backward propagating modes}
    \label{fig:beta_FP_BP_diff}
\end{figure}

As it can be seen, the range of propagation constants available for propagating modes is much smaller than the absolute value of the propagation constant itself ($\approx 1\%$). Furthermore, based on the figure in section \ref{ssec:discussion_delta_beta}, appreciable coupling occurs for a detunement range of $\approx(-1000, 1000) \rm{\frac{rad}{m}}$, which is another 2 orders of magnitude lower. Thus, coupling is extremely specific in terms of phase matching compared to the range of propagation constants within a set of fibre modes.

Before discussing which right-hand side terms in equation \ref{eq:fullform} contribute to which left-hand side term or not at all, a discussion is required on what each right-hand side term represents physically. 

The terms which contain $\kappa_{dc}$ components (and implicitly $\overline{\Delta n}$) represent DC coupling generated by the DC component of the grating. Similarly, terms which contain $\kappa_{ac}$ components (and implicitly $\Delta n$) represent AC coupling generated by the AC component of the grating. 

In terms of propagation direction, the terms which contain a phase term of $-\beta_h$ represent coupling generated by FP modes, while terms which contain a phase term of $+\beta_h$ represent coupling generated by BP modes.

Lastly, AC terms are also split based on the sign of the harmonic order of the grating in question. Thus, terms with a phase evolution term of $+ \frac{2 \pi N}{\Lambda}$ represent a positive grating harmonic, while terms with a phase evolution term of $- \frac{2 \pi N}{\Lambda}$ represent a negative grating harmonic.

Table \ref{tab:eq_lines_coupling} provides a concise description of all of the right-hand side terms. The first column shows an index being assigned to each term. This will become useful in the discussion below. The last two columns indicate which terms contribute to coupling into which mode, which is discussed below. 

For a term on the right-hand side of equation \ref{eq:fullform} to represent coupling into the forward propagating mode $m$, it needs to phase match with $\exp \left[ - iz \beta_m \right]$. Similarly, for a right-hand side term to represent coupling into the backward propagating mode $-m$, it needs to phase match $\exp \left[ iz \beta_m \right]$.

\begin{table}[h]
    \centering
    \resizebox{\textwidth}{!}{%
    \begin{tabular}{|c|c|c|c|c|c|c|}
         \hline

          & Right-hand side term & DC/AC & FP/BP & $N$ sign & $A_m$ term & $A_{-m}$ term \\
         
         \hline
         
         1 & $- i \sum_{h} A_h \kappa_{dc,(h,m)} \exp \left[ iz \left( - \beta_h \right) \right]$ &  
         DC & FP & N/A & Yes & No \\ 
         
         \hline
         
         2 & $- i \sum_{h} A_h \kappa_{ac,(h,m)} \exp \left[ iz \left( - \beta_h + \frac{2 \pi N}{\Lambda} \right) \right]$ & 
         AC & FP & + & No & Yes \\ 
         
         \hline
         
         3 & $- i \sum_{h} A_h \kappa_{ac,(h,m)} \exp \left[ iz \left( - \beta_h - \frac{2 \pi N}{\Lambda} \right) \right]$ & 
         AC & FP & - & No & No \\ 
         
         \hline
         
         4 & $- i \sum_{-h} A_{-h} \kappa_{dc,(-h,m)} \exp \left[ iz  \beta_h \right]$ & 
         DC & BP & N/A & No & Yes \\ 
         
         \hline
         
         5 & $- i \sum_{-h} A_{-h} \kappa_{ac,(-h,m)} \exp \left[ iz \left( \beta_h + \frac{2 \pi N}{\Lambda} \right) \right]$ & 
         AC & BP & + & No & No \\ 
         
         \hline
         
         6 & $- i \sum_{-h} A_{-h} \kappa_{ac,(-h,m)} \exp \left[ iz \left( \beta_h - \frac{2 \pi N}{\Lambda} \right) \right]$ & AC & BP & - & Yes & No \\ 
         
         \hline
    \end{tabular}}
    \caption{Table describing the effect of each line of the right hand side of the equation \ref{eq:fullform}. The first column assigns an index to each term. The second column shows the right-hand side term in question. The ``DC/AC'' column describes if the term is generated by DC or AC coupling. The ``FP/BP'' column indicates if the term represents coupling generated by FP or BP modes. The ``$N$ sign'' column describes if the term represents a coupling generated by a positive or negative sign harmonic order. This does not apply for DC terms (it could also be considered a harmonic order of 0). The ``$A_m$ term'' and ``$A_{-m}$ term'' columns indicate if the term represents coupling into the $A_m$ and $A_{-m}$ modes, respectively.}
    \label{tab:eq_lines_coupling}
\end{table}

Due to the large difference in propagation constants (see fig. \ref{fig:beta_FP_BP_diff}), DC coupling cannot occur between counter-propagating modes (i.e. a FP mode cannot DC couple into a BP mode and vice-versa). Thus, term 1 only contributes to coupling into the $A_m$ mode (if  $-\beta_h \approx  -\beta_m$), while term 4 only contributes to coupling into the $A_{-m}$ mode (if  $\beta_h \approx  \beta_m$).

For co-propagation AC coupling to occur (terms 2 and 3 to contribute to coupling into $A_m$ mode and terms 5 and 6 to contribute to coupling into $A_{-m}$ mode), the requirement is: $\beta_h \approx \beta_m \pm \frac{2 \pi N}{\Lambda}$ for terms 5 and 6 and $-\beta_h \approx -\beta_m \pm \frac{2 \pi N}{\Lambda}$ for terms 2 and 3. 

The diagram in figure \ref{fig:beta_AC_co_prop} provides a visual representation of the requirement for co-propagation AC coupling.

\begin{figure} [!htbp]
    \centering
    \includegraphics[width=\textwidth]{images/beta_line_AC_co_prop.png}
    \caption{Diagram of co-propagation AC coupling requirement. Middle line shows the entire range of propagation factors for both FP and BP fibre modes. The top and bottom line is a zoomed-in view on the propagation constants range of the BP and FP modes, respectively. The top line illustration shows the requirement for the AC component of the grating to couple a BP mode with propagation constant $\beta_m$ into other BP modes with propagation constants $\beta_{h1}$ and $\beta_{h2}$ which are co-propagating. The bottom line shows the same requirement for FP modes.}
    \label{fig:beta_AC_co_prop}
\end{figure}

It is clear that for co-propagation AC coupling to occur, $\frac{2 \pi N}{\Lambda}$ must be small, thus, the grating period must be long. For the fiber parameters given at the beginning of the subsection, in order for such coupling to even be possible the grating period would need to be: $\Lambda \geq \frac{2\pi}{(n_{core} - n_{clad}) k_0} \approx 28 \rm{\mu m}$.

For counter-propagation AC coupling to occur (terms 2 and 3 to contribute to coupling into $A_{-m}$ mode and terms 5 and 6 to contribute to coupling into $A_m$ mode), the requirement is: $-\beta_h \approx \beta_m \pm \frac{2 \pi N}{\Lambda}$ for terms 2 and 3 and $\beta_h \approx -\beta_m \pm \frac{2 \pi N}{\Lambda}$ for terms 5 and 6.

The diagram in figure \ref{fig:beta_AC_counter_prop} provides a visual representation of the requirement for AC counter-propagation coupling for gratings with a small period.

The first observation is that FP mode into BP mode AC coupling is only provided by positive sign harmonics, thus, only term 2 couples into the $A_{-m}$ mode and term 3 does not. Similarly, BP mode into FP mode AC coupling is only provided by negative sign harmonics, thus, only term 6 couples into the $A_m$ mode and term 5 does not.

In order for this type of coupling to occur, either the grating period needs to be small, or the harmonic order large.
In the case of a small grating period, where a low-order harmonic provides counter-propagation coupling, co-propagation AC coupling can be assumed to be impossible (terms 2 and 3 do not contribute to the $A_m$ mode, while terms 5 and 6 do not contribute to $A_{-m}$ mode). Furthermore, only one harmonic provides such coupling, as other orders would fall outside fibre mode propagation range.

\begin{figure} [!htbp]
    \centering
    \includegraphics[width=\textwidth]{images/beta_line_AC_counter_prop.png}
    \caption{Diagram of counter-propagation AC coupling requirement. The top line shows the requirement for the AC component of the grating to couple a FP mode with propagation factor $-\beta_m$ into a BP mode with propagation factor $\beta_{h2}$. Technically, a mode with propagation constant $-\beta_{h1}$ could also be coupled into, however, there are no propagation fibre modes with such a propagation factor as they would be outside the propagation range of the fibre. The same is shown with the second line for the case of a BP mode.}
    \label{fig:beta_AC_counter_prop}
\end{figure}

For both co-propagation and counter-propagation AC coupling to occur, a large grating period is required, however, the grating also needs to have very-high order harmonic components.

For the fibre parameters given at the beginning of the subsection and given the limitation of $\Lambda \geq 28 \rm{\mu m}$ required for co-propagation Ac coupling, then counter-propagation coupling would occur at harmonics orders of: $N > 82$. In such a case, it is also important to note, that multiple harmonics could provide counter-propagation coupling (i.e. if the harmonic of order $N$ provides counter-propagation coupling, then the harmonics of orders \{... $N-2$, $N-1$, $N+1$, $N+2$ ...\} could also provide coupling, as the step difference between them of $\frac{2\pi}{\Lambda}$ is smaller than the range of propagation constants of the fibre modes of one type of propagation: BP or FP).

My work focuses on small grating periods, where co-propagation AC coupling does not occur and only a single harmonic provides counter-propagation AC coupling, each term contribution to coupling into either the FP or BP modes being indicated in the last two columns of table \ref{tab:eq_lines_coupling}.

Thus, the simplified and separated form of the equation \ref{eq:fullform} is:

\begin{align} \label{eq:simple_F}
    \frac{\partial A_m}{\partial z} + cc & =  \nonumber \\
    & - i \sum_{h} A_h \kappa_{dc,(h,m)} \exp \left[ iz \left( \beta_m - \beta_h \right) \right] \nonumber \\
    & - i \sum_{-h} A_{-h} \kappa_{ac,(-h,m)} \exp \left[ iz \left( \beta_m + \beta_h - \frac{2 \pi N}{\Lambda} \right) \right] + cc
\end{align}

\begin{align} \label{eq:simple_B}
    \frac{\partial A_{-m}}{\partial z} + cc & = \nonumber \\
    & + i \sum_{h} A_h \kappa_{ac,(h,-m)} \exp \left[ iz \left( - \beta_m - \beta_h + \frac{2 \pi N}{\Lambda} \right) \right] \nonumber \\
    & + i \sum_{-h} A_{-h} \kappa_{dc,(-h,-m)} \exp \left[ iz \left( - \beta_m + \beta_h) \right) \right] + cc
\end{align}

%%%%
\subsection{Generating a linear state matrix}
The purpose of this section is to reduce the forms of equations \ref{eq:simple_F} and \ref{eq:simple_B} to a state matrix form $\mathbf{x}' = \mathbf{M} \mathbf{x}$. If matrix $\mathbf{M}$ is not dependent on $z$, then an eigenvalue approach would allow for decoupling the equations and writing them into a simple set of first-order ordinary differential equations.

The first target is to eliminate the $e^{iz}$ dependence in the state matrix. The trick is to have the exponential terms within the new coefficients in the column vectors, rather than in the matrix:
\begin{equation}
    B_m = A_m e^{i(-\beta_m + \frac{\pi N}{\Lambda}) z} \Leftrightarrow A_m = B_m e^{i (\beta_m - \frac{\pi N}{\Lambda}) z}
\end{equation}
\begin{equation}
    B_{-m} = A_{-m} e^{i(\beta_m - \frac{\pi N}{\Lambda}) z} \Leftrightarrow A_{-m} = B_{-m} e^{i (-\beta_m + \frac{\pi N}{\Lambda}) z}
\end{equation}

First, let's introduce $B$, which is a function of $z$ as was $A$, to equation \ref{eq:simple_F}. The terms on the left hand side will become:

\begin{align}
    \frac{\partial A_m}{\partial z} & = \frac{\partial}{\partial z} \left( B_m \exp \left[ i (\beta_m - \frac{\pi N}{\Lambda}) z \right] \right) \nonumber \\
    & = \frac{\partial B_m}{\partial z} \exp \left[ i (\beta_m - \frac{\pi N}{\Lambda}) z \right] + i \left( \beta_m - \frac{\pi N}{\Lambda} \right) B_m \exp \left[ i (\beta_m - \frac{\pi N}{\Lambda}) z \right]
\end{align}

The first term on the right hand side of equation \ref{eq:simple_F} will become:
\begin{align}
    & - i \sum_{h} A_h \kappa_{dc,(h,m)} \exp \left[ iz \left( \beta_m - \beta_h \right) \right] = \nonumber \\
    = & - i \sum_{h} B_h \kappa_{dc,(h,m)} \exp \left[ iz \left(\beta_h - \frac{\pi N}{\Lambda} \right) \right] \exp \left[ iz \left( \beta_m - \beta_h \right) \right] \nonumber \\
    = & - i \sum_{h} B_h \kappa_{dc,(h,m)} \exp \left[ iz \left(\beta_m - \frac{\pi N}{\Lambda} \right) \right]
\end{align}

Lastly, the second term on the right hand side of equation \ref{eq:simple_F} will become:
\begin{align}
    &   \sum_{-h} -i\kappa_{ac,(m,-h)}^\ast \exp \left[ i \left(\beta_m + \beta_h - \frac{2\pi N}{\Lambda}\right) z\right] A_{-h} = \nonumber \\
    & = \sum_{-h} -i\kappa_{ac,(m,-h)}^\ast \exp \left[ i \left(\beta_m + \beta_h - \frac{2\pi N}{\Lambda} \right) z \right] B_{-h} \exp \left[i \left(-\beta_h + \frac{\pi N}{\Lambda} \right) z \right] \nonumber \\
    & = \sum_{-h} -i\kappa_{ac,(m,-h)}^\ast \exp \left[ i \left( \beta_m - \frac{\pi N}{\Lambda} \right) z \right] B_{-m}
\end{align}

Combining the new form of the terms in equation \ref{eq:simple_F}, the $e^{iz}$ factor can be eliminated:
\begin{gather}
    \frac{\partial B_m}{\partial z} = -i \left( \beta_m - \frac{\pi N}{\Lambda} \right) B_m + \sum_{h} - i \kappa_{dc,(h,m)} B_h + \sum_{-h} -i\kappa_{ac,(m,-h)}^\ast  B_{-h}
\end{gather}

Applying the same approach for $\frac{\partial A_{-m}}{\partial z}$ leads to:

\begin{equation}
    \frac{\partial B_{-m}}{\partial z} = i \left( \beta_m - \frac{\pi N}{\Lambda} \right) B_{-m} + \sum_{-h} i \kappa_{dc,(-h,-m)} B_{-h} + \sum_{h} i\kappa_{ac,(h,-m)}  B_{h} 
\end{equation}

Thus, the equations above can be written as a linear state matrix. An example for two FP and two BG modes is given below:
\begin{equation} \label{eq:explicit_state_system}
\frac{\partial}{\partial z}
\begin{pmatrix}
        B_1 \\
        B_2 \\
        B_{-1} \\
        B_{-2}
\end{pmatrix}
= 
\bm{M_B}
\begin{pmatrix}
        B_1 \\
        B_2 \\
        B_{-1} \\
        B_{-2}
\end{pmatrix}
\end{equation}

\noindent where:

\begin{equation}
\bm{M_B} = \nonumber
\end{equation}

\begin{equation}
\noindent\resizebox{\columnwidth}{!}{%
\(
\begin{pmatrix}
    -i \left( \kappa_{dc,(1,1)} + \beta_1 - \frac{\pi N}{\Lambda} \right) & -i \kappa_{dc,(2,1)} & -i \kappa_{ac,(-1,1)} & -i\kappa_{ac,(-2,1)} \\
    -i \kappa_{dc,(1,2)} & -i \left( \kappa_{dc,2} + \beta_2 - \frac{\pi N}{\Lambda} \right) & -i \kappa_{ac,(-1,2)} & -i\kappa_{ac,(-2,2)} \\
    i \kappa_{ac,(1,-1)} & i \kappa_{ac,(2,-1)} & i \left( \kappa_{dc,(-1,-1)} + \beta_1 - \frac{\pi N}{\Lambda} \right) & i \kappa_{dc,(-2,-1)} \\
    i \kappa_{ac,(1,-2)} & i \kappa_{ac,(2,-2)} & i \kappa_{dc,(-1,-2)} & i \left( \kappa_{dc,(-2,-2)} + \beta_2 - \frac{\pi N}{\Lambda} \right) \\
\end{pmatrix}
\) }
\end{equation}


This result shows that the $z$ dependence can be eliminated from the matrix, allowing for the system to be treated as an usual state matrix system.

%%
\subsubsection{Solving the linear state system}

The result in the last section can be summarised by the equation below, where $\mathbf{b}$ represents the vectors of $B_m$ (as was shown explicitly in equation \ref{eq:explicit_state_system} for a fibre with two modes).
\begin{equation}
    \frac{\partial}{\partial z} \bm{b} = \bm{M_B} \bm{b}
\end{equation}

If $\bm{M_B}$ is diagonalisable, then the system of differential equation can be decoupled, where $\mathbf{D}$ is a diagonal matrix of the eigenvalues of $\mathbf{M_B}$, $\mathbf{V}$ is the matrix of eigenvectors of $\mathbf{M_B}$ and the notation $\mathbf{p} = \mathbf{V}^{-1}\mathbf{b}$ is used.

\begin{equation}
    \bm{M_B} = \bm{ V D V^{-1}}
\end{equation}

\begin{equation}
    \frac{\partial}{\partial z} (\bm{V^{-1} b}) = \bm{D V^{-1} b}
\end{equation}

\begin{equation}
    \frac{\partial}{\partial z} \bm{p} = \bm{D p}
\end{equation}

This is a set of decoupled ordinary differential equations, which can be solved through direct integration.

\begin{equation}
    \frac{d p_h(z)}{d z} = \lambda_h p_h(z) \Leftrightarrow \frac{d p_h(z)}{p_h(z)} = \lambda_h dz
\end{equation}

\begin{equation}
    p_h(z) = p_h(0) e^{\lambda_h z}
\end{equation}

\begin{equation}
    \bm{p}(z) = e^{\bm{D}z} \bm{p}(0)
\end{equation}

\begin{align} \label{eq:diag_form_scaterring}
    \bm{V^{-1}b}(z) & = e^{\bm{D}z} \bm{V^{-1}b}(0) \nonumber \\
    \bm{b}(z) & = \bm{V} e^{\bm{D}z} \bm{V^{-1}b}(0)
\end{align}

Before focusing on initial conditions, a discussion on the relationship between $\bm{a}$, $\bm{b}$ and the total transverse electric field needs to be made.

In the analytical approach, the total transverse electric field was taken as:
\begin{equation}
    E_t(z) = \sum_{h=1}^{N_{modes}} A_h e^{i(\omega t - \beta_h z)} +  \sum_{-h=-1}^{-N_{modes}} A_{-h} e^{i(\omega t + \beta_h z)}
\end{equation}

Thus, it is clear that $\bm{A}$ does not take into account the phase shift due to propagation. However, using the relationship between $\bm{A}$ and $\bm{B}$ from last section, the transverse electric field can be written as:

\begin{gather}
    E_t(z) = \sum_{h=1}^{N_{modes}} B_h e^{i (\beta_h - \frac{\pi N}{\Lambda}) z} e^{i(\omega t - \beta_h z)} +  \sum_{-h=-1}^{-N_{modes}} B_{-h} e^{i (-\beta_h + \frac{\pi N}{\Lambda}) z} e^{i(\omega t + \beta_h z)} \nonumber \\
    E_t(z) = \sum_{h=1}^{N_{modes}} B_h e^{i (\omega t - \frac{\pi N}{\Lambda} z)} + \sum_{-h=-1}^{-N_{modes}} B_{-h} e^{i (\omega t + \frac{\pi N}{\Lambda} z)}
\end{gather}

Thus, $\bm{b}$ represents the mode coefficients with propagation taken into account, but with an extra opposite phase shift for FP and BP modes. This allows for applying initial conditions directly to the $\bm{b}$ coefficients.

The initial conditions applied will be setting $E_{t,FP}(z=0)$ to the field at the beginning of the grating and $E_{t,BP}(z=L)$ the back-propagating field at the end of the grating. In practice, the latter can be set to 0, assuming no other reflections occur after the grating and that no light is being inputted from the distal facet.

In other words, the initial conditions are $B_h(z=0)$ and $B_{-h}(z=L)$ (which could be 0).

Let the parameters be split between BP (-) and FP (+) modes:
\begin{equation}
    \bm{b} = 
    \begin{pmatrix}
        \bm{b_+} \\
        \bm{b_-}
    \end{pmatrix}
\end{equation}

\begin{equation}
    \bm{V} e^{\bm{D}L} \bm{V^{-1}} = 
    \begin{pmatrix}
        \bm{S_{++}} & \bm{S_{+-}} \\
        \bm{S_{-+}} & \bm{S_{--}} \\
    \end{pmatrix}
\end{equation}

Then:
\begin{equation}
    \begin{cases}
        \bm{b_+}(z=L) = \bm{S_{++}b_+}(z=0) + \bm{S_{+-}b_-}(z=0) \\
        \bm{b_-}(z=L) = \bm{S_{-+}b_+}(z=0) + \bm{S_{--}b_-}(z=0)
    \end{cases}
\end{equation}
Which leads to:
\begin{equation}
    \begin{cases}
        \bm{b_-}(z=0) = \bm{S_{--}}^{-1} \left[ \bm{b_-}(z=L) - \bm{S_{-+}b_+}(z=0) \right] \\
        \bm{b_+}(z=L) = \bm{S_{++}b_+}(z=0) + \bm{S_{+-}b_-}(z=0)
    \end{cases}
\end{equation}

%%%%
% Subsection discussing the use of a non-sinusoidal grating
\subsection{Discussion on non-sinusoidal gratings} \label{ssec:discussion_delta_beta}
The approach in section \ref{ssec:bragg_theory_old} as well as the approach detailed in subsection \ref{ssec:multi_modal_approach} focus on sinusoidal gratings. In practice, such gratings would be hard to manufacture. However, any periodical grating can be described through the use of Fourier series. Or to be more specific:

\begin{equation}
    \delta n = a_0 + \sum_{k=1}^\infty a_k \cos \left( \frac{2 \pi k}{\Lambda} z \right) + \sum_{k=1}^\infty b_k \sin \left( \frac{2 \pi k}{\Lambda} z \right)
\end{equation}

If the assumption is made that only the DC component and a single harmonic order provide coupling, due to choosing the grating to specifically filter a wavelength, then the grating will act as if its a sinusoid:

\begin{equation}
    \delta n = a_0 + a_N \cos \left( \frac{2 \pi N}{\Lambda} z \right) + b_N \sin \left( \frac{2 \pi N}{\Lambda} z \right)
\end{equation}

This can be further rewritten as:
\begin{equation}
    \delta n = a_0 + \sqrt{a_N^2 + b_N^2} \cos \left( \frac{2 \pi N}{\Lambda} z + \phi \right)  
\end{equation}

where:
\begin{equation}
    \sin (\phi) = \frac{-b_N}{\sqrt{a_N^2 + b_N^2}}
\end{equation}

If the assumption is made that the core refractive index is one of the limits of the refractive index variation (or similarly 0 is either the upper or lower limit of $\delta n$) and that the variation is symmetric with respect to z, then $\delta n$ is an even function, where $z=0$ is the start of the grating, thus, the sine terms in the Fourier series become 0, leading to:

\begin{equation} \label{eq:delta_N_fourier_even}
    \delta n = a_0 + a_N \cos \left( \frac{2 \pi N}{\Lambda} z \right)
\end{equation}

The form in equation \ref{eq:delta_N_fourier_even} is identical to equation \ref{eq:delta_N_original}.

%%%%
\subsection{Implementation of the theory in a simulation}

%%%
\subsubsection{Focus simulation} \label{ssec:bragg_focus_simulation}
The theory described has been implemented in Matlab. 
In terms of the simulation parameters, the grating chosen was a line across the fibre core as shown in figure \ref{fig:sample_16_base}. While the discussion in the previous chapters have used reflectors very small in size, the first simulation has been done on a grating with a different geometry. The reason for such change is that a line grating is closer to what can already be more easily achieved in practice, thus, to provide and initial investigation into a grating which has a higher change of being tested experimentally. However, this chapter will also investigate later gratings which are small in size. The target of the first simulation was to test the approach, that it gives a reasonable answer and that it can be used for localising the focus, whether it is on the grating or not.

The grating parameters were:
\begin{itemize}
    \item Pitch: $534.9$ nm, designed for a Bragg wavelength of $515$ nm (third order - $N=3$)
    \item Cross-sectional area: $\sim 30 \;\mathrm{\mu m}$ by $ \sim 2 \;\mathrm{\mu m}$
    \item Length of the grating: $\sim 5 \;\mathrm{mm}$
\end{itemize}

The grating profile was chosen to be a two-dimensional Gaussian distribution, setting the maximum change in refractive index at 0.02. The relationship between the DC refractive index change and the amplitude of the AC refractive index change was set as: $\overline{\Delta n} = - \Delta n$, such that the change in refractive index $\delta n$ at the beginning of the grating at $z=0$ is 0 (see equation \ref{eq:delta_N_original}).

\begin{figure}
    \centering
    \includegraphics[width = 0.5\textwidth]{images/sample_16_refr_sim.png}
    \caption{The refractive index change profile of the simulated sample. The colours in the figure represent: dark blue - the cladding of the fibre, cyan to yellow - the amplitude of the sinusoidal change in refractive index profile ($\Delta n$) from 0 (cyan) to 0.02 (yellow). The value of the DC change is the opposite ($\overline{\Delta n} = -\Delta n$) in order for the total refractive index change at z=0 to be 0.}
    \label{fig:sample_16_base}
\end{figure}

\begin{figure}
    \centering
    \includegraphics[width = \textwidth]{images/sample_16_result.png}
    \caption{Initial simulation result for the Bragg grating, given by focusing input light at a location away from the grating.}
    \label{fig:sample_16_sim_res}
\end{figure}

The simulation was done after generating the set of fibre modes for a light wavelength identical to the Bragg wavelength. Then the procedure described in the previous subsections was followed.

Only the portion of the fibre containing the grating was simulated by specifying the optical field at the beginning of the grating. The first set of initial conditions, the forward input optical field at $z=0$, was generated by focusing input light at different positions. The second set of initial conditions, the backward optical field at $z=L$, was chosen such that there was no light entering the grating from the output end. The target result of the simulation is the light exiting the grating at both ends of the fibre section.

The initial result for an input focus is given in figure \ref{fig:sample_16_sim_res}. The ``Forward-in'' panel shows the location of the input focus. The ``Forward-out'' and ``Backward-in'' panels show the light exiting at the output end and the input end, respectively. The last quadrant (``Backward-out'') represents the light entering in the backward direction at the output end, which is shown to be set to 0.

It is clear that there is a problem. The input focus intensity is in the range of $2\cdot 10^{-3}$, while the forward light is in the range of $10^{215}$. The reflected light is also several orders of magnitude larger than the input light. Thus, the simulation shows an increase in total energy, which is not possible in reality.

In order to understand why this effect occurs, a look has been given at the eigenvalues of the matrix $\mathbf{M_B}$, the results being shown in figure \ref{fig:eigenvalue_MB}.

The blue plot shows the cosine of the argument of the eigenvalue. This illustrates which eigenvalues are purely imaginary (value of 0), purely real (value of $\pm 1$), or in-between. While most eigenvalues are purely imaginary, some of them are not.

However, in order to understand better where which eigenvalues contribute more to the exponential amplification, the orange plot shows on a logarithmic scale the real part of the eigenvalue. The absolute value of the real part has been used in order to include negative value as well. While a negative real part might seem to imply exponential decay, the forward and backward propagations are in opposite directions, thus the component of one will decay while the other will amplify as it propagates.

It is clear that the eigenvalues are split into two groups: the first group has real values in the $10^{-10}$ range, which are within numerical errors of a purely imaginary eigenvalue, while the second group are in the range of $10^0$ to $10^5$. The second group and more specifically the eigenvalues which have a larger real part are the problem in this situation.  

\begin{sidewaysfigure}
    \centering
    \includegraphics[width = \textwidth]{images/eigenvalue_MB.png}
    \caption{A plot of the eigenvalues of the matrix $\mathbf{M_B}$ for the simulated grating. The horizontal axis represents the index of the eigenvalue. The left vertical axis (blue) represents the cosine of the argument of the eigenvalue. The right vertical axis (orange) represents a logarithmic plot of the absolute value of the real part of the eigenvalue.}
    \label{fig:eigenvalue_MB}
\end{sidewaysfigure}

An important note to mention is that the sum of the eigenvalues is 0, as the trace of $\mathbf{M_B}$ is 0. This implies that the determinant of $\mathbf{S}=e^{\mathbf{M_B}L}$ is 1. In practice, the resulting $\mathbf{S}$ had a determinant of $+\infty$, which is the value allocated by Matlab for a very large determinant. This might mean that there is some numerical error which occurs due to the large real parts of the eigenvalues when calculating the matrix exponential. Thus, it could be that working in the scattering modal space might be a better approach to have greater care in using the very large eigenvalues.

An interesting idea was observed in \cite{scattering_approach}, where the propagation of evanescent fields in stratified media was encountered. The main observation was that due to the exponential evolution of such fields, a numerical stable approach was to consider light fields exiting out of the system at both ends as a function of light fields entering into the system at both ends. This is different to the approach described which considers the light fields, entering and exiting, at one end of the optical system as a function of the light fields, entering and exiting, at the other end. Their approach focused on treating each field individually, however, in the Bragg grating case this would be very complex to compute, as each period of the grating would represent a section of the system. However, the idea of relating inputs to outputs can still be used.

Taking the equation \ref{eq:diag_form_scaterring} and using the notations $\mathbf{q} = \mathbf{V}^{-1} \mathbf{b}(z=L)$ and $\mathbf{p} = \mathbf{V}^{-1} \mathbf{b}(z=0)$ leads to:

\begin{equation} \label{eq:scat_pq}
    \mathbf{q} = e^{\mathbf{D}L} \mathbf{p}
\end{equation}

If the rows of the equation are reordered and separated into positive (index +) and negative (index -) real parts of the eigenvalues' exponent, equation \ref{eq:scat_pq} can be written as:

\begin{equation}
    \begin{pmatrix}
        \bm{q_+} \\
        \bm{q_-}
    \end{pmatrix} = 
    \begin{pmatrix}
        e^{\mathbf{D_+}L} & 0 \\
        0 & e^{\mathbf{D_-}L}
    \end{pmatrix} 
    \begin{pmatrix}
        \bm{p_+} \\
        \bm{p_-}
    \end{pmatrix}
\end{equation}

The individual exponent terms in the matrix are diagonal matrices, thus, the positive real part rows can be inverted to only work with negative real parts and avoid exponential gain (i.e. a constituent equation in the form $q_{k} = e^{D_k L} p_{k}$ was changed to $p_{k} = e^{-D_k L} q_{k}$, where $D_k > 0$).

\begin{equation} \label{eq:scat_pq_neg}
    \begin{pmatrix}
        \bm{p_+} \\
        \bm{q_-}
    \end{pmatrix} = 
    \begin{pmatrix}
        e^{-\mathbf{D_+}L} & 0 \\
        0 & e^{\mathbf{D_-}L}
    \end{pmatrix} 
    \begin{pmatrix}
        \bm{q_+} \\
        \bm{p_-}
    \end{pmatrix}
\end{equation}

The next step is to expand $\mathbf{p}$ and $\mathbf{q}$, however, beforehand, the matrix $\mathbf{V}^{-1}$ needs to be reordered to match the new order of $\mathbf{p}$ and $\mathbf{q}$. The order of the columns of $\mathbf{V}^{-1}$ was kept the same, as the order of the modes ($\mathbf{b}(L)$ and $\mathbf{b}(0)$) was not changed. A reminder is given that they were split into two sets based on the propagation direction of the fibre modes they are related to ($b+$ for FP modes and $b-$ for BP modes). The order of the rows of $\mathbf{V}^{-1}$ was changed to separate and split the eigenvectors they relate to into two sets based on the sign of the real part of their associated eigenvalue exponent ($r+$ for positive real part and $r-$ for negative real part). This allows for the new reordered $\mathbf{V}^{-1}$ to be separated into four sub-matrices:

\begin{equation}
    \mathbf{q}_{reordered} = \mathbf{V}_{reordered}^{-1} \mathbf{b}(L) \;\;\; \mathbf{p}_{reordered} = \mathbf{V}_{reordered}^{-1} \mathbf{b}(0) \nonumber
\end{equation}
\begin{equation}
    \begin{pmatrix}
        \bm{q_+} \\
        \bm{q_-}
    \end{pmatrix} = 
    \begin{pmatrix}
        \bm{V}_{r+b+}^{-1} & \bm{V}_{r+b-}^{-1} \\
        \bm{V}_{r-b+}^{-1} & \bm{V}_{r-b-}^{-1}
    \end{pmatrix} 
    \begin{pmatrix}
        \bm{b_+}(L) \\
        \bm{b_-}(L)
    \end{pmatrix} \nonumber
\end{equation}
\begin{equation}
    \begin{pmatrix}
        \bm{p_+} \\
        \bm{p_-}
    \end{pmatrix} = 
    \begin{pmatrix}
        \bm{V}_{r+b+}^{-1} & \bm{V}_{r+b-}^{-1} \\
        \bm{V}_{r-b+}^{-1} & \bm{V}_{r-b-}^{-1}
    \end{pmatrix} 
    \begin{pmatrix}
        \bm{b_+}(0) \\
        \bm{b_-}(0)
    \end{pmatrix}
\end{equation}

Writing $\mathbf{p}_-$, $\mathbf{p}_+$, $\mathbf{q}_-$ and $\mathbf{q}_+$ in terms of $\mathbf{b}_-(0)$, $\mathbf{b}_+(0)$, $\mathbf{b}_-(L)$ and $\mathbf{b}_+(L)$  and introducing them into equation \ref{eq:scat_pq_neg} leads to:

\begin{equation}
    \begin{pmatrix}
        \bm{p_+} \\
        \bm{q_-}
    \end{pmatrix} = 
    \begin{pmatrix}
        \bm{V}_{r+b+}^{-1} & \bm{V}_{r+b-}^{-1} & 0 & 0 \\
        0 & 0 & \bm{V}_{r-b+}^{-1} & \bm{V}_{r-b-}^{-1}         
    \end{pmatrix} 
    \begin{pmatrix}
        \bm{b_+(0)} \\
        \bm{b_-(0)} \\
        \bm{b_+(L)} \\
        \bm{b_-(L)}
    \end{pmatrix} 
    = \mathbf{M}_1
    \begin{pmatrix}
        \bm{b_+(0)} \\
        \bm{b_-(0)} \\
        \bm{b_+(L)} \\
        \bm{b_-(L)}
    \end{pmatrix} \nonumber
\end{equation}

\begin{equation} 
    \begin{pmatrix}
        \bm{q_+} \\
        \bm{p_-}
    \end{pmatrix} = 
    \begin{pmatrix}
        0 & 0 & \bm{V}_{r+b+}^{-1} & \bm{V}_{r+b-}^{-1}  \\
        \bm{V}_{r-b+}^{-1} & \bm{V}_{r-b-}^{-1} & 0 & 0       
    \end{pmatrix} 
    \begin{pmatrix}
        \bm{b_+}(0) \\
        \bm{b_-}(0) \\
        \bm{b_+}(L) \\
        \bm{b_-}(L)
    \end{pmatrix}
    = \mathbf{M}_2
    \begin{pmatrix}
        \bm{b_+}(0) \\
        \bm{b_-}(0) \\
        \bm{b_+}(L) \\
        \bm{b_-}(L)
    \end{pmatrix} \nonumber
\end{equation}

\begin{equation} \label{eq:scat_big_form}
    \mathbf{M}_1 
    \begin{pmatrix}
        \bm{b_+(0)} \\
        \bm{b_-(0)} \\
        \bm{b_+(L)} \\
        \bm{b_-(L)}
    \end{pmatrix} = 
    \begin{pmatrix}
        e^{-\mathbf{D_+}L} & 0 \\
        0 & e^{\mathbf{D_-}L}
    \end{pmatrix} 
    \mathbf{M}_2
    \begin{pmatrix}
        \bm{b_+(0)} \\
        \bm{b_-(0)} \\
        \bm{b_+(L)} \\
        \bm{b_-(L)}
    \end{pmatrix}
\end{equation}

In equation \ref{eq:scat_big_form}, taking the right hand side to the left leads to:
\begin{equation}
    \left[
    \mathbf{M}_1 - 
    \begin{pmatrix}
        e^{-\mathbf{D_+}L} & 0 \\
        0 & e^{\mathbf{D_-}L}
    \end{pmatrix} 
    \mathbf{M}_2
    \right]
    \begin{pmatrix}
        \bm{b_+(0)} \\
        \bm{b_-(0)} \\
        \bm{b_+(L)} \\
        \bm{b_-(L)}
    \end{pmatrix}
    =0 \nonumber
\end{equation}
\begin{equation}
    \bm{M}
    \begin{pmatrix}
        \bm{b_+(0)} \\
        \bm{b_-(0)} \\
        \bm{b_+(L)} \\
        \bm{b_-(L)}
    \end{pmatrix}
    =0
\end{equation}

The matrix $\mathbf{M}$ can be calculated using only the eigenvalues and eigenvectors of $\mathbf{M}_B$ and, due to only containing exponents with negative real parts, it is avoiding any numerical issues occurring due to exponential gain. Splitting $\mathbf{M}$ along the columns into 4 sub-matrices, each associated with a component of the column vector and taking $\mathbf{b}_+(0)$ and $\mathbf{b}_-(L)$ as known initial conditions leads to the solution in equation \ref{eq:scat_solution}. In terms of dimensions, $\mathbf{M}$ is a matrix of $2N_m$ rows by $4N_m$ columns, where $N_m$ is the number of fibre modes. Thus, each sub-matrix will be of size $2N_m$ by $N_m$ and two sub-matrices together form a square matrix, allowing for the inversion to occur, assuming the determinant is not 0.

\begin{equation}
    \begin{pmatrix}
        \bm{M}_{b+(0)} & \bm{M}_{b-(0)} & \bm{M}_{b+(L)} & \bm{M}_{b-(L)}  
    \end{pmatrix}
    \begin{pmatrix}
        \bm{b_+(0)} \\
        \bm{b_-(0)} \\
        \bm{b_+(L)} \\
        \bm{b_-(L)}
    \end{pmatrix}
    =0 \nonumber
\end{equation}
\begin{equation}
    \begin{pmatrix}
        \bm{M}_{b+(0)} & \bm{M}_{b-(L)} 
    \end{pmatrix}
    \begin{pmatrix}
        \bm{b_+}(0) \\
        \bm{b_-}(L)
    \end{pmatrix}
    +
    \begin{pmatrix}
        \bm{M}_{b-(0)} & \bm{M}_{b+(L)} 
    \end{pmatrix}
    \begin{pmatrix}
        \bm{b_-}(0) \\
        \bm{b_+}(L)
    \end{pmatrix}
    =0 \nonumber
\end{equation}

\begin{equation} \label{eq:scat_solution}
    \begin{pmatrix}
        \bm{b_-}(0) \\
        \bm{b_+}(L)
    \end{pmatrix} = 
    \begin{pmatrix}
        \bm{M}_{b-(0)} & \bm{M}_{b+(L)} 
    \end{pmatrix} ^ {-1}
    \begin{pmatrix}
        -\bm{M}_{b+(0)} & -\bm{M}_{b-(L)} 
    \end{pmatrix}
    \begin{pmatrix}
        \bm{b_+}(0) \\
        \bm{b_-}(L)
    \end{pmatrix}
\end{equation}

The approach described above has been implemented into code and the previously described grating was simulated, the results being shown in figure \ref{fig:5mm_sim}. Two cases were simulated, where the focus was either on or away from the grating to observe the ability of the grating to localise the reflection. It is clear that for a 5 mm grating there is no localisation occurring, in both cases the amount of light reflected being similar to the amount of light passing through. However, the grating seems to be confining the light within or outside its cross-sectional area. This implies that the grating behaves as a waveguide within the fibre core.

\begin{figure}
    \centering
    \begin{subfigure}[t]{0.75\textwidth}
        \includegraphics[width = \textwidth]{images/5mm_ONgr_lin.png}
        \caption{The simulation result for a focus on the grating.}
        \label{sfig:5mm_ON}
    \end{subfigure}
    \begin{subfigure}[t]{0.75\textwidth}
        \includegraphics[width = \textwidth]{images/5mm_OFFgr_lin.png}
        \caption{The simulation result for a focus away from the grating.}
        \label{sfig:5mm_AWAY}
    \end{subfigure}
    \caption{Simulation results for focusing light on and away from the 5 mm long grating.}
    \label{fig:5mm_sim}
\end{figure}

To illustrate the reason why a longer grating cannot localise the reflection, figure \ref{fig:focus_evo} shows the evolution of the light field inside the core of a fibre given a focus at a specific location within the core. It is clear that, while the light is localised close to the focus location, it expands to almost the entire core within 200 $\mathrm{\mu m}$ of the focus plane. Thus, a grating 5 mm long would interact with the light field no matter where the light is focused as, eventually, within its length, the light will reach the grating.

\begin{sidewaysfigure}
    \centering
    \begin{subfigure}[t]{\textwidth}
        \includegraphics[width = \textwidth]{images/focus_evolution_centre.png}
        \caption{The axial evolution of the light field for a focus in the centre of the core of the fibre.}
        \label{sfig:focus_evo_centre}
    \end{subfigure}
\end{sidewaysfigure}
\begin{sidewaysfigure}
    \ContinuedFloat
    \begin{subfigure}[t]{\textwidth}
        \includegraphics[width = \textwidth]{images/focus_evolution_side.png}
        \caption{The axial evolution of the light field for a focus to the side of the centre of the core of the fibre.}
        \label{sfig:focus_evo_side}
    \end{subfigure}
    \caption{Simulation results for light field evolution of a focus generated within the core of a fibre. The simulation shows two cases: in the centre and off to the side of the centre of the fibre core. Firstly, the left figures show the axial field evolution of the field intensity along the x-z plane, where z is the axial dimension. Secondly, the right figures show the transversal (x-y) field intensity at different depths (z = -150 $\mathrm{\mu m}$, -50 $\mathrm{\mu m}$, 0 $\mathrm{\mu m}$ and 50 $\mathrm{\mu m}$, where z = 0 $\mathrm{\mu m}$ is the focus plane).}
    \label{fig:focus_evo}
\end{sidewaysfigure}

However, for the grating to be used as a reflector, the grating must be short in order to localise the reflection axially. Thus, in figure \ref{fig:5um_sim} the same grating was simulated, but with its length set to 5 $\mu m$. It is clear that in this case the reflection amount is dependent on the focus location. In the case of focusing on the grating, the reflection has about $15\%$ of the input intensity and it is localised. The light which passes through the grating is slightly affected, the change will be required to be investigated if it affects the excitation focus quality. 

In the case of focusing away from the grating, the light passes almost unchanged. There is still a small amount of reflection, though 4 to 5 orders of magnitude smaller than the input focus. This is due to the focus quality being limited by the fibre modes, thus, some light still exists in the background even with a perfect focus. Thus, a short grating will be able to localise the reflection both in the transversal plane and axially.  

\begin{figure}
    \centering
    \begin{subfigure}[t]{0.75\textwidth}
        \includegraphics[width = \textwidth]{images/5um_ONgr_lin.png}
        \caption{The simulation result for a focus on the grating.}
        \label{sfig:5um_ON}
    \end{subfigure}
    \begin{subfigure}[t]{0.75\textwidth}
        \includegraphics[width = \textwidth]{images/5um_OFFgr_lin.png}
        \caption{The simulation result for a focus away from the grating.}
        \label{sfig:5um_OFF}
    \end{subfigure}
    \caption{Simulation results for focusing light on and away from the 5 $\mu m$ grating.}
    \label{fig:5um_sim}
\end{figure}

%%%
\subsubsection{Reflection spectrum simulation}
One of the most important properties of a Bragg grating is its ability to reflect light within a limited range of wavelengths. Thus, one important property to simulate is the spectrum of the reflected light.

The Bragg grating simulation approach described in this chapter follows a number of steps:
\begin{enumerate}
    \item Choose the fibre parameters.
    \item Compute the set of fibre modes for a specific wavelength of light.
    \item Choose grating parameters (excluding length).
    \item Compute the $\mathbf{M_B}$ matrix and diagonalise it.
    \item Choose a length for the grating and compute the scattering matrix defined in equation \ref{eq:scat_solution}.
    \item Choose a set of incoming fields (both forward and backwards propagating) and compute outgoing fields.
\end{enumerate}

In order for a spectrum to be simulated, as the wavelength of the light changes, all of the steps described above, except the first, must be passed through, as the number of modes and their profile will change. This implies that the process to simulate the spectrum is time consuming, depending on the available hardware. However, the advantage of the approach is that as long as some parameters are not changed, steps can be skipped. For example, changing the grating parameters does not require for the fibre modes to be generated again as long as the same fibre is used and the same wavelength range is desired. This is important if a specific grating performance is required and the profile needs to be optimised, while having the wavelength of light and the fibre parameters already be specified.

The first simulation was to compute the spectrum of the grating described in the previous section (\ref{ssec:bragg_focus_simulation}. The grating period was chosen for a wavelength of 515 nm. The result shown in figure \ref{fig:sample_16_spectrum}, which clearly illustrates that this choice is illustrated in the spectrum of the grating, with a bandwidth of approximately 3 nm.

\begin{sidewaysfigure}
    \centering
    \includegraphics[width = \textwidth]{images/Spectrum_on_grating_sample_16.png}
    \caption{The simulation spectrum result for the grating defined in section \ref{ssec:bragg_focus_simulation}. The input was normalised to a total power of 1 W.}
    \label{fig:sample_16_spectrum}
\end{sidewaysfigure}

Another simulation was done for a different Bragg grating. The profile of the grating was changed to be in the centre of the fibre core and the height of the grating increased to cover the entire facet, as shown in figure \ref{fig:grating_1550_profile}. The parameters of the new grating were:
\begin{itemize}
    \item Grating period of 1070.3 nm, for Bragg wavelength of 1550 nm (2nd order harmonic)
    \item Cross-sectional area: height extending more than the core ($>50 \,\mu m$), width $\approx 2 \,\mu m$
    \item Length of 5 mm
\end{itemize}

\begin{figure}
    \centering
    \includegraphics[width = 0.6\textwidth]{images/grating_profile_1550.png}
    \caption{Profile of the refractive index change of the Bragg grating defined for 1500 nm. Red circle indicates the core of the fibre. The simulation window is of size $128 \times 128$ pixels, indicated on the x and y axis.}
    \label{fig:grating_1550_profile}
\end{figure}

While the grating was kept the same, the two different cases for the input light field were chosen: a focus at the centre of the fibre core at the start of the grating and a focus at the centre of the core in the middle of the fibre grating. The results are shown in figure \ref{fig:grating_1550_spectrum}. As expected, the spectrum changes between the cases, but it is still around the 1500 nm chosen wavelength. The reason for the change is due to modes having a different coefficient. While in both cases it is a focus in the centre of the core, thus, the amplitude of the modal coefficients is identical, the phase components are not. This leads to a difference in the interaction with the grating. This is an important observation to keep in mind when optimising for a desired grating, that the input light is also of importance.

\begin{sidewaysfigure}
    \centering
    \begin{subfigure}[t]{\textwidth}
        \includegraphics[width = \textwidth]{images/reflected_power_focus_in_centre_at_fibre_input_1550.png}
        \caption{The simulation result for a focus in the centre of the core at the start of the grating.}
        \label{sfig:grating_1550_spectrum_focus_in}
    \end{subfigure}
\end{sidewaysfigure}
\begin{sidewaysfigure}
    \ContinuedFloat
    \begin{subfigure}[t]{\textwidth}
        \includegraphics[width = \textwidth]{images/reflected_power_focus_in_centre_at_grating_1550.png}
        \caption{The simulation result for a focus in the centre of the core in the middle of the grating.}
        \label{sfig:grating_1550_spectrum_focus_middle}
    \end{subfigure}
\end{sidewaysfigure}
\begin{figure}
    \ContinuedFloat
    \begin{subfigure}[t]{\textwidth}
        \includegraphics[width = \textwidth]{images/Grating_spectrum_from_paper_GRIN_Uniform_modes.png}
        \caption{The simulation result, figure 6-5 (b) in \cite{Bragg_master_thesis} for a graded index profile fibre. More info in the text body.}
        \label{sfig:grating_1550_grin_example}
    \end{subfigure}
    \caption{Simulated spectrum of the line Bragg grating defined for a wavelength of 1550 nm.}
    \label{fig:grating_1550_spectrum}
\end{figure}

To provide some confirmation to the truthfulness of the simulation, a plot from \cite{Bragg_master_thesis} is shown in figure \ref{sfig:grating_1550_grin_example}. The parameters are not identical, as in the cited thesis the grating covers the entire facet, is in a grade index fibre, the input is a field of modes of identical coefficients and the change in refractive index profile is in the order of $2 \times 10^{-4}$. However, a similar spectrum is achieved, the thesis specifically indicating that a large refractive index change leads to a wider and more uniform spectrum.

%%%%
\subsection{Conclusion}
In conclusion, this chapters presents a novel approach in simulating Bragg gratings. The advantage of the described approach is that it solves the mode-coupling equation simultaneously, without the resort to more numerical and iterative approaches, which can take a longer time. This allows for a quicker approach in simulating and optimising for the desired Bragg grating depending on the application.

In terms of the main objective of this thesis, the program allows for the computation of the optical field of a given Bragg grating inside a multimode optical fibre as well. This is important if Bragg gratings are to be used as reflectors for a correction approach to bent optical fibres, specifically within the scope of \textit{in vivo} imaging, where one of the ends of the fibre is not accessible during the correction procedure.

A number of simulation results have been shown, however, a more in depth analysis backed by experimental data is required in order to truly validate the approach.

\end{document}