\section{Introduction}\label{sec:introduction}

Music is often described as structured notes in time. Musical notations are systems that visually communicate this definition of music. The earliest known scores date back to 1250-1200 BC in Babylonia \cite{west1994babylonian}. Since then, many notation systems have emerged in different eras and different locations. Common Western Music Notation (CWMN) has become one of the most frequently used systems. This notation has evolved from the mensural music notation used before the seventeenth century. Current work in Optical Music Recognition focuses on the CWMN; nonetheless, studies are also carried out for old notations, including mensural, as shown in Table \ref{table:music-notation}.

\begin{table*}
    \centering

    \begin{tabular}{|p{9.1cm}|p{1.2cm}|p{0.8cm}|p{1.2cm}|p{1.7cm}|}
    \hline
    References & CWMN & Old & Typeset & Handwritten \\
    \hline
    \hline
    Fujinaga \cite{fujinaga1988optical},
    Co{\"u}asnon et al. \cite{couasnon1995using},
    Ng and Boyle \cite{ng1996recognition},
    Chen et al. \cite{chen2010math},
    Vidal \cite{vidal2012staffline},
    Bui et al. \cite{bui2014staff},
    Huang et al. \cite{huang2019state} & \tikz\fill[scale=0.4](0,.35) -- (.25,0) -- (1,.7) -- (.25,.15) -- cycle; & & \tikz\fill[scale=0.4](0,.35) -- (.25,0) -- (1,.7) -- (.25,.15) -- cycle; & \\
    \hline
    \hline
    Ng et al. \cite{ng1999embracing},
    Bainbridge and Bell \cite{Bainbridge2003AMN},
    Gocke \cite{gocke_building_nodate},
    Rebelo et al. \cite{shortest-path},
    Forn{\'e}s et al.\cite{fornes_use_2009},
    Pinto et al. \cite{pinto2011music},
    Haji{\v{c}} and Pecina \cite{hajivc2017muscima++},
    Roy et al. \cite{roy2017hmm},
    Pacha et al. \cite{pacha_handwritten_2018},
    Tuggener et al. \cite{watershed2018deep},
    Bar{\'o} et al. \cite{baro2019optical, baro2016towards} &
    \tikz\fill[scale=0.4](0,.35) -- (.25,0) -- (1,.7) -- (.25,.15) -- cycle; & & & \tikz\fill[scale=0.4](0,.35) -- (.25,0) -- (1,.7) -- (.25,.15) -- cycle; \\
    \hline
    Calvo-Zaragoza and Rizo \cite{calvo2018end},
    Wen et al. \cite{wen2015new},
    Pacha and Eidenberger \cite{pacha_self, pacha2017towards},
    Calvo-Zaragoza et al. \cite{calvo2018camera} & \tikz\fill[scale=0.4](0,.35) -- (.25,0) -- (1,.7) -- (.25,.15) -- cycle; & & \tikz\fill[scale=0.4](0,.35) -- (.25,0) -- (1,.7) -- (.25,.15) -- cycle; & \tikz\fill[scale=0.4](0,.35) -- (.25,0) -- (1,.7) -- (.25,.15) -- cycle; \\
    \hline
    Calvo-Zaragoza et al. \cite{calvo2015staff-avoiding, mensural-calvo2017handwritten},
    Huang et al. \cite{huang2015automatic},
    Tard{\'o}n et al. \cite{tardon2009optical} & & \tikz\fill[scale=0.4](0,.35) -- (.25,0) -- (1,.7) -- (.25,.15) -- cycle; & & \tikz\fill[scale=0.4](0,.35) -- (.25,0) -- (1,.7) -- (.25,.15) -- cycle; \\

    \hline
    \hline
    \end{tabular}
    \footnotesize
    \caption{Studies conducted in CWMN (Common Western Music Notation), and old notations (used before the CWMN, mostly mensural notations}
    \label{table:music-notation}
\end{table*}

Classifying music based on its difficulty is highly subjective. Nevertheless, Byrd and Simonsen \cite{byrd2015towards} in their attempt to have a standardised test-bed for OMR, name four categories based on the complexity of the score \cite{byrd2015towards} (see Figure \ref{fig:monophonic}):
\begin{enumerate}
 \item Monophonic: music in one staff with one note at a time;
 \item Polyphonic: multiple voices in one staff;
 \item Homophonic: multiple notes can occur at the same time to build up a chord, but only as a single voice;
 \item Pianoform: music in multiple staffs and multiple voices with significant structural interactions.
\end{enumerate}

\begin{figure*}[ht]
    \centerline{
    \includegraphics[width=11cm]{monophonic-pianophonic.png}}
    \caption{A visual representation of the four categories of music notations \cite{byrd2015towards, calvo-zaragoza_understanding_2019}}
    \label{fig:monophonic}
\end{figure*}

OMR has been researched for the last five decades; nonetheless, a unified definition of the problem is yet to emerge. However, Calvo-Zaragoza \cite{calvo-zaragoza_understanding_2019} offers the following definition of OMR.

\newtheorem{defn}{Definition}[section]
\begin{defn}
``Optical Music Recognition is a field of research that investigates how to computationally read music notation in documents.''
\end{defn}

The importance of OMR is evident both in the abundance of sheet music in archives and libraries, much of this is yet to be digitised, and in the common practice of musicians. Paper remains the first medium authors use to write music. By taking a picture of a score, OMR would enable us to later modify, play, add missing voices and share music using ubiquitous digital technologies. It also enables search capabilities, which are especially crucial for long pieces or large catalogues in music information retrieval and digital musicology. Other advantages of OMR include conversions to different sheet music formats (e. g. Braille music notation) and the ability to archive musical heritage \cite{jones2008optical}.

Fundamentally, OMR's goal is to interpret musical symbols from images of sheet music. The output would be a transcribed version of the sheet, which is also machine-readable, i.e., musical symbols can be interpreted and manipulated computationally. The usual output formats are MusicXML and MIDI. These formats will include musical attributes and information such as pitches, duration, dynamics and notes.

OMR has previously been referred to as Optical Character Recognition (OCR) for music. However, music scores carry information in a more complex structure, with ordered sequence of musical symbols together with their spatial relationships. In contrast, OCR deals with sequences of characters and words that are one-dimensional.

Recently, the success deep learning has had in improving text and speech recognition has triggered a paradigm shift in OMR as well. One of the most comprehensive reviews on OMR was written in 2012 by Rebelo et al. \cite{rebelo_optical_2012}. However, at that time, the field had not yet seen the emergence of deep learning approaches. This position paper aims to update on these approaches.

State of the art works in OMR perform well with digitally written monophonic music, but there is plenty of room for improvement when it comes to reading handwritten music and complex pianoform scores \cite{calvo2018end, wen2015new, pacha_self}. The difficulty thus increases with the complexity of the music notation.

\begin{figure*}[ht]
 \centerline{
 \includegraphics[width=13cm]{paper-pipeline.png}}
 \caption{Conventional OMR pipeline}
 \label{fig:OMR-framework}
\end{figure*}

\section{OMR Pipeline}

The standard OMR pipeline given by Rebelo et al. \cite{rebelo_optical_2012} is depicted in Figure \ref{fig:OMR-framework}:
\begin{enumerate}
 \item Image preprocessing;
 \item Music symbol recognition;
 \item Musical information reconstruction;
 \item Construction of a musical notation model.
\end{enumerate}
In the first stage, images of sheet music are subject to techniques such as noise removal, binarisation, de-skewing and blurring in order to make the rest of the OMR processes more robust. Subsequently, reference lengths, such as staff lines thickness and distances between them are calculated. Typically, the next stage is musical symbol recognition. This stage consists of staff line processing and musical symbol processing and ends with classification. Primitives of musical symbols  will be used in the third stage in order to reconstruct semantic meaning. Finally, all retrieved information should be embedded in an appropriate output file. A summary of these stages and the particular image processing and machine learning techniques employed in each stage are summarised in Table \ref{summary-of-studies}.

\begin{table*}
    \centering

    \begin{tabular}{|p{2.6cm}|p{13.5cm}|}
    \hline
    Stage & Related Work \\
    \hline
    \hline
    Image preprocessing &
    Fujinga \cite{fujinaga1988optical, fujinaga2004staff},
    Ng and Boyle \cite{ng1996recognition},
    Forn{\'e}s et al \cite{fornes2008writer, fornes_use_2009},
    Tard{\'o}n et al. \cite{tardon2009optical},
    Pinto et al. \cite{pinto2011music},
    Calvo-Zaragoza et al. \cite{calvo2015staff-avoiding},
    Huang et al. \cite{huang2015automatic},
    Wen et al. \cite{wen2015new},
    Ridler et al. \cite{ridler1978picture},
    Gocke \cite{gocke_building_nodate},
    Ballard \cite{ballard1981generalizing},
    Bainbridge and Bell \cite{Bainbridge2006},
    Cardoso et al. \cite{dos_santos_cardoso_staff_2009},
    Dalitz et al. \cite{dalitz_comparative_2008} \\
    \hline

    Symbol Recognition &
    Mahoney \cite{mahoney1982automatic},
    Prerau \cite{prerau-grammar},
    Tard{\'o}n et al. \cite{tardon2009optical},
    Pacha \cite{incremental-apacha},
    Rebelo et al. \cite{shortest-path},
    Ng and Boyle \cite{ng1996recognition},
    Choudhury et al. \cite{choudhury2000optical},
    Bainbridge and Bell \cite{Bainbridge2003AMN},
    Forn{\'e}s et al. \cite{fornes2008writer, fornes_use_2009},
    Huang et al. \cite{huang2015automatic}
    Fujinaga \cite{fujinaga1988optical},
    Wen et al. \cite{wen2015new},
    Pacha et al. \cite{pacha_handwritten_2018, pacha2017towards},
    Chen et al. \cite{chen2010math},
    Gocke \cite{gocke_building_nodate},
    Miyao and Nakano \cite{miyao1996note} \\
    \hline

    Musical Information Reconstruction &
    Prerau \cite{prerau-grammar},
    Pacha et al. \cite{pacha_learning_2019, mensural_pacha},
    Roy et al. \cite{roy2017hmm},
    Bainbridge and Bell \cite{Bainbridge2003AMN},
    Co{\"u}asnon et al. \cite{couasnon1995using},
    Ng and Boyle \cite{ng1996recognition},
    Bar{\'o} et al. \cite{baro2017optical, baro2019optical}
    Calvo-Zaragoza et al. \cite{calvo2018end, calvo2018camera, mensural-calvo2017handwritten} \\
    \hline

    Musical Notation Model &
    Droettboom et al. \cite{droettboom2002optical},
    Chen et al. \cite{chen2010math},
    Choudhury et al. \cite{choudhury2000optical}, Ng et al. \cite{ng1999embracing},
    Tard{\'o}n \cite{tardon2009optical}, Bainbridge and Bell \cite{Bainbridge2003AMN},
    Huang et al. \cite{huang2015automatic} \\

    \hline
    \hline
    \end{tabular}
    \footnotesize
    \caption{Summary of the studies carried in each of the OMR pipeline stages}
    \label{summary-of-studies}
\end{table*}

\section{Image preprocessing}

Image preprocessing is a fundamental step in many computer vision tasks. The primary outcome of this stage is an adjusted image that is easier to manipulate. Most common image manipulations include enhancement, de-skewing, blurring, noise removal and binarisation \cite{fujinaga1988optical, ng1996recognition, fornes2008writer, fornes_use_2009, tardon2009optical, pinto2011music, calvo2015staff-avoiding, huang2015automatic, wen2015new}. Image enhancement can include filters and adjusting the contrast or brightness for optimal object detection. De-skewing eliminates skewness and helps in obtaining a more appropriate view in the object detection stage. Most of the digital images during the acquisition, transmission or processing are subject to noise. Both colour and brightness contain signals that carry random noise. Depending on the features of the image, different types of filters are used to remove some of the noise.
During the process of binarisation, images are analysed to decide what is noise and what constitutes useful information for the task. Techniques to choose a binarisation threshold include global and adaptive methods. A global threshold is typically determined for the whole image, while for the adaptive threshold, local information in the image should be considered. Ng et al. \cite{ng1996recognition} adapt the global threshold proposed by Ridler and Calvard \cite{ridler1978picture}. While adaptive threshold is used in several more recent OMR studies too \cite{fornes_use_2009, fornes2008writer, gocke_building_nodate}.

Gocke's \cite{gocke_building_nodate} pipeline starts with a Gaussian filter, thenceforth a histogram on each colour channel is built. Then, the image is rotated to find the best angle that maximises horizontal projections. The image is then segmented into smaller 30x30 pixels, and a local threshold is found for each tile. Following threshold selection, all elements smaller than 4 pixels in diameter are removed, making the image clearer. The image is finally ready for staff-removal and symbol recognition. Local thresholding in this case yielded better results than the global one.

Similarly, in Forn{\'e}s et al. \cite{fornes_use_2009} binarisation is followed by de-skewing using the Hough Transform \cite{ballard1981generalizing}. A coarse approximation of the staff lines is obtained using median filters with horizontal masks to reconstruct the staff lines later. However, in this process, some residual colour information is retained, especially where the lines intersect with musical symbols, hence, some noise is still left.  This approach is not robust to damaged paper.

Pinto et al. \cite{pinto2011music} propose a content-aware binarisation method for music scores. The model captures content-related information during the process from a greyscale image. It also extracts the staff line thickness and the vertical line distance in staff to guide binarisation. This algorithm tries to find a threshold that maximises the extracted content information from images. However, the performance hugely depends on the document characteristics, limiting performance across different documents.

Calvo-Zaragoza and Gallego \cite{calvo2019selectional, gallego2017staff} propose using selectional auto-encoders \cite{masci2011stacked} to learn an end-to-end transformation for binarisation. The network activation nodes indicate the likelihood of whether pixels are foreground or background pixels. Ensuing training, documents are parsed through the model and binarised using an appropriate global threshold. This approach performs better than the conventional binarisation methods in some document types. Nonetheless, errors happen around foreground strokes and are emphasised along edges of the input windows, due to the lack of context in the neighbourhood.

\section{Music Symbol Recognition}

The next stage typically constitutes dealing with musical symbol recognition. Here, the three main steps are staff processing, isolating musical symbols and finally, classification. Usually, staff lines are first detected and then removed from the images. The model then isolates the remaining notations as primitive elements. These are later used to extract features and feed those features to train the classifier.

\subsection{Staff processing}

Staff lines are a set of five horizontal lines from one side of the music score to the other. Each line and gap represent a different pitch. For better object detection, the question of staff line removal has been of prime importance. Researchers take two different approaches; one is only detecting and isolating them, while the other approach goes one step further in removing them.

While in printed sheet music, staff lines are straight, parallel and horizontal, in handwritten scores, these lines might be tilted, curved and may not be parallel at all. These lines might also look curved or skewed depending on the image skew angle \cite{shortest-path} or the degradation of the paper. The model needs to separate staff lines from actual music objects. Since the lines overlap with musical objects, simply cutting and removing them degrades the notes and make them harder to recognise, further limiting performance.

Consequently, an increasing number of studies take the approach of removing the staff lines in a more intelligent fashion \cite{ng1996recognition, choudhury2000optical, blostein1992critical,Bainbridge2003AMN, fornes2008writer, fornes_use_2009, tardon2009optical, huang2015automatic, wen2015new}. In this section, we outline typical staff line processing approaches. Blostein and Baird \cite{blostein1992critical} suggests using horizontal projections of the black pixels and finding their maxima. The drawback is that the method only considers horizontal straight lines. In order to deal with non-horizontal, the process is followed with image rotations and choosing an angle with a higher maxima.

Rebelo et al. (2007) \cite{shortest-path} consider staff lines to be the shortest path between two horizontal page margins if those paths have black pixels throughout the entire path. The height between every two lines is first estimated and later used as a reference length for the following operations. Upon choosing an estimation, using the Dijkstra algorithm \cite{dijkstra1959note}, the shortest path between the leftmost pixel and the rightmost pixel is found. Their method is robust to lines with some curvature and discontinuity since it follows continuous paths connecting line ends from both sides. However, this algorithm may sometimes retain paths that do not follow the staff line. This happens when there is a higher density of beamed notes, and the estimated path follows the beams or when the staff lines are very curved.

Cardoso et al. \cite{dos_santos_cardoso_staff_2009} propose stable paths, considering the sheet music image as a graph. The staff lines in the graph are the less costly paths between the left and right margins. Subsequently, the model should differentiate between score pixels and staff line pixels. This model is robust to discontinuities, skewness, curvature in staff lines and one-pixel thin staff lines. Both the shortest path and stable paths give a similar false detection rate in test set of 32 ideal score images. This set is subject to different deformations, resulting in 2,688 total images. However, the stable path approach is five times faster. This technique is often used in the preprocessing stage \cite{mensural-calvo2017handwritten}.

Another study \cite{bui2014staff} uses stable paths approach to extract staff line skeletons. Then, the line adjacency graph (LAG) \cite{iliescu1996proposed} is used to cluster pixel runs generated from run-length encoding (RLE) of the image \cite{tsukiyama1986method}. The last step involves removing clusters lying on the staff line. This step has two passes; the first step estimates the height line for each staff by averaging the section height being cut with the staff lines. The second pass filters out the noise left from the last pass. This method takes a similar approach with \cite{carter1992automatic} grouping staff line pixels into segments.

Other studies follow the approach of keeping the staff lines during the next stages \cite{ng1999embracing, gocke_building_nodate, pugin2006optical, calvo2015staff-avoiding, roy2017hmm, calvo2016early}. They argue that the staff line removal task is very complex and often ends up being inaccurate and passes errors to the following stages. These studies usually detect and isolate staff lines ahead of object processing. Recent object detection studies show that removing staff lines does not add much improvement to this stage \cite{pacha_handwritten_2018}.

 \begin{figure*}[ht]
 \centerline{
 \includegraphics[width=13cm, ]{dl-framework.png}}
 \caption{Typical OMR pipeline using deep neural networks}
 \label{fig:dl-learning-framework}
\end{figure*}
A more recent work \cite{incremental-apacha} investigates how incremental learning can assist staff line detection using convolutional neural networks (CNNs) and human annotation. To begin with, a CNN model is fed a small amount of data with available annotations for training. Using this training, the model makes predictions on a larger dataset, and a human annotator rejects or accepts the predictions. The accepted predictions are added to the training dataset to repeat the process. This method enables the creation of a more extensive dataset. After four iterations, the dataset contains 70\% annotated scores of the original set. One drawback of incremental learning is that if the annotator accepts samples with imperfect annotations, the error accumulates in each iteration, introducing inaccuracy, while it also needs a human annotator. This yields similar results with \cite{dos_santos_cardoso_staff_2009, carter1992automatic}, however, different evaluation metrics are used.

Despite the substantial research effort put into staff line removal, it is still far from being accurate in handwritten sheet music. Handwritten scores exhibit a wide variety in line length and distance, thickness, curvatures of staff lines and also the quality of the image.

\subsection{Music symbol processing}

The next step after removing the staff lines is to isolate the musical symbols. Staff line removal will strongly affect this step as it can cause fragmentation in the parts where staff lines and musical objects are tangent to each other. One widely used approach is hierarchical decomposition \cite{rebelo_optical_2012}, where staff lines split a music sheet and then extract noteheads, rests stems and other notation elements \cite{choudhury2000optical, droettboom2002optical, gocke_building_nodate, miyao1996note, ng1996recognition}. Some approaches consider, for instance, a half-note instead of its primitives for the classification step. Mahoney \cite{mahoney1982automatic} uses descriptors to choose the matching candidate between a set of candidates of symbol types. Carter \cite{carter1989automatic} uses the line-adjacency graph (LAG) of an image for both removing the staff lines and providing a structural analysis of symbols. This technique helps in obtaining more consistent image sectioning, but it is limited to a small range of symbols as well as a potentially severe break-up of symbols.

Some studies skip segmentation and staff line removal \cite{pugin2006optical, calvo2016early, roy2017hmm} and use Hidden Markov Models (HMM). HMMs work on low-level features that are robust to poor quality images and can detect early topographic prints and handwritten pieces. Calvo-Zaragoza \cite{calvo2016early} split sheet music pages into staves following preprocessing. All staves are normalised and later represented as a sequence of feature vectors. This approach is very similar to \cite{pugin2006optical}, however, this study goes one step further and supports the HMM with a statistical N-gram model and achieve a 30\% error rate. This performance could be further improved if lyrics are removed, light equalisation is performed and data variations are statistically modelled.

\subsection{Music symbol classification}

After the segmentation of musical primitives, the subsequent process is classification. Objects are classified based on their shapes and similarities. However, since these objects are very often densely packed and overlapping their shapes can become very complex. Therefore, this step is very sensitive to all possible variations in music notations.
Fujinaga \cite{fujinaga1988optical} uses projection profiles for classification, Gocke \cite{gocke_building_nodate} uses template matching to classify the objects. Other methods used are support vector machines (SVMs), k-nearest neighbour (kNN), neural networks (NN) and hidden Markov models (HMM). A comparative study of the four methods \cite{rebelo2010optical}, finds SVM performs better than HMMs.

Considering the success of deep neural networks (DNN) in many machine learning tasks, recent studies take this approach in music object recognition and classification. A typical pipeline is shown in Figure~\ref{fig:dl-learning-framework}. These networks have many layers with activation functions employed before information propagates to the next layer. The deeper the model, the more complicated it gets and is able to detect hidden nonlinear relationships between the data, in this case, music objects. The problem with using DNNs in OMR is that they require a significant amount of labelled data for supervised training.

Object detection in images is a very active research field. Regional CNNs (R-CNNs), Faster R-CNN \cite{ren2015faster}, U-nets \cite{hajic2018towards}, deep watershed detectors \cite{watershed2018deep} and Single-shot detectors \cite{dai2016r, liu2016ssd} are among some of the approaches proposed recently. Pacha et al. \cite{pacha_handwritten_2018} use Faster R-CNN networks with pre-trained models fine-tuned with data from MUSCIMA++ (see Sect.~\ref{sec:datasets} for a summary of OMR datasets). They achieve a mean average precision of up to 80 \%. However, such performance is achieved with cropping the image into individual staff lines.

Tuggener et al. \cite{watershed2018deep} use deep watershed detectors in the whole image. It is faster than Faster R-CNN approach in image snippets, and it allows some shift in the data distribution. Nonetheless, it does not perform well on underrepresented classes.

Going further into the pipeline, we should be able to capture and reconstruct the right positions, relationships between notes, and relevant musical semantic information such as duration, onsets, pitch.

\section{Notation Reconstruction}

After classifying and recognising musical objects, the next block should extract musical semantics and structure. As mentioned earlier, OMR is two-dimensional, meaning that recognising the note sequence as well as their spatial relationships are essential. Hence, a model should identify the information about the spatial relationship between the recognised objects.
Ng et al. \cite{ng1999embracing} believe that domain knowledge is key to improving OMR tasks and especially music object recognition, similarly to a trained copyist or engraver, to decipher poorly written scores, building on the authors' previous research on printed scores \cite{ng1996recognition}. A multi-stage process is adopted, in which the first search is for essential features helping the interpretation of the score, verified by their mutual coherence, followed by a more intelligent search for more ambiguous features. Key and time signatures are detected after low-level processing and classification, using these global high-level features to test the earlier results.

Ng and Boyle \cite{ng1996recognition} base their study on three assumptions: i) foreknowing the time signature, ii) key signature, and iii) that the set of the primitive feature set under examination is limited to ten. The first and second assumptions are overcome by geometrically predicting a limited symbol set such as numbers, flats and sharps. The input image goes through binarisation using a threshold, image rotation for de-skew, then the staff lines are detected and erased. Now the image has blocks of pixels, music object primitives and groups of primitives. Further segmentation based on some rules is needed for a group of primitives. After the segmentation process, a classifier uses only the width and the height of the bounding box for recognition based on a sampled training set. The recognised primitives are grouped to reconstruct their semantic meaning. The reconstruction consists of overlaying an ellipse and counting the number of foreground pixels, finding the pitch, search the neighbourhood for other features that might belong to the object and identifying the possible accidents using a nearest neighbourhood (NN) classifier. Music knowledge related to bars, time, and key signatures is applied at this stage. During segmentation, the process relies on straight edges of the objects, therefore is not robust to handwritten scores. The method fails if the symbols are skewed, for instance, when a stem is not perpendicular to a stave line.

Similar to the method mentioned above, another approach is formalising musical knowledge and/or encoding knowledge into grammar rules that explain, for instance, how primitives are to be processed or how graphical shapes are to be segmented \cite{Bainbridge2003AMN, couasnon1995using}.

Prerau\cite{prerau-grammar} proposes two levels of grammar. One being notational grammar while the other is a higher-level grammar for music. The first allows the recognition of symbol relationships, the second deals with larger music units.
Many other techniques use musical rules to create grammar rules for OMR. Such rules can be exemplified as \cite{rebelo_optical_2012}:
\begin{itemize}
\item An accidental is placed before a notehead and at the same height;
\item A dot is placed after or above a notehead in a variable distance;
\item Between any pair of symbols: they cannot overlap.
\end{itemize}

The issue with music rules and heuristics is that these rules are very often violated, especially in handwritten music. Furthermore, it is challenging to create rules for many different variations and notations with a high level of complexity. As a result, this approach would not perform well with both typeset and handwritten complex notations, and it is difficult to scale to a broad range of notation and engraving styles.

Pacha et al. \cite{pacha_learning_2019} propose using graphs to move towards a universal music representation. Considering that in music notations, the relationship between primitives contains the semantic meaning of each primitive; they suggest that OMR should employ a notation assembly stage to represent this relationship. Instead of using grammar and rules mentioned earlier, they use a machine learning approach to assemble a set of detected primitives. The assembly is similar to a graph containing syntactic relationships among primitives capturing the symbol configuration. The robustness of the model regarding variations in bounding boxes leaves room for improvement and so does the notation assembly stage, due to the lack of broader hypotheses on the detected objects.

Bar{\'o} et al. \cite{baro2017optical} consider monophonic scores as sequences and use Long Short-Term Memory (LSTM) Recurrent Neural Networks (RNNs) for reading such sequences to retrieve pitch and duration. For evaluation they use Symbol Error Rate (SER) defined as the minimum number of edit operation to convert an array to another. This approach shows to work well with simple scores such as monophonic scores, but fundamental remodelling is needed for more complex scores \cite{baro2017optical}.

This stage is concerned with reconstructing relationships from the detected musical objects. A challenge in this stage is to model a musical output representation that encodes sheet music both a similar rendering of the original image and the semantics (e.g. onsets, duration, pitch).

\vspace{-5pt}
\section{Music Notation Encoding}

The output from the previous steps is used to construct a semantic model or data model. This model should represent a re-encoding of the score in the input. The output model should be expressible in a machine-readable format. Usual OMR output formats include MIDI, MusicXML, MEI, NIFF, Finale, and in some software, the music is even rendered into WAVE files. Musical Instrument Digital Interface (MIDI) \cite{rothstein1995midi} is an interchange medium between the computer and digital instruments. At the basic level, MIDI includes the temporal position when a note starts, stops, how loud the note is, the pitch of the note, instrument and channel. The main drawback of MIDI is that it cannot represent the relationships between musical symbols, or produce a re-encoded structured file, limiting the output to replayability only.

Notable formats that allow a structured encoding and storing notations include MusicXML \cite{good2001musicxml,good2006lessons} and MEI \cite{roland2002music, hankinson2011music}. Both allow further editing in a music notation software. MusicXML is more focused on encoding notation layout. It is designed for archiving and for sharing sheet music between applications. There is ongoing research in the W3C Music Notation Community Group on improving MusicXML format to handle more specific tasks and applications.

The Music Encoding Initiative (MEI) \cite{roland2002music} claims to be comprehensive, declarative, explicit and hierarchical. MEI has not been widely used as the final output of OMR systems yet. However, based on the characteristics mentioned above, MEI is able to capture and retain musical semantics better, e.g. relationships between voices, which may benefit music engraving.

There is also work converting OMR output into Semantic Web formats. Jones et. al.~\cite{jones2017musicowl} propose the use of Linked Data to annotate and improve discovery of music scores using the Resource Description Framework (RDF). The captured information is limited to the number of voices, movements and melodies. Further extensions are needed to store more sophisticated music semantics that support harmony or melody analysis. Nevertheless, the use of Linked Data compatible formats may benefit OMR applications in multiple ways. Linking scores to other music related data on the Web~\cite{fazekas2010an} or even features of the audio of a performance~\cite{allik2016ontology} could support interactive applications such as score following or large catalogue navigation~\cite{turchet2019cloud}. The ontologies governing these formats may be used to encode musical or engraving rules to complement probabilistic inference in machine learning models.

To decide which of the encodings to use, we have to think of what an application may require. Using the knowledge obtained in the previous steps and from different studies would assist this stage in its standardisation. Currently there is little research in OMR dealing with encoding, however, many works in other fields focus on encoding formats that better represent music and its structure.

\vspace{-10pt}
\section{Datasets}
\label{sec:datasets}
Depending on the OMR task to be performed and the nature of the application, different datasets may be suitable. Existing datasets contain handwritten or copyright-free printed music sheets in mensural or CWMN notations.
Calvo-Zaragoza et al. \cite{calvo2014recognition} introduced a new dataset called HOMUS (Handwritten Online Musical Symbols). This contains 15200 samples of 32 types of musical symbols from 100 different musicians.
Universal Music Symbol Collection is a dataset of 90000 tiny handwritten and typeset music symbols from 79 classes that can be used to train classifiers.

As for staff line removal, a commonly used dataset is CVC-MUSCIMA \cite{fornes2012cvc}. It contains 1000 music sheets written by 50 different musicians. Each musician was asked to transcribe the same given 20 pages of music using the same pen and same style of sheet music paper. These pages include monophonic and polyphonic music, consisting of scores for solo instruments and music scores for choir and orchestra.

A derived version of CVC-MUSCIMA dataset is MUSCIMA++ \cite{hajivc2017muscima++}. This dataset is more suitable for musical symbol detection. It has 91255 symbols with both notation primitives and higher-level notation objects, key signatures or time signatures. Notes are captured using the annotated relationships of the primitives, having this way both low and high-level symbols.
DeepScores is a collection that contains 300k annotated images of written music mainly for object classification, detection, and segmentation \cite{tuggener2018deepscores}. This dataset has large images containing tiny objects.

There are also datasets for an end-to-end recognition such as the Printed Images of Music Staves (PrIMuS) \cite{calvo2018end}, or the extended version of this with distorted images to simulate imperfections Camera-PrIMuS \cite{calvo2018camera}. These datasets have 87678 real-music scripts in five different formats: PNG, MIDI, MEI, semantic and agnostic encoding which is a sequence that contains the graphical symbols and their positions without any musical meaning.

Given that the performance of the deep learning methods usually depends on the amount of the data the model is fed, for future work, we propose creating a universal dataset that facilitates the intermediate stages but also an end-to-end system. We want to start by generating music files using a music notation software such as Dorico \cite{dorico} or Rosegarden \cite{rosegarden}. This work will be harmonized with the before-mentioned MUSCIMA++ and DeepScores datasets.

\vspace{-10pt}
\section{Open Issues and Conclusions}

Low-quality images of sheet music, complex scores, handwritten music and alternate notations are still challenging for OMR, while most of the work focuses on monophonic scores. CWMN notation is highly complex, having dense scores, overlapping symbols, structural complexity, semantic rules that are sometimes violated. For a deep learning approach, in particular, class imbalance is one of the most significant issues; some note types are persistent while some others are rare. An further open issue is the lack of a large labelled dataset with a broad variety of image quality and balanced classes \cite{novotny2015introduction}.

We can observe a shift in OMR from using conventional image processing and object detection to using neural networks, as shown in Figure \ref{fig:dl-learning-framework}. Recently published papers take novel approaches and use deep learning methods in all stages of the OMR pipeline. These stages are not necessarily in the order presented above or exhibit all the steps described.

Despite the introduction of deep learning, the field leaves space for improvement in all stages of the pipeline. New opportunities include creating more diverse and better balanced datasets, improving the detection of music objects and staff lines, the reconstruction of semantic meaning, and, perhaps most importantly, standardising the evaluation metrics and the output of the pipeline. A possible final goal is end-to-end learning that would not need intermediate steps. Neural networks are already applied to problems like text and speech recognition and machine translation in this manner. However, these systems are still not adapted to a two-dimensional output sequence such as music \cite{calvo-zaragoza_understanding_2019}.

This paper summarised seminal and influential studies conducted in the field of OMR. We discussed different methods and approaches in prominent stages of the OMR pipeline. Our review aims to identify important older works and current state-of-the-art approaches, which can be used as a reference by researchers to begin further work in OMR. It also represents a position in several aspects of the field, including the need for incorporating more prior knowledge, theory and musical information in the processing pipeline, the need for finding new methods to incorporate these priors into statistical learning models such as deep neural networks and a need for more standardisation in OMR evaluation.

\begin{acknowledgments}
 The authors acknowledge the support of the AI and Music CDT, funded by UKRI and EPSRC under grant agreement no. EP/S022694/1 and our industry partner Steinberg Media Technologies GmbH for their continuous support.
\end{acknowledgments}


\section{Introduction}

\subsection{Statement of the Problem}
Let $T>0$ be fixed, and let $(M,g)$ be a connected, compact, smooth Riemannian manifold of dimension $m\geq2$ with boundary $\d M$. In what follows, we shall additionally assume that $(M,g)$ is non-trapping. Consider a trivial Hermitian vector bundle $E=M\times\C^n$ equipped with the Hermitian inner product $\pair{\cdot,\cdot}_E$.\\

We say that a connection $\nabla:C^\infty(M;E)\rightarrow C^\infty(M;E\otimes T^\ast M)$ is compatible with the Hermitian structure of $E$ if for any sections $u,v\in C^\infty(M;E)$ it holds that\begin{equation}\label{HIP}d\pair{u,v}_E=\pair{\nabla u,v}_E+\pair{u,\nabla v}_E,\end{equation}where both sides of the above are regarded as sections of the cotangent bundle.\\

Such a connection has the form $\nabla=d+A$, where $A=A_idx^i$ and each $A_i(x)$ is given by an $n\times n$ skew-Hermitian matrix. In what follows, we allow the connection form $A$ to also depend smoothly on time, and write $\nabla^A:C^\infty((0,T)\times M;E)\rightarrow C^\infty((0,T)\times M;E\otimes T^\ast M)$ for the time-dependent connection corresponding to the connection form $A$. In other words, $\nabla^{A(t)}$ is a connection on $C^\infty(M;E)$ for each $t\in[0,T]$, and each $\nabla^{A(t)}$ is compatible with the Hermitian metric on $E$.\\

We can define a natural $L^2$-inner product on $C^\infty((0,T)\times M;E)$ via \[\pair{u,v}_{L^2((0,T)\times M;E)}=\int_0^T\int_M\pair{u,v}_EdVdt,\]where $dV$ denotes the usual Riemannian volume measure of $(M,g)$. We can similarly define a natural $L^2$-inner product on $C^\infty((0,T)\times M;E\otimes T^\ast M)$. For $E$-valued $1$-forms $\alpha=\alpha_jdx^j$ and $\beta=\beta_jdx^j$, we set\[\pair{\alpha,\beta}_{L^2((0,T)\times M;E\otimes T^\ast M)}=\int_0^T\int_Mg^{ij}\pair{\alpha_i,\beta_j}_EdVdt,\]where $g^{ij}$ denotes the inverse of the metric tensor.\\

We let $(\nabla^A)^\ast$ denote the adjoint of $\nabla^A$ with respect to the above inner products. We can then define the connection Laplacian $\Delta_A=-(\nabla^A)^\ast\nabla^A$, which corresponds to the connection form $A$.\\

We can compute local expressions for $(\nabla^A)^\ast$ and $\Delta_A$. Consider a section $u\in C^\infty(M;E)$ and an $E$-valued $1$-form $\beta=\beta_jdx^j$ supported on a local trivialisation. Since $A$ is skew-Hermitian, it holds that\[\pair{Au,\beta}_{L^2((0,T)\times M;E\otimes T^\ast M)}=\int_0^T\int_Mg^{ij}\pair{A_iu,\beta_j}_EdVdt=-\int_0^T\int_M\pair{u,g^{ij}A_i\beta_j}_EdVdt.\]

Letting $(A,\beta)_g=g^{ij}A_i\beta_j$, we see that $(\nabla^A)^\ast=d^\ast-(A,\cdot)_g$. Therefore, we have\[\Delta_Au=-d^\ast du-d^\ast(Au)+(A,du)_g+(A,Au)_g.\]Recall that $d^\ast\alpha=-\abs{g}^{-\frac{1}{2}}\d_i(\abs{g}^{\frac{1}{2}}g^{ij}\alpha_j)$, where $\abs{g}$ is the determinant of the metric tensor. Hence it holds that $d^\ast(Au)=(d^\ast A)u-(A,du)_g$. Thus, we conclude that\begin{equation}\Delta_Au=-d^\ast du+2(A,du)_g-(d^\ast A)u+(A,Au)_g.\label{CLL}\end{equation}

Lastly, we say that a section $V\in C^\infty((0,T)\times M;\C^{n\times n})$ is a potential if $V$ is Hermitian or, equivalently, for any sections $u,v\in C^\infty((0,T)\times M;E)$ it holds that\[\pair{Vu,v}_E=\pair{u,Vv}_E.\]

Let $\Delta_A$ and $V$ be as above and consider the following initial and boundary value problem for sections $u\in C^\infty((0,T)\times M;E)$.\begin{equation}\begin{split}\label{1.1}i\d_tu(t,x)+\Delta_{A(t)}u(t,x)+V(t,x)u(t,x)&=0\textrm{ in }(0,T)\times M,\\u(t,x)&=f\textrm{ on }(0,T)\times\d M,\\u(0,x)&=0\textrm{ in }M,\end{split}\end{equation}where the inhomogeneous Dirichlet data is given by $f\in C^\infty((0,T)\times\d M;E)$ satisfying $f\vert_{t=0}=\d_tf\vert_{t=0}=0$. We can then define the associated Dirichlet-to-Neuman map via\[\Lambda_{A,V}f={\nabla^A_\nu u}\Big\vert_{(0,T)\times\d M}\ ,\]where $\nu$ denotes the outward pointing unit normal vector field on $\d M$.\\

There is a natural gauge group associated with the equation above. Let $G:(0,T)\times M\rightarrow U(n)$ be a smooth map such that $G(t)\vert_{\d M}=\Id$, and choose $A_2=G^{-1}A_1G+G^{-1}dG$ and $V_2=G^{-1}V_1G+iG^{-1}\d_tG$. It then holds that $\nabla^{A_2}=G^{-1}\nabla^{A_1}G$, and hence that $\Delta_{A_2}=G^{-1}\Delta_{A_1}G$. We observe that if $u$ solves (\ref{1.1}) with $A=A_2$ and $V=V_2$, then $Gu$ solves the equation with $A=A_1$ and $V=V_1$, since\[(i\d_t+\Delta_{A_1}+V_1)Gu=G(i\d_t+\Delta_{A_2}+V_2)u=0.\]Furthermore, we observe that when the pairs $(A_1,V_1)$ and $(A_2,V_2)$ are as above, it holds that $\Lambda_{A_1,V_1}=\Lambda_{A_2,V_2}$. Therefore, we can only hope to recover the pair $(A,V)$ up to a gauge transform. The aim of the present work is to establish unique recovery of the connection form and potential from the knowledge of the Dirichlet-to-Neumann map, modulo gauge invariance.\\

\subsection{History of the Problem}

Literature dealing with the recovery of space- and time-dependent potentials of the dynamic Schr\"odinger equation is limited, even in the scalar case. For Euclidean domains, it was shown in \cite{11} that the time-dependent electromagnetic potentials are uniquely determined by the Dirichlet-to-Neumann map. Logarithmic-stable determination was shown for the electric potential in \cite{9}, and this result was extended to the full electromagnetic potential in \cite{7}, provided that the time-independent part of the magnetic potential is sufficiently small. Indeed, it was only recently shown in \cite{22} that time-dependent electromagnetic potentials in a Euclidean domain can be H\"older-stably recovered from the knowledge of the D-to-N map. We also mention here the recent work of \cite{BbF}, which establishes logarithmic and double-logarithmic stability estimates for the same problem with partial data.\\

In the Riemannian setting, \cite{3} and \cite{4} establish, respectively, H\"older-stable recovery of a time-independent magnetic and electric potential of the dynamic Schr\"odinger equation on a simple manifold. These results were extended to simultaneous recovery of both electromagnetic potentials in \cite{2}. In the case of time-dependent potentials in the Riemannian context, the only result is that of \cite{KT}, establishing, on a simple manifold, the H\"older-stable recovery of both potentials from the knowledge of the Dirichlet-to-Neumann map.\\

In the case of the vector-valued dynamic Schr\"odinger equation, there are, to the best of the author's knowledge, no results establishing unique recovery even for time-independent coefficients. However, such results do exist for the related case of the stationary Schr\"odinger equation. In particular, for the stationary Schr\"odinger equation on a trivial vector bundle over a Euclidean domain, \cite{Esk1} establishes unique recovery of a connection form and potential from the knowledge of the Dirichlet-to-Neumann map. Additionally, for the stationary Schr\"odinger equation on a Hermitian vector bundle over a two-dimensional Riemann surface, \cite{AGTU} uniquely recovers the coefficients from Cauchy data at the boundary.\\

Let us also mention the paper \cite{MC1}, where it is conjectured that the Dirichlet-to-Neumann maps for two connection Laplacians coincide in the case of the stationary Schr\"odinger equation if and only if the associated connection forms are gauge equivalent. The present work solves this conjecture in the non-stationary case. More precisely, we show that the Dirichlet-to-Neumann map uniquely determines, up to gauge invariance, the space- and time-dependent connection form and potential appearing in the dynamic Schr\"odinger equation on a trivial vector bundle over a Riemannian manifold, provided that the manifold in question satisfies certain geometric conditions.\\

Finally, we note the works \cite{MC2}, \cite{KOP}, and \cite{CLOP}, where various inverse problems for partial differential equations involving connections are considered. In particular, \cite{MC2} establishes unique recovery of a time-independent unitary Yang-Mills connection on a Hermitian vector bundle, as well as recovering the bundle structure. Similarly, in the case of the wave equation, \cite{KOP} uses techniques from the boundary control method to reconstruct a Riemannian manifold and Hermitian vector bundle with a time-independent compatible connection from the knowledge of the associated hyperbolic Dirichlet-to-Neumann map, and the work \cite{CLOP} considers the non-linear inverse problem of recovering a time-independent connection from a cubic wave equation on a Hermitian vector bundle over the Minkowski space $\R^{1+3}$.\\

Lastly, let us mention the recent work of \cite{MV}, where the authors recover a time-dependent potential of the wave equation on a trivial vector bundle over a Euclidean domain from knowledge of the input-output operator on the partial boundary.\\

\subsection{Geodesics and Parallel Transport}

Let us assume that $(M,g)$ is non-trapping, which is to say that every geodesic in $M$ reaches the boundary in finite time. We now take a moment to recall certain key facts relating to the geodesics in $M$.\\

Given $x\in M$ and $\theta\in T_xM$, we denote by $\gamma_{x,\theta}$ the geodesic with initial point $x$ and initial direction $\theta$. We define the sphere bundle of $M$ via\[SM=\{(x,\theta)\in TM:\abs{\theta}_g=1\}.\]

Likewise, we define the submanifold of inner vectors $\d_+SM$ via\[\d_+SM=\{(x,\theta)\in SM:x\in\d M,\ \pair{\theta,\nu(x)}_{g(x)}<0\},\] where $\nu(x)$ is the outward pointing unit normal vector at $x\in\d M$, and we define the submanifold of outer vectors $\d_-SM$ via\[\d_-SM=\{(x,\theta)\in SM:x\in\d M,\ \pair{\theta,\nu(x)}_{g(x)}>0\}.\] Then, for $\gamma_{x,\theta}$ such that $(x,\theta)\in\d_+SM$, we can define the exit-time $\rho_+(x,\theta)$ of the geodesic in $M$ by\[\rho_+(x,\theta)=\min\{s>0:\gamma_{x,\theta}(s)\in\d M\}.\]

Given the above, we further recall the parallel transport equations associated to a connection $\nabla^A$. For any geodesic $\gamma_{x,\theta}$ with $(x,\theta)\in\d_+SM$, and any initial vector $w\in E_x$, we consider the parallel transport equation along $\gamma_{x,\theta}$, given by\begin{equation}\begin{split}\label{PT1}\Big[\d_r+A\big(\gamma'_{x,\theta}(r)\big)\Big]W=&\ 0\\W(0)=&\ w.\end{split}\end{equation}The transport of $w\in E_x$ along $\gamma_{x,\theta}$ is thus given by $W(r)$.\\

It is frequently helpful to consider the fundamental matrix solution $U_A:\big[0,\rho_+(x,\theta)\big]\rightarrow U(n)$ of the parallel transport equation:\begin{equation}\begin{split}\label{PT2}\big[\d_r+A\big(\gamma'_{x,\theta}(r)\big)\big]U_A=&\ 0\\U_A(0)=&\ \Id.\end{split}\end{equation}

It is clear from the above that the transport of $w\in E_x$ along $\gamma_{x,\theta}$ is given by $W(r)=U_A(r)\cdot w$.\\

Given $(x,\theta)\in\d_+SM$, we define the scattering data for the connection as the map $C_A:\d_+SM\rightarrow U(n)$ given by\begin{equation}C_A(x,\theta):=U_A\big(\rho_+(x,\theta)\big).\label{ScatteringData}\end{equation}

Since elements of the gauge group must satisfy $G(t)\vert_{\d M}=\Id$, we note that the scattering data $C_A$ is gauge invariant. We are now in a position to state the main results of the present work.\\

\subsection{Main Results}
\begin{theorem}\label{t1}Suppose that for $j=1,2$, we have connection forms $A_j\in C^\infty((0,T)\times M;\C^{n\times n}\otimes T^\ast M)$, and potentials $V_j\in C^\infty((0,T)\times M;\C^{n\times n})$. Then $\Lambda_{A_1,V_1}=\Lambda_{A_2,V_2}$ implies that $C_{A_1}=C_{A_2}$.
\end{theorem}

\begin{theorem}\label{t2}
Assume the conditions of Theorem \ref{t1} hold, and assume further that $M$ is either i) $2$-dimensional and simple, or ii) of dimension $m\geq3$ with strictly convex boundary, and admits a smooth strictly convex function. Then $(A_1,V_1)$ is gauge equivalent to $(A_2,V_2)$.
\end{theorem}

These results are, as far as the author is aware, the first dealing with the recovery of coefficients appearing in the dynamic Schr\"odinger equation on a vector bundle. In fact, the above results are the first showing recovery of time-dependent coefficients of any linear second-order partial differential equation with variable coefficients of leading order, in the vector-valued case. The proof of these results relies on the construction of Gaussian beam solutions which allow recovery of the scattering data corresponding to the connection form. This data is then used to recover the connection form and potential of the Schr\"odinger equation via the inversion of attenuated ray-transforms with matrix-weights corresponding to the connection form and potential we wish to recover. This last step relies on the results of \cite{C3} and \cite {C1}, which guarantee that the appropriate attenuated ray-transform is invertible when the base manifold is either i) two-dimensional and simple or ii) of higher dimension with strictly convex boundary and admits a smooth strictly convex function.\\

Here follows an outline of the present work. In section 2, we give some regularity results for the forward problem and for the Neumann trace. In section 3 we construct special Gaussian beam solutions for the Schr\"odinger equation. The proofs of Theorems \ref{t1} and \ref{t2} are given in sections 4 and 5 respectively.\\

\section{The Forward Problem}

Before we proceed, let us define for all $r,s\in(0,\infty)$ and for $X=M$ or $X=\d M$ the energy spaces $H^{r,s}((0,T)\times X;E)=H^r(0,T;L^2(X;E))\cap L^2(0,T;H^s(X;E))$, together with the associated norm\[\norm{u}^2_{H^{r,s}((0,T)\times X;E)}=\norm{u}^2_{H^r(0,T;L^2(X;E))}+\norm{u}^2_{L^2(0,T;H^s(X;E))}.\]

The main result of this section is the following well-posedness result for the Dirichlet-to-Neumann map:

\begin{proposition}\label{edit_p1}
The IBVP (\ref{1.1}) has a unique solution $u\in C^\infty((0,T)\times M;E)$. Further, there exists a constant $C>0$ such that the associated Dirichlet-to-Neumann map satisfies the estimate\begin{equation}\label{edit_DN}\norm{\Lambda_{A,V}f}_{L^2((0,T)\times\d M;E)}\leq C\norm{f}_{H^{\frac{9}{4},\frac{3}{2}}((0,T)\times\d M;E)}.\end{equation}
\end{proposition}

The proof of the above result essentially reduces to proving a suitable energy estimate for the source problem for the Schr\"odinger equation. Thus, for $F\in C^\infty((0,T)\times M;E)$ we consider the solution of the source problem\begin{equation}\begin{split}\label{s2.1}\big(i\d_t+\Delta_A+V\big)u&=F(t,x)\textrm{ in }(0,T)\times M,\\u(t,x)&=0\textrm{ on }(0,T)\times\d M,\\u(0,x)&=0\textrm{ in }M.\end{split}\end{equation}

\begin{proposition}\label{edit_p2}
The source problem (\ref{s2.1}) satisfies the energy estimates\begin{align*}\norm{u}_{L^\infty(0,T;L^2(M;E))}\leq \norm{F}_{L^2((0,T)\times M;E)} \\ \norm{u}_{L^\infty(0,T;H^1(M;E))}\leq \norm{F}_{H^{1,0}((0,T)\times M;E)} \\ \norm{\d_tu}_{L^\infty(0,T;L^2(M;E))}\leq \norm{F}_{H^{1,0}((0,T)\times M;E)} \\ \norm{u}_{H^{1,2}((0,T)\times M;E)}\leq C\norm{F}_{H^{1,0}((0,T)\times M;E)}.\end{align*}
\end{proposition}

\begin{proof}
By taking the inner product of (\ref{s2.1}) with $u$ and integrating by parts, we deduce that\begin{equation}\label{s2.2}i\int_M\pair{\d_tu,u}_EdV_g-\int_M\pair{\nabla^Au,\nabla^Au}_EdV_g+\int_M\pair{Vu,u}_EdV_g=\int_M\pair{F,u}_EdV_g\end{equation}Taking the imaginary part of (\ref{s2.2}) yields\[\frac{d}{dt}\Big(\norm{u(t)}_{L^2(M;E)}^2\Big)\leq C\Big(\norm{F(t,\cdot)}_{L^2(M;E)}\norm{u(t)}_{L^2{(M;E)}}+\norm{u(t)}_{L^2(M;E)}^2\Big).\]Then Gr\"onwall's inequality tells us that\begin{equation}\label{s2.3}\norm{u}_{L^\infty(0,T;L^2(M;E))}\leq C\norm{F}_{L^2((0,T)\times M;E)}.\end{equation}

On the other hand, taking the inner product of (\ref{s2.1}) with $\d_tu$, we can integrate by parts to deduce that\begin{equation}\label{s2.4}i\int_M\pair{\d_tu,\d_tu}_EdV_g-\int_M\pair{\nabla^Au,\nabla^A\d_tu}_EdV_g+\int_M\pair{Vu,\d_tu}_EdV_g=\int_M\pair{F,\d_tu}_EdV_g.\end{equation}Then, by setting\begin{gather*}\alpha(t;u,v)=\int_M\pair{\nabla^Au,\nabla^Av}_EdV_g-\int_M\pair{Vu,v}_EdV_g\end{gather*}and\begin{gather*}\alpha'(t;u,u)=\int_M\pair{(\d_tA)u,\nabla^Au}_EdV_g+\int_M\pair{\nabla^Au,(\d_tA)u}_EdV_g-\int_M\pair{(\d_tV)u,u}_EdV_g,\end{gather*}we can take the real part of (\ref{s2.4}) to conclude that\[\frac{d}{ds}\alpha(s;u,u)=\alpha'(s;u,u)-2\,\Real\int_M\pair{F(s,\cdot),\d_tu(s)}_EdV_g.\]

By integrating, we can rewrite the above as\begin{equation}\label{s2.5}\alpha(t;u,u)=\int_0^t\alpha'(s;u,u)ds-2\,\Real\int_0^t\pair{F(s,\cdot),\d_tu(s)}_{L^2(M;E)}ds,\end{equation}and since $\int_0^t\pair{F(s,\cdot),\d_tu(s)}_{L^2(M;E)}ds=\pair{F(t,\cdot),u(t)}_{L^2(M;E)}-\int_0^t\pair{\d_tF(s,\cdot),u(s)}_{L^2(M;E)}ds$, we can rewrite the identity (\ref{s2.5}) in the form\begin{equation}\label{s2.6}\alpha(t;u,u)\leq\int_0^t\alpha'(s;u,u)ds+2\norm{F(t)}_{L^2(M;E)}\norm{u(t)}_{L^2(M;E)}+2\int_0^t\norm{\d_tF(s,\cdot)}_{L^2(M;E)}\norm{u(s)}_{L^2(M;E)}ds.\end{equation}

Further, by expanding the inner products appearing in $\alpha(t;u,u)$ and using Cauchy's inequality with $\varepsilon=\frac{1}{2}$, we can deduce that\[\begin{split}\frac{1}{2}\norm{\nabla u(t)}_{L^2(M;E\otimes T^\ast M)}^2\leq\alpha(t;u,u)+\norm{V(t)}_{L^\infty(\C^{n \times n}\otimes T^\ast M)}\norm{u(t)}_{L^2(M;E)}^2\\+\norm{A(t)}_{L^\infty(\C^{n \times n}\otimes T^\ast M)}^2\norm{u(t)}_{L^2(M;E)}^2\end{split}\]whence\[\alpha(t;u,u)+\lambda\norm{u(t)}_{L^2(M;E)}^2\geq\frac{1}{2}\norm{u(t)}_{H^1(M;E)}^2,\] for $\lambda=\frac{1}{2}+\norm{V(t)}_{L^\infty(M;\C^{n\times n}\otimes T^\ast M)}+\norm{A(t)}_{L^\infty(M;\C^{n\times n}\otimes T^\ast M)}^2$. Then, combining the above with identity (\ref{s2.6}) and the definition of $\alpha'(t;u,v)$, we may deduce that $\norm{u(t)}_{H^1(M;E)}^2$ is bounded above by\begin{equation}\label{edit_2.8}2\lambda\norm{u(t)}_{L^2(M;E)}^2+4\norm{F(t)}_{L^2(M;E)}\norm{u(t)}_{L^2(M;E)}+C\Big(\int_0^t\norm{u(s)}_{H^1(M;E)}^2+\norm{\d_tF(s,\cdot)}_{L^2(M;E)}^2ds\Big).\end{equation}Let us briefly recall the inequality\[\norm{F(t,\cdot)}_{L^2(M;E)}^2=2\Real\int_0^t\pair{F(s,\cdot),\d_tF(s,\cdot)}_{L^2(M;E)}ds\leq\int_0^t\Big(\norm{F(s,\cdot)}_{L^2(M;E)}^2+\norm{\d_tF(s,\cdot)}_{L^2(M;E)}^2\Big)ds.\] Using the above in (\ref{edit_2.8}), together with Cauchy's inequality and the estimate (\ref{s2.3}), we deduce that\begin{equation}\label{s2.7}\norm{u(t)}_{H^1(M;E)}^2\leq C\Big(\int_0^t\norm{u(s)}_{H^1(M;E)}^2ds+\norm{F}_{H^{1,0}((0,T)\times M;E)}^2\Big).\end{equation}Then an application of Gr\"onwall's inequality tells us that\begin{equation}\label{s2.8}\norm{u}_{L^\infty(0,T;H^1(M;E))}\leq C\norm{F}_{H^{1,0}((0,T)\times M;E)}.\end{equation}

For the next estimate, we begin by applying $\d_t$ to (\ref{s2.1}). Using the expression (\ref{CLL}) for the connection Laplacian, we deduce that\begin{equation}\label{s2.9}(i\d_t+\Delta_A+V)\d_tu=\d_tF-2(\d_tA,du)_g+(\d_td^\ast A)u-(\d_tA,Au)_g-(A,(\d_tA)u)_g-(\d_tV)u.\end{equation}We now apply the estimate (\ref{s2.3}) to $\d_tu$, replacing $F$ appearing in (\ref{s2.3}) by the right-hand side of (\ref{s2.9}). We deduce, therefore, that\[\begin{split}&\norm{\d_tu}_{L^\infty(0,T;L^2(M;E))}\\\leq& C\norm{\d_tF-2(\d_tA,du)_g+(\d_td^\ast A)u-(\d_tA,Au)_g-(A,(\d_tA)u)_g-(\d_tV)u}_{L^2((0,T)\times M;E)}.\end{split}\]Using the estimates (\ref{s2.3}) and (\ref{s2.8}), on the right-hand side of the above, we observe that\begin{equation}\label{s2.10}\norm{\d_tu}_{L^\infty(0,T;L^2(M;E))}\leq C\norm{F}_{H^{1,0}((0,T)\times M;E).}\end{equation}

Lastly, we rearrange (\ref{s2.1}) to obtain \begin{equation}\begin{split}(\Delta_A+V)u=&\ F-i\d_tu\textrm{ in }(0,T)\times M\\u=&\ 0\textrm{ on }(0,T)\times\d M.\end{split}\end{equation}Then the bounds (\ref{s2.3}), (\ref{s2.8}) and (\ref{s2.10}) immediately imply the desired energy estimate\begin{equation}\label{NRGest}\norm{u}_{H^{1,2}((0,T)\times M;E)}\leq C\norm{F}_{H^{1,0}((0,T)\times M;E).}\end{equation}\end{proof}

We now turn to the proof of Proposition \ref{edit_p1}.

\begin{proof}[Proof of Proposition \ref{edit_p1}]

The unique solvability of the IBVP (\ref{1.1}) can be established in a similar manner to the scalar valued case, using the energy estimates of Proposition \ref{edit_p2}. Existence and uniqueness can then be proven using, for example, the Galerkin approach (see e.g. \cite[Theorem 2.3]{22} or \cite[Section 3, Theorem 10.1]{23}).\\

We now turn to establishing the bound (\ref{edit_DN}) for the Dirichlet-to-Neumann map. Recall the initial and boundary value problem (\ref{1.1}):\[\begin{split}(i\d_t+\Delta_{A}+V)u&=0\textrm{ in }(0,T)\times M,\\u(t,x)&=f\textrm{ on }(0,T)\times\d M,\\u(0,x)&=0\textrm{ in }M,\end{split}\]where the inhomogeneous Dirichlet data is given by $f\in C^\infty((0,T)\times\d M;E)$ satisfying $f\vert_{t=0}=\d_tf\vert_{t=0}=0$.\\

Note that we can find $\Phi\in C^\infty((0,T)\times M;E)$ such that\[\Phi(0,\cdot)=\d_t\Phi(0,\cdot)=0\textrm{ in M},\qquad\Phi=f\textrm{ on }\d M,\]and\begin{equation}\norm{\Phi}_{H^{3,2}((0,T)\times M;E)}\leq C\norm{f}_{H^{\frac{9}{4},\frac{3}{2}}((0,T)\times \d M;E)},\end{equation} for some $C>0$, depending only on $M$ and $T$. See \cite[Chapter 4, Section 2]{23} for a proof of this fact in the scalar case. The proof for vectors is analogous and, therefore, omitted. From the above, it holds that\begin{equation}\label{s2.14}F:=-(i\d_t+\Delta_A+V)\Phi\end{equation}satisfies $F(0,\cdot)=0$ in $M$. Then, letting $v$ be the solution of (\ref{s2.1}) corresponding to the source term $F$ defined in (\ref{s2.14}), we see that $u=\Phi+v$ is a solution to (\ref{1.1}). Then, it follows by (\ref{NRGest}) that\[\norm{u}_{H^{1,2}((0,T)\times M;E)}\leq C\norm{f}_{H^{\frac{9}{4},\frac{3}{2}}((0,T)\times \d M;E)}\]and applying this estimate with $f=0$ implies that such a solution $u$ is unique. Finally, we observe that\[\norm{\Lambda_{A,V}f}_{L^2((0,T)\times\d M;E)}\leq\norm{u}_{H^{1,2}((0,T)\times M;E)}\leq C\norm{\Phi}_{H^{3,2}((0,T)\times M;E)}\leq C\norm{f}_{H^{\frac{9}{4},\frac{3}{2}}((0,T)\times \d M;E)}.\]\end{proof}

\section{Construction of Gaussian Beam Solutions}
In this section, we shall construct Gaussian beam solutions to the Schr\"odinger equation which concentrate along geodesics in the high frequency limit.\\

Let $(\widehat{M},g)$ be a closed manifold. Recall that for a geodesic segment $\gamma:(a,b)\rightarrow\widehat{M}$ with no closed loops, there exist only finitely many values of $r\in(a,b)$ for which $\gamma$ self-intersects at $\gamma(r)$. We begin by recording the following system of Fermi coordinates near a geodesic, which we shall later use to construct our Gaussian beam solutions.\\

\begin{nobreak}\begin{lemma}\label{Fermi}
Let $(\widehat{M},g)$ be a compact $m$-dimensional manifold without boundary, $m\geq2$, and assume that $\gamma:(a,b)\rightarrow\widehat{M}$ is a unit-speed geodesic with no closed loops. Given a closed sub-interval $[a_0,b_0]$ of $(a,b)$ such that $\gamma\vert_{[a_0,b_0]}$ self-intersects only at $\gamma(r_j)$ with $a_0<r_1<\cdots<r_K<b_0$, and setting $r_0=a_0$, $r_{N+1}=b_0$, there exists an open cover $\{U_j,\phi_j\}_{j=0}^{K+1}$ of $\gamma([a_0,b_0])$ consisting of coordinate neighbourhoods with the following properties:
\begin{itemize}
\item $\phi_j(U_j)=I_j\times B$, where $I_j$ are open intervals and $B=B(0,\delta')$ is an open ball in $\R^{m-1}$, where $\delta'$ can be taken arbitrarily small.
\item $\phi_j(\gamma(r))=(r,0)$ for $r\in I_j$
\item $r_j$ only belongs to $I_j$ and $\overline{I_j}\cap\overline{I_k}=\emptyset,$ unless $\abs{j-k}\leq1$
\item $\phi_j=\phi_k$ on $\phi_j^{-1}((I_j\cap I_k)\times B)$
\end{itemize}
Furthermore, the metric in these coordinates satisfies $g^{jk}\vert_{\gamma(r)}=\delta^{jk}$ and $\d_ig^{jk}\vert_{\gamma(r)}=0$.
\end{lemma}\end{nobreak}

\begin{proof}
See e.g. \cite[Lemma 3.5]{C4} for details.
\end{proof}

We now turn to the construction of the Gaussian beam solutions. We consider here a non-tangential unit-speed geodesic in $M$ given by $\gamma:[0,L]\rightarrow M$. That is, $\gamma'(0)$ and $\gamma'(L)$ are both non-tangential to $\d M$, and $\gamma(r)\in M^{int}$ for $0<r<L$. Note that this implies that the geodesic $\gamma$ is not a closed loop in $M$.\\

We may then embed $(M,g)$ in some closed manifold $\widehat{M}$, and extend $\gamma$ to $\widehat{M}$ as a unit-speed geodesic $\gamma:[-\varepsilon,L+\varepsilon]\rightarrow\widehat{M}$. Our aim is to construct a Gaussian beam solution near $\gamma([0,L])$. We fix a point $x_0$ on $\gamma$, and apply Lemma \ref{Fermi} on $\widehat{M}$ with $a_0<0$ and $b_0>L$ chosen so that $\gamma(a_0)$ and $\gamma(b_0)$ are in the interior of $\widehat{M}\setminus M$. This gives us a system of coordinates $(r,y)$ around $x_0=(r_0,0)$, defined in a set $U=\{(r,y):\abs{r-r_0}<\delta,\ \abs{y}<\delta'\}$ such that the geodesic near $x_0$ is given by $\Gamma=\{(r,0):\abs{r-r_0}<\delta\}$.\\

The main aim of the present section is to establish the following result.

\begin{proposition}\label{edit_p3}
Let $(M,g)$ be non-trapping, let $\gamma:[0,L]\rightarrow M$ be a non-tangential geodesic, and let $s\gg 1$. There exists a function $v\in C^\infty((0,T)\times M;E)$, supported in a tubular neighbourhood of $\gamma$, which is an approximate solution of the Schr\"odinger equation in the sense that\begin{equation}\label{edit_p3est}\norm{(i\d_t+\Delta_{A}+V)v}_{L^2((0,T)\times M;E)}=O(s^{-1}),\qquad\norm{v}_{L^2((0,T)\times M;E)}=O(1).\end{equation}Away from the self-intersections of $\gamma$, this approximate solution is given by $v(t,r,y)=e^{is(\Psi(r,y)-st)}a(t,r,y)$, where the phase function $\Psi$ has the form $\Psi(r,y)=r+\frac{1}{2}H(r)y\cdot y+O(\abs{y}^3)$ and the amplitude satisfies\[a(t,r,0)=c_0\tilde{\chi}(t)e^{-\frac{1}{2}\int_{r_0}^r\tr H(\tilde{r})d\tilde{r}}\ U_{A}w+O(s^{-1})\]for some choices of arbitrary constant $c_0$, compactly supported smooth function $\tilde{\chi}$ and initial vector $w$.
\end{proposition}

\begin{proof}

The first step is the local construction of such an approximate solution. That is, we wish to construct a solution $v$ of the Schr\"odinger equation in $U$, with the form\[v=e^{is(\Psi(r,y)-st)}a(s;t,r,y),\]where $\Psi\in C^\infty(M;\C)$, $a\in C^\infty((0,T)\times M;E)$ are given near $\Gamma$, with $a$ supported in $\{\abs{y}<\delta'/2\}$. For convenience, we shall supress the dependence on $s$ of the amplitude function $a$.\\

In practice, we will determine the functions $\Psi,a$ by solving certain Eikonal and transport equations up to $N$th order on $\Gamma$, much as one would in the analogous construction for the stationary Schr\"odinger equation (see e.g. \cite[Theorem 5.4]{MC1}). As a result of using the phase $e^{is(\Psi-st)}$, we obtain the same Eikonal equation for $\Psi$ as in \cite{MC1}, but derive easier transport equations for the amplitude $a$ (compared to \cite[Theorem 5.4]{MC1}, where $a$ satisfies transport equations of $\overline{\d}$-type).\\

Therefore, let us begin by deriving these Eikonal and transport equations. Recalling the expression (\ref{CLL}) for the connection Laplacian, we first compute the Schr\"odinger operator applied to $v$:\begin{equation}\begin{split}\label{WKB}(i\d_t+\Delta_A+V)v=&e^{is(\Psi-st)}(i\d_t+\Delta_A+V)a+s^2e^{is(\Psi-st)}\Big(1-(d\Psi,d\Psi)_g\Big)a\\&+2ise^{is(\Psi-st)}\Big((d\Psi,\nabla^Aa)_g+\frac{1}{2}(\Delta_g\Psi)a\Big).\end{split}\end{equation}

In light of the above, we seek $\Psi$ satisfying the Eikonal equation\begin{equation}\label{Eikonal}(d\Psi,d\Psi)_g-1=0\textrm{ to $N$th order on }\Gamma,\end{equation}exactly as in the case of the stationary Schr\"odinger equation. On the other hand, the amplitude $a$ should satisfy, up to a small error, the transport equation\[s\Big((d\Psi,\nabla^Aa)_g+\frac{1}{2}(\Delta_g\Psi)a\Big)-\frac{i}{2}(i\d_t+\Delta_A+V)a=0\quad\textrm{to $N$th order on }\Gamma.\]\\

Thus, we begin by seeking a solution $\Psi$ of (\ref{Eikonal}) having the form $\Psi=\sum_{j=0}^{N}\Psi_j$, where\[\Psi_j(r,y)=\sum_{\abs{\alpha}=j}\frac{\Psi_{j,\alpha}(r)}{\alpha!}y^\alpha.\]Let us also write the metric in the form $g^{jk}=\sum_{l=0}^Ng^{jk}_l+r^{jk}_{N+1}$, where\[g^{jk}_l(r,y)=\sum_{\abs{\beta}=l}\frac{g^{jk}_{l,\beta}(r)}{\beta!}y^\beta,\quad r^{jk}_{N+1}=O(\abs{y}^{N+1}).\]

By the properties of the Fermi coordinates, we observe that $g^{jk}_0=\delta^{jk}$ and $g^{jk}_1=0$. Thus, we can immediately choose $\Psi_0(r)=r$ and $\Psi_1(r,y)=0$. Then, for $j,k=1...m$ and $\alpha,\beta=2...m$, we have\begin{equation}\begin{split}\label{G2}g^{jk}\d_j\Psi\d_k\Psi-1&=(1+g^{11}_2+\cdots)(1+\d_r\Psi_2+\cdots)(1+\d_r\Psi_2+\cdots)\\&+2(g^{1\alpha}_2+\cdots)(1+\d_r\Psi_2+\cdots)(\d_{y^\alpha}\Psi_2+\cdots)\\&+(\delta^{\alpha\beta}+g^{\alpha\beta}_2+\cdots)(\d_{y^\alpha}\Psi_2+\d_{y^\alpha}\Psi_3+\cdots)(\d_{y^\beta}\Psi_2+\d_{y^\beta}\Psi_3+\cdots)-1\\&=[2\d_r\Psi_2+\nabla_y\Psi_2\cdot\nabla_y\Psi_2+g^{11}_2]\\&+\sum_{p=3}^N\Big[2\d_r\Psi_p+2\nabla_y\Psi_2\cdot\nabla_y\Psi_p+\sum_{l=0}^pg^{11}_l \sum_{\substack{{j+k=p-l}\\{0\leq j,k<p}}}\d_r\Psi_j\d_r\Psi_k\\&+2\sum_{l=2}^pg^{1\alpha}_l\sum_{\substack{{j+k=p+1-l}\\{2\leq k<p}\\{0\leq j<p}}}\d_r\Psi_j\d_{y^\alpha}\Psi_k+\sum_{l=0}^{p-2}g^{\alpha\beta}_l\sum_{\substack{{j+k=p+2-l}\\{2\leq j,k<p}}}\d_{y^\alpha}\Psi_j\d_{y^\beta}\Psi_k\Big]+O\big(\abs{y}^{N+1}\big).\end{split}\end{equation}
In the last equality, we have chosen to collect the terms into homogeneous polynomials in $y$ (so that the first term is the second degree part of the right-hand side, and the rest are the parts of degree $p=3,\dots,N$). We first choose $\Psi_2$ such that the second-degree term $[2\d_r\Psi_2+\nabla_y\Psi_2\cdot\nabla_y\Psi_2+g^{11}_2]$ vanishes.\\

To this end, we choose $\Psi_2(r,y)=\frac{1}{2}H(r)y\cdot y$, where $H$ is a smooth, symmetric, complex matrix solving the matrix Riccati equation\begin{equation}\label{Riccati}H'(r)+H^2(r)=F(r),\end{equation}and $F(r)$ is the symmetric matrix such that $g^{11}_2(r,y)=-F(r)y\cdot y$. If we impose some initial condition $H(r_0)=H_0$ on this equation, where $H_0$ is  chosen to be a complex symmetric matrix with $\Imag(H_0)$ positive definite, then \cite[Lemma 2.56]{C5} implies that the matrix Riccati equation above has a unique smooth symmetric solution $H(r)$, for which $\Imag(H(r))$ is positive definite.\\

We now choose $\Psi_3$ so that the term corresponding to $p=3$ in the right-hand side of (\ref{G2}) vanishes. We obtain the equation\[2\d_r\Psi_3+2\nabla_y\Psi_2\cdot\nabla_y\Psi_3=F(r,y),\]where $F$ is a third-order polynomial in $y$ which only depends on $\Psi_2$ and $g$. This gives us a linear system of first-order ODEs for the Taylor coefficients $\Psi_{3,\alpha}(r)$, which can be solved uniquely if we prescribe some initial conditions at $r_0$. We may, then, repeat this argument in order to obtain $\Psi_4,\dots,\Psi_N$ by solving ODEs on $\Gamma$, given inital conditions at $r_0$.\\

Thus, we have $\Psi(r,y)=r+\frac{1}{2}H(r)y\cdot y+\tilde{\Psi}$, where $\tilde{\Psi}=O(\abs{y}^3)$. We now turn to finding the amplitude $a$ such that, up to a small error, we have\[s\Big((d\Psi,\nabla^Aa)_g+\frac{1}{2}(\Delta_g\Psi)a\Big)-\frac{i}{2}(i\d_t+\Delta_A+V)a=0\quad\textrm{to $N$th order on }\Gamma.\]

We choose $a$ of the form\[a=s^{\frac{m-1}{4}}(a_0+s^{-1}a_{1}+\cdots+s^{-N}a_{N})\chi(y/\delta'),\]where $\chi$ is a smooth function such that $\chi=1$ for $\abs{y}\leq\nicefrac{1}{4}$ and $\chi=0$ for $\abs{y}\geq\nicefrac{1}{2}$. Letting $\eta=\Delta_g\Psi$, it is enough to find $a_j$ such that\begin{equation}\begin{split}\label{transport}(d\Psi,\nabla^Aa_0)_g+\frac{1}{2}\eta a_0&=0\quad\textrm{to $N$th order on $\Gamma$}\\(d\Psi,\nabla^Aa_1)_g+\frac{1}{2}\eta a_1-\frac{i}{2}(i\d_t+\Delta_A+V)a_0&=0\quad\textrm{to $N$th order on $\Gamma$}\\\vdots\\(d\Psi,\nabla^Aa_N)_g+\frac{1}{2}\eta a_N-\frac{i}{2}(i\d_t+\Delta_A+V)a_{N-1}&=0\quad\textrm{to $N$th order on $\Gamma$}.\end{split}\end{equation}

We can write $\eta=\sum_{l=0}^N\eta_l+r_{N+1}$ and $a_0=a_{00}+\cdots+a_{0N}$, where each $\eta_l$, $a_{0l}$ is a homogenenous polynomial of order $j$ in $y$, and the remainder $r_{N+1}$ is $O(\abs{y}^{N+1})$. Thus, writing $A=A(\gamma')dr+A(\d_{y^\alpha})dy^\alpha$, we can rewrite the transport equation for $a_0$ in the system (\ref{transport}) above as\begin{equation}\begin{split}\label{edit_transport}\Big(1+g^{11}_2+\cdots\Big)\Big(1+\d_r\Psi_2+\cdots\Big)\Big(\d_r+A(\gamma')\Big)\Big(a_{00}+a_{01}+\cdots\Big)\\+\Big(g^{1\alpha}_2+\cdots\Big)\Big(1+\d_r\Psi_2+\cdots\Big)\Big(\d_{y^\alpha}+A(\d_{y^\alpha})\Big)\Big(a_{00}+a_{01}+\cdots\Big)\\+\Big(g^{\alpha1}_2+\cdots\Big)\Big(1+\d_{y^\alpha}\Psi_2+\cdots\Big)\Big(\d_r+A(\gamma')\Big)\Big(a_{00}+a_{01}+\cdots\Big)\\+\Big(\delta^{\alpha\beta}+g^{\alpha\beta}_2+\cdots\Big)\Big(\d_{y^\alpha}\Psi_2+\d_{y^\alpha}\Psi_3+\cdots\Big)\Big(\d_{y^\beta}+A(\d_{y^\beta})\Big)\Big(a_{00}+a_{01}+\cdots\Big)\\+\frac{1}{2}(\eta_0+\eta_1+\cdots)(a_{00}+a_{01}+\cdots)\\=0.\end{split}\end{equation}

We can then write $A(\theta)=\sum_{l=0}^NA_l(\theta)+R_{N+1}(\theta)$, where the entries of each $A_l(\theta)$ are homogeneous polynomials in $y$ of order $l$, and the remainder $R_{N+1}(\theta)$ is $O(\abs{y}^{N+1})$. As a result, we deduce that the left-hand side of (\ref{edit_transport}) becomes

\[\begin{split}\d_ra_{00}+A_0(\gamma')a_{00}+\frac{1}{2}\eta_0a_{00}\\+\d_ra_{01}+A_0(\gamma')a_{01}+A_1(\gamma')a_{00}+\nabla_y\Psi_2\cdot\nabla_ya_{01}+\nabla_y\cdot A_0(\gamma')a_{00}+\frac{1}{2}\eta_0a_{01}+\frac{1}{2}\eta_1a_{00}\\+\cdots\end{split}\]

We wish to find $a_{00}$ such that the first line of the above expression vanishes. To this end, we note that $\eta_0(r)=\Delta_g\Psi(r,0)=g^{\alpha\beta}\d_{y^\alpha}[H_{\alpha\beta}y^\alpha]=\tr H(r)$. We therefore choose $a_{00}$ such that\[\d_ra_{00}+A_0\big(\gamma'(r)\big)a_{00}+\frac{1}{2}\big(\tr H(r)\big)a_{00}=0.\]Since $A_0(t,r)=A(t,r,0)$, this equation has the solution\[a_{00}(t,r)=c_0\tilde{\chi}(t)e^{-\frac{1}{2}\int_{r_0}^r\tr H(\tilde{r})d\tilde{r}}\cdot U_Aw,\]where $\tilde{\chi}\in C_0^\infty((\tau,T-\tau))$ satisfies $\tilde{\chi}=1$ on $[2\tau, T-2\tau]$, $0\leq\tilde{\chi}\leq1$ and $\norm{\tilde{\chi}}_{W^{k,\infty}(\R)}\leq C_k\tau^{-k}$ with $C_k$ independent of $\tau$, and where $w$ is some arbitrary initial vector and $c_0$ is some initial constant which determines the value of $a_{00}(t,r_0)$. Since it will simplify later calculations, we choose the value\begin{equation}\label{c_0}c_0=\frac{\sqrt[4]{\det\Imag(H(r_0))}}{\sqrt{\Big(\int_{\R^{m-1}}e^{-\abs{y}^2}dy\Big)}}=\frac{\sqrt[4]{\det\Imag(H(r_0))}}{\pi^{\frac{m-1}{4}}}.\end{equation}We note that the constructed $a_{00}$ is smooth in time, since the connection form $A$ is smooth in time and the transport equation for $a_{00}$ is well-posed. We can then obtain the rest of $a_{01},\dots,a_{0N}$ by solving linear first-order ODEs. The sections $a_1,\dots,a_N$ may then be determined in much the same manner as $a_0$, so that the transport equations in (\ref{transport}) are satisfied to $N$th degree on $\Gamma$, and the smooth time-dependence of the amplitude $a$ follows from the smooth time-dependence of the connection form and the well-posedness of the transport equations that determine $a$.\\

Thus, we have constructed a function $v=e^{is(\Psi-st)}a$ in $U$ such that:\begin{align*}&\Psi(r,y)=r+\frac{1}{2}H(r)y\cdot y+\tilde{\Psi},\ \textrm{where }\tilde{\Psi}=O(\abs{y}^3),\\&a(t,r,y)=s^{\frac{m-1}{4}}(a_0+s^{-1}a_1+\cdots+s^{-N}a_N)\chi(y/\delta'),\\&a_0(t,r,0)=c_0\tilde{\chi}(t)e^{-\frac{1}{2}\int_{r_0}^r\tr H(\tilde{r})d\tilde{r}}\ U_{A}w.\end{align*}We now turn to establishing the bounds (\ref{edit_p3est}) for the function $v$ in $U$. Henceforth, we shall choose $N=5$. To begin with, note that\[\abs{e^{is(\Psi-st)}}=e^{-s\Imag\Psi}=e^{-\frac{1}{2}s\Imag(H(r))y\cdot y}e^{-sO(\abs{y}^3)}.\]

Note also that $\Imag(H(r))y\cdot y\geq c\abs{y}^2$ for $(r,y)\in U$, where the constant $c>0$ depends on $H_0$ and the value of $\delta$ appearing in the definition of $U$. Thus, it follows for $(r,y)\in U$ that, provided $y$ is sufficiently small, we have\[\abs{e^{is(\Psi(r,y)-st)}}\leq Ce^{-\frac{1}{4}cs\abs{y}^2}.\]From this fact, together with the definition of $\chi(y/\delta')$ and $\tilde{\chi}(t)$, we can deduce that for $r$ in a compact interval, possibly after decreasing $\delta'$, we have\[\abs{v(t,r,y)}\leq Cs^{\frac{m-1}{4}}e^{-\frac{1}{4}cs\abs{y}^2}\chi(y/\delta'),\quad s\gg1.\]

As a result, we have that for $s\gg1$,\begin{equation}\label{bnd1}\norm{v}_{L^2((0,T)\times U;E)}\leq C\norm{s^{\frac{m-1}{4}}e^{-\frac{1}{4}cs\abs{y}^2}}_{L^2((0,T)\times U;E)}=O(1).\end{equation}

Further, since $\Psi,a$ satisfy (\ref{Eikonal}) and (\ref{transport}) to $5$th order on $\Gamma$, we deduce from (\ref{WKB}) that\[\abs{(i\d_t+\Delta_A+V)v}\leq Cs^{\frac{m-1}{4}}e^{-\frac{1}{4}cs\abs{y}^2}\chi(y/\delta')\Big(s^2\abs{y}^6+s^{-5}\Big).\]It follows, then, that\begin{equation}\label{bnd3}\norm{(i\d_t+\Delta_A+V)v}_{L^2((0,T)\times U;E)}\lesssim\norm{s^\frac{m-1}{4}e^{-\frac{1}{4}cs\abs{y}^2}\big(s^2\abs{y}^6+s^{-5}\big)}_{L^2((0,T)\times U;E)}=O\big(s^{-1}\big),\end{equation}and by the same method we may further deduce that\begin{equation}\label{bnd4}\norm{(i\d_t+\Delta_A+V)v}_{H^{1,0}((0,T)\times U;E)}\lesssim\norm{s^\frac{m-1}{4}e^{-\frac{1}{4}cs\abs{y}^2}\big(s^2\abs{y}^6+s^{-3}\big)}_{L^2((0,T)\times U;E)}=O\big(s\big).\end{equation}

Thus, we have established the result of Proposition \ref{edit_p3} for $v$ in $U$. It remains to show that we can construct an approximate Gaussian beam solution $v$ in $M$ by gluing together the approximate solutions which we have constructed in $U$.\\

To this end, recall the definition of the sets $U_j$ from Lemma \ref{Fermi}. Fix $\delta'$ and choose an open cover $U_0,\dots,U_{K}$ of $\gamma([a_0,b_0])$, with each $U_j$ corresponding to an interval $I_j$, as in Lemma \ref{Fermi}. We first find a function $v^{(0)}=e^{is(\Psi^{(0)}-st)}a^{(0)}$ in $U_0$ following the method above, with some fixed initial conditions at $r_0$ for the ODEs that determine $\Psi^{(0)}$ and $a^{(0)}$. We continue by choosing some $\tilde{r}_1\in I_0\cap I_1$ so that $\gamma(\tilde{r}_1)\in U_0\cap U_1$, and construct $v^{(1)}=e^{is(\Psi^{(1)}-st)}a^{(1)}$ in $U_1$ again by the above method, choosing our initial conditions for $\Psi^{(1)}$, $a^{(1)}$ at $\tilde{r}_1$ to be, respectively, the values of $\Psi^{(0)}$ and $a^{(0)}$ at $\tilde{r}_1$. In this manner, we can proceed to determine $v^{(K)}$. We choose a partition of unity $\{\rho_j(r)\}$ for $[a_0,b_0]$ corresponding to the family of intervals $\{I_j\}$, and let $\tilde{\rho}_j(r,y)=\rho_j(r)$ in $U_j$. We can then define\[v=\sum_{j=0}^K\tilde{\rho}_jv^{(j)}.\]

Since the ODEs for the phases and amplitudes have the same initial value in $U_j$ as in $U_{j+1}$, we can deduce that $v^{(j)}=v^{(j+1)}$ in $U_j\cap U_{j+1}$. Therefore, we conclude that the $L^2$-scale bounds (\ref{bnd1}) for $v$ and (\ref{bnd3})-(\ref{bnd4}) for $(i\d_t+\Delta_A+V)v$ follow with $U=M$ from the corresponding bounds in $U_l$ for each $v^{(l)}$.
\end{proof}

Before we proceed, we note also the following partition, which shall be useful later. Suppose that $p_1,\dots,p_{K'}$ are the distinct points where the geodesic self-intersects, $0<r_1<\cdots<r_K<L$ are the times where the geodesic self-intersects, and $V_1,\dots,V_{K'}$ are balls centered at $p_1,\dots,p_{K'}$ respectively. Note that in each $V_j$, the function $v$ is given by\[v\vert_{V_j}=\sum_{\gamma(t_l)=p_j}v^{(l)}.\] By considering the individual sets $U_l\setminus(\cup_{j=1}^{K'}V_j)$, we can then choose a finite cover $\{W_1,\cdots,W_{K''}\}$ of the remaining points of $(\cup_{l=1}^KU_l)\setminus(\cup_{j=1}^{K'}V_j)$, where each $W_k$ is a subset of $U_{l_k}$ for some $l_k$. In particular, we can choose $W_k$ which remain away from $\gamma(r_{l_k})$, so that for small enough $\delta'$ it follows that\[v\vert_{W_k}=v^{(l_k)}.\] Together, the $V_j$ and $W_k$ cover form a cover\begin{equation}\label{cover}\textrm{supp}(v)\cap M\subset \Big(\cup_{j=1}^{K'} V_j\Big)\cup\Big(\cup_{k=1}^{K''}W_k\Big).\end{equation}

\section{Determination of the Scattering Data}

Let us begin by establishing the following useful Lemma, which will aid us in the proof of Theorem \ref{t1}.

\begin{lemma}\label{edit_gaugefix}
Given any connection form $A$ in $M$, there is a gauge-equivalent connection form $\tilde{A}$ with the additional property that $\tilde{A}(\nu)\vert_{\d M}=0$.
\end{lemma}

\begin{proof}
Let us use boundary normal coordinates in a neighbourhood of $\d M$ (see e.g. \cite[Section 2.1]{C5}) given by $(y,r)\in\d M\times[0,\varepsilon]$, with $\varepsilon>0$. Then for $u\in C^\infty(M)$ it holds that $\d_ru\vert_{\d M}=-\d_\nu u\vert_{\d M}$. Further, in these coordinates, we can decompose the connection form as $A=A(\nu)dr+A(\d_{y^\alpha})dy^\alpha$. We fix some $\mu\in C^\infty([0,\varepsilon];\R)$ such that $\mu(r)=1$ for $r$ near $0$ and $\mu(r)=0$ for $r$ near $\varepsilon$. We may then define a local gauge $G\in C^\infty((0,T)\times\d M\times[0,\varepsilon];\C^{n\times n})$ via \[G(r,\cdot)=e^{-r\mu(r)A(\nu)\, (r,\cdot)},\] where the dependence on $t$ and $y$ is left implicit. Note that, since $A$ is skew-Hermitian, it follows that $G$ is unitary, and since $G(r)=\Id$ for $r$ near $\varepsilon$, we can extend $G$ to an element of $C^\infty((0,T)\times M;U(n))$. Moreover, observe that\[G\vert_{\d M}=\Id,\quad\d_rG\vert_{\d M}=-A(\nu)\vert_{\d M}.\] Let $\tilde{A}$ denote the gauge transform of $A$ by $G$, that is\[\tilde{A}=G^{-1}AG+G^{-1}dG\]and observe that\[\tilde{A}(\nu)\vert_{\d M}=A(\nu)\vert_{\d M}-A(\nu)\vert_{\d M}=0.\]\end{proof}

Since the data $\Lambda_{A_j,V_j}$ and $C_{A_j}$ are gauge invariant, Lemma \ref{edit_gaugefix} implies that we may assume without any loss of generality that $A_1(\nu)\vert_{\d M}=A_2(\nu)\vert_{\d M}=0$ in the proofs of Theorems \ref{t1} and \ref{t2}. With this in mind, we proceed to the proof of Theorem \ref{t1}.\\

For $j=1,2$, we construct Gaussian beam solutions $u_j$ of the Schr\"odinger equations\begin{equation}\begin{split}\big(i\d_t+\Delta_{A_j}+V_j\big)u_j(t,x)&=0\textrm{ in }(0,T)\times M,\\u_1(0,\cdot)=u_2(T,\cdot)&=0\textrm{ in }M.\end{split}\end{equation}To this end, we fix some geodesic $\gamma_{x,\theta}$ in $\widehat{M}$ for $x\in\d M$, and choose a system of Fermi coordinates along $\gamma_{x,\theta}$ as in Lemma \ref{Fermi}. Using the work of the previous section, we construct approximate solutions $v_j$ in $M$ having the form $e^{is(\Psi-st)}a^{(A_j)}$ away from the self-intersections of $\gamma_{x,\theta}$. We can turn these $v_j$ into exact solutions $u_j=v_j+R_j$ by solving\[\begin{split}\big(i\d_t+\Delta_{A_j}+V_j\big)R_j&=-(i\d t+\Delta_{A_j}+V_j)v_j\textrm{ in }(0,T)\times M,\\R_j&=0\textrm{ on }(0,T)\times\d M,\\R_1(0,\cdot)&=R_2(T,\cdot)=0\textrm{ in }M.\end{split}\]

Note that for $s\gg1$, (\ref{bnd3}) and the energy estimate (\ref{s2.3}) yield\begin{equation}\label{d1est}\norm{R_j}_{L^2((0,T)\times M;E)}\leq C\norm{(i\d_t+\Delta_{A_j}+V_j)v_j}_{L^2((0,T)\times M;E)}=O\big(s^{-1}\big),\end{equation}whereas (\ref{bnd4}) and the energy estimate (\ref{NRGest}) yield\[\norm{R_j}_{H^{0,2}((0,T)\times M;E)}\leq C\norm{(i\d_t+\Delta_{A_j}+V_j)v_j}_{H^{1,0}((0,T)\times M;E)}=O(s).\]

Interpolating between (\ref{d1est}) and the above, we conclude that for $s\gg1$ we have\begin{equation}\label{Rjbnd}\norm{R_j}_{H^{0,1}((0,T)\times M;E)}+s\norm{R_j}_{L^2((0,T)\times M;E)}=O(1).\end{equation}

We then set $\phi_j=u_j$ on $(0,T)\times\d M$ and consider $\omega\in H^{1,2}((0,T)\times M;E)$ the solution of
\begin{equation}\begin{split}\label{3.1}(i\d_t+\Delta_{A_2(t)}+V_2(t,x))\omega(t,x)&=0\textrm{ in }(0,T)\times M,\\\omega(t,x)&=\phi_1\textrm{ on }(0,T)\times\d M,\\\omega(0,\cdot)&=0\textrm{ in }M.\end{split}\end{equation}

We observe that the difference $\omega-u_1$ solves the following Schr\"odinger equation:

\begin{equation}\begin{split}\label{3.2}\big(i\d_t+\Delta_{A_2}+V_2\big)(\omega-u_1)&=2(A_1-A_2,du_1)_g+Qu_1\textrm{ in }(0,T)\times M,\\\omega-u_1&=0\textrm{ on }(0,T)\times\d M,\\\omega(0,x)-u_1(0,x)&=0\textrm{ in }M,\end{split}\end{equation}
where $Qu_1=(V_1-V_2)u_1+(A_1,A_1u_1)_g-(A_2,A_2u_1)_g-(d^\ast A_1)u_1+(d^\ast A_2)u_1$.\\

Taking the Hermitian inner product of the above equation with $u_2$, we deduce that \begin{equation}\label{3.3}\int_0^T\int_M\pair{2(A_1-A_2,du_1)_g+Qu_1,u_2}_EdV_gdt=\int_0^T\int_{\d M}\pair{\d_\nu(\omega-u_1),u_2}_Ed\sigma_gdt.\end{equation}

As a result of Lemma \ref{edit_gaugefix}, we may assume without loss of generality that $A_1(\nu)\vert_{\d M}=A_2(\nu)\vert_{\d M}=0$. Therefore, the right-hand side of the above is bounded by \begin{equation}\begin{split}\label{3.4}\abs{\int_0^T\int_{\d M}\pair{\d_\nu(\omega-u_1),u_2}_Ed\sigma_gdt}&\leq C\norm{(\Lambda_{A_1,V_1}-\Lambda_{A_2,V_2})\phi_1}_{L^2((0,T)\times\d M;E)}\norm{\phi_2}_{L^2((0,T)\times\d M;E)}\\&\leq C\norm{\Lambda_{A_1,V_1}-\Lambda_{A_2,V_2}}\norm{\phi_1}_{H^{\frac{9}{4},\frac{3}{2}}((0,T)\times\d M;E)}\norm{\phi_2}_{L^2((0,T)\times\d M;E)},\end{split}\end{equation}which vanishes when $\Lambda_{A_1,V_1}=\Lambda_{A_2,V_2}$. On the other hand, since $du_1=is(d\Psi)u_1+e^{is(\Psi-st)}da^{(A_1)}+dR_1$, the left-hand side of (\ref{3.3}) can be written as

\[\begin{split}\int_0^T\int_M\pair{2(A_1-A_2,du_1)_g+Qu_1,u_2}_EdV_gdt=is\int_{(0,T)\times M}\pair{2\big((A_1-A_2),(d\Psi)u_1\big)_g,u_2}_EdV_gdt\\+\int_{(0,T)\times M}\pair{2\big((A_1-A_2),e^{is(\Psi-st)}da^{(A_1)}+dR_1\big)_g+Qu_1,u_2}_EdV_gdt.\end{split}\]

We can divide the above by $s$ and use the bounds (\ref{Rjbnd}) and (\ref{bnd1}) to deduce that\begin{equation}\begin{split}\label{3.5}&\abs{\int_{(0,T)\times M}\pair{\big((A_1-A_2),(d\Psi)u_1\big)_g,u_2}_EdV_gdt}\\\leq &s^{-1}\abs{\int_{(0,T)\times M}\pair{\big((A_1-A_2),du_1\big)_g+Qu_1,u_2}_EdV_gdt}+O\big(s^{-1}\big).\end{split}\end{equation}

By combining (\ref{3.3}), (\ref{3.4}), and (\ref{3.5}), we conclude that when $\Lambda_{A_1,V_1}=\Lambda_{A_2,V_2}$ it follows that\[\abs{\int_{(0,T)\times M}\pair{\Big((A_1-A_2),(d\Psi)u_1\Big)_g\,,u_2}_EdV_gdt}\leq O(s^{-1}).\]Thus, we can let $s\rightarrow\infty$ in the right-hand side of the above to conclude that \begin{equation}\label{3.6}\lim_{s\rightarrow\infty}\int_{(0,T)\times M}\pair{\Big((A_1-A_2),(d\Psi)u_1\Big)_g\,,u_2}_EdV_gdt=0.\end{equation}

We now make use of the following Lemma:

\begin{lemma}\label{SPC}
\[\begin{split}&\lim_{s\rightarrow\infty}\int_{(0,T)\times M}\pair{\Big((A_1-A_2),(d\Psi)u_1\Big)_g\,,u_2}_EdV_gdt\\=\int_0^T&\tilde{\chi}^2\int_{0}^{\rho_+(x,\theta)}\pair{\Big(A_1\big(\gamma'_{x,\theta}(r)\big)-A_2\big(\gamma'_{x,\theta}(r)\big)\Big)U_{A_1}w_1,U_{A_2}w_2}_Edrdt.\end{split}\]
\end{lemma}

Before proving the above result, we first conclude the proof of Theorem \ref{t1}. Since $\tau\in(0,T/4)$ in the definition of $\tilde{\chi}$ is arbitrary, we deduce that\begin{equation}\label{3.8}\int_{0}^{\rho_+(y,\theta)}\pair{\Big(A_1\big(\gamma'_{x,\theta}(r)\big)-A_2\big(\gamma'_{x,\theta}(r)\big)\Big)U_{A_1}w_1,U_{A_2}w_2}_Edr=0.\end{equation}

Note that the scattering data $C_{A_j}$ takes values in $U(n)$, so that we can define $C_{A_j}^{-1}$ as the matrix inverse of $C_{A_j}$. Making use of (\ref{HIP}), (\ref{PT2}), and (\ref{ScatteringData}), it can be shown that\begin{equation}\begin{split}\label{3.9}\pair{\big(C_{A_2}^{-1}C_{A_1}-\Id\big)w_1,w_2}_E=&\int_{0}^{\rho_+(y,\theta)}\d_r\pair{U_{A_1}w_1,U_{A_2}w_2}_Edr\\=&\int_{0}^{\rho_+(y,\theta)}\pair{\nabla_{\gamma'_{x,\theta}}^{A_2}U_{A_1}w_1,U_{A_2}w_2}_Edr+\int_{0}^{\rho_+(y,\theta)}\pair{U_{A_1}w_1,\nabla_{\gamma'_{x,\theta}}^{A_2}U_{A_2}w_2}_Edr\\=&\int_{0}^{\rho_+(y,\theta)}\pair{\Big(A_2\big(\gamma'_{x,\theta}(r)\big)-A_1\big(\gamma'_{x,\theta}(r)\big)\Big)U_{A_1}w_1,U_{A_2}w_2}_Edr.\end{split}\end{equation}

However, (\ref{3.8}) and (\ref{3.9}) imply that $\pair{(C_{A_2}^{-1}C_{A_1}-\Id)w_1,w_2}_E=0$. Since $w_1,w_2$ were arbitrary, we conclude that $C_{A_2}^{-1}C_{A_1}=\Id$, whence $C_{A_1}=C_{A_2}$ and the proof of Theorem \ref{t1} is complete. It remains for us to provide a proof of Lemma \ref{SPC}.

\begin{proof}[Proof of Lemma \ref{SPC}]
By considering a partition of unity subordinate to the open cover (\ref{cover}), it is sufficient to consider $A_1-A_2$ compactly supported in $(V_j\cap M)$ or $(W_k\cap M)$. We shall first prove the latter case. Recall that $a^{(A_j)}=s^{\frac{m-1}{4}}(a_0+O(s^{-1}))\chi(y/\delta')$. We can observe that
\[\begin{split}\lim_{s\rightarrow\infty}&\int_{(0,T)\times M}\pair{\Big((A_1-A_2),(d\Psi)u_1\Big)_g\,,u_2}_EdV_gdt\\=\lim_{s\rightarrow\infty}&\int_0^T\int_{0}^{\rho_+(x,\theta)}\int_{\R^{m-1}}e^{-s\Imag(H(r))y\cdot y}e^{sO(\abs{y}^3)}s^{\frac{m-1}{2}}\chi^2(y/\delta')\times\\&\Big[\pair{\Big(A_1(t,r,y)-A_2(t,r,y),d\Psi(t,r,y) a_0^{(A_1)}(t,r,y)\Big)_g\,,a_0^{(A_2)}(t,r,y)}_E+O(s^{-1})\Big]\abs{g(r,y)}^\frac{1}{2}dydrdt.\end{split}\]

We make the substitution $y\mapsto s^{-\frac{1}{2}}y$ above, and recall that $d\Psi\vert_{y=0}=dr$ is dual via the Riemannian metric to $\gamma'_{x,\theta}(r)$. Then, since $\Imag(H)$ is positive definite and $\delta'$ is small, we note that the exponential term involving $\Imag (H)$ dominates the others. Thus, we can conclude that the right-hand side of the above is given by\[\int_0^T\int_{0}^{\rho_+(x,\theta)}\Big(\int_{\R^{m-1}}e^{-\Imag(H(r))y\cdot y}dy\Big)\pair{\big(A_1(\gamma'_{x,\theta})-A_2(\gamma'_{x,\theta})\big)a_0^{(A_1)}(t,r,0),a_0^{(A_2)}(t,r,0)}_E\abs{g(r,0)}^{\frac{1}{2}}drdt.\]

Evaluating the integral in $y$ and using that $\abs{g(r,0)}=1$, we rewrite the above integral in the form

\[\int_0^T\tilde{\chi}^2\Big(\int_{\R^{m-1}}e^{-\abs{y}^2}dy\Big)\int_{0}^{\rho_+(x,\theta)}\frac{\abs{c_0}^2e^{-\int_{r_0}^r\tr\Real(H(s))ds}}{\sqrt{\det\Imag(H(r))}}\pair{\Big(A_1(\gamma'_{x,\theta})-A_2(\gamma'_{x,\theta})\Big)U_{A_1}w_1,U_{A_2}w_2}_Edrdt.\]

We now use the fact that, according to \cite[Lemma 2.58]{C5}, matrices which solve the Riccati equation (\ref{Riccati}) have the property\[\det\Imag(H(r))=\det\Imag(H(r_0))e^{-2\int_{r_0}^r\tr\Real(H(s))ds}.\]

This fact, together with our choice of constant $c_0$ in (\ref{c_0}) is sufficient to conclude that\[\begin{split}\lim_{s\rightarrow\infty}&\int_{(0,T)\times M}\pair{\big((A_1-A_2),(d\Psi)u_1\big)_g\,,u_2}_EdV_gdt\\=\int_0^T\tilde{\chi}^2&\int_{0}^{\rho_+(y,\theta)}\pair{\Big(A_1\big(\gamma'_{x,\theta}(r)\big)-A_2\big(\gamma'_{x,\theta}(r)\big)\Big)U_{A_1}w_1,U_{A_2}w_2}_Edrdt,\end{split}\]when $A_1-A_2$ is compactly supported in $W_k\cap M$. For $A_1-A_2$ compactly supported in $V_j\cap M$, we can write $v=\sum_{\gamma(t_l)=p_j}v^{(l)}$. Thus, the limit has the form\[\begin{split}\lim_{s\rightarrow\infty}\int_{(0,T)\times M}\pair{\big((A_1-A_2),(d\Psi)u_1\big)_g\,,u_2}_EdV_gdt\\=\sum_{\gamma(t_l)=p_j}\lim_{s\rightarrow\infty}\int_{(0,T)\times M}\pair{\big((A_1-A_2),(d\Psi)v_1^{(l)}\big)_g\,,v_2^{(l)}}_EdV_gdt\\+\sum_{\substack{l\neq l'\\{\gamma(t_l)=\gamma(t_{l'})=p_j}}}\lim_{s\rightarrow\infty}\int_{(0,T)\times M}\pair{\big((A_1-A_2),(d\Psi)v_1^{(l)}\big)_g\,,v_2^{(l')}}_EdV_gdt.\end{split}\]

We observe that the first sum converges to the required limit by the same compuations used to prove the limit in  $W_k\cap M$, and the second sum vanishes via stationary phase arguments as in the proof of \cite[Proposition 3.1]{C4}, thus completing the proof of Lemma \ref{SPC}.\end{proof}

\section{Proof of Gauge Equivalence}

In what follows, we denote by $\varphi_r(x,v)$ the geodesic flow given by $\varphi_r(x,v)=(\gamma_{x,v}(r),\gamma'_{x,v}(r))\in TM$. Additionally, we denote by $X$ the geodesic vector field on $M$, which satisfies $\d_r(F\circ\varphi_r)=X(F)\circ\varphi_r$, for any function $F:SM\mapsto\C^n$.\\

For a function $F:SM\rightarrow\C^n$ and a connection $1$-form $B:TM\mapsto\C^{n\times n}$, we can consider the transport equation:
\begin{equation}\begin{split}\label{ray}Xw+Bw=&-F\textrm{ in }SM\\w\vert_{\d_-SM}=&0,\end{split}\end{equation}where $B$ acts on $w:SM\rightarrow\C^n$ by multiplication at each point. In the present work, we shall consider only the cases $F(x,\theta)=f(x)$ or $F(x,\theta)=\alpha_j(x)\theta^j$, where $f,\alpha_j:M\rightarrow\C^n$.  Using the transport equation (\ref{ray}), we can define the attenuated ray transform of $F$ with attenuation due to $B$ (see e.g. \cite[Section 1]{C1}), via:\begin{equation}\label{ARTDefn}I_BF=w\vert_{\d_+SM}.\end{equation}

We can now proceed to the proof of Theorem \ref{t2} below.

\begin{proof}[Proof of Theorem \ref{t2}]
In order to prove Theorem \ref{t2}, we further assume that $M$ is either i) $2$-dimensional and simple, or ii) of dimension $m\geq3$ with strictly convex boundary and admits the existence of a smooth strictly convex function $\phi:M\rightarrow\R$. We note that this second condition is conjectured to be true for all simple manifolds (amongst others, see e.g. \cite[Section 2]{C1} for discussion), although the question is still open at present.\\

Consider the candidate gauge $G(t)=U_{A_1(t)}(U_{A_2(t)})^{-1}$. Note that $G(t)$ is unitary, $G(t):SM\rightarrow U(n)$, and depends smoothly on time. Further, since $C_{A_1}=C_{A_2}$, it holds that $G(t)\vert_{\d_+SM}=\Id$, and therefore that $G(t)\vert_{\d SM}=\Id$. We can observe that\begin{equation}\begin{split}\label{6.1}XG+A_1G-GA_2&=0\\G\vert_{\d SM}&=\Id.\end{split}\end{equation}

It remains to show that $G(t)$ depends only on the base-point $x\in M$, and further that $G(t)$ is smooth in $M$. Note that (\ref{6.1}) is equivalent to\begin{equation}\begin{split}\label{6.2}X(G-\Id)+A_1(G-\Id)-(G-\Id)A_2=A_2-A_1\\(G-\Id)\vert_{\d SM}=0.\end{split}\end{equation}

We can henceforth fix some $t\in(0,T)$ and define a new connection form $B$ via $B(W)=A_1W-WA_2$ for $W:SM\rightarrow\C^{n\times n}$. Then (\ref{6.2}) becomes\begin{equation}\begin{split}\label{ARTGauge}X(G-\Id)+B(G-\Id)=A_2-A_1\\(G-\Id)\vert_{\d SM}=0.\end{split}\end{equation}

We can interpret the above equation as a ray transform (see e.g. \cite[Section 1]{C1}), as $I_B(A_1-A_2)=0$, where $I_B$ denotes the ray transform with attenuation due to the connection form $B$. We emphasize here that $A_1$, $A_2$ and $B$ are smooth, whereas we have not yet shown that $G(t)$ in (\ref{ARTGauge}) is smooth.\\

We now apply either \cite[Theorem 1.3]{C3} if $M$ is $2$-dimensional and simple, or \cite[Theorem 1.6]{C1} if $M$ is of dimension $m\geq3$ with strictly convex boundary and admits a smooth strictly convex function. Thus, we conclude that $A_2(t)-A_1(t)=\nabla^{B(t)}p(t)$ for some smooth function $p(t):M\rightarrow \C^{n\times n}$, and (\ref{ARTGauge}) then implies that $G(t)=\Id+p(t)$. Hence, we have shown that $G(t)$ depends only on the base-point $x\in M$, and the smoothness of $G(t)$ in $M$ follows from the smoothness of $p(t)$. Therefore, $G$ satisfies all the necessary conditions to be a gauge, and it follows from (\ref{6.1}) that $A_1$ is gauge-equivalent to $A_2$ via the gauge-transform $A_2=G^{-1}dG+G^{-1}A_1G$, as required.\\

We now turn to showing unique determination of the potential. Note that since $A_2=G^{-1}dG+G^{-1}A_1G$, it remains only to show that $V_2=G^{-1}V_1G+iG^{-1}\d_tG$. We define $V_3=G^{-1}V_1G+iG^{-1}\d_tG$. By gauge invariance, it holds that $\Lambda_{A_1,V_1}=\Lambda_{A_2, V_3}$. Thus, by the assumption $\Lambda_{A_1,V_1}=\Lambda_{A_2,V_2}$ it follows that\begin{equation}\label{6.3}\Lambda_{A_2,V_3}=\Lambda_{A_2,V_2}.\end{equation}

It remains to show that condition (\ref{6.3}) implies that $V_2=V_3$. Note that we can take $A_1=A_2=A$ in (\ref{3.3}), and use the fact that $\Lambda_{A,V_3}=\Lambda_{A,V_2}$ to deduce that \[\int_{(0,T)\times M}\pair{(V_3-V_2)u_1,u_2}_E\circ\varphi_r(y,\theta)dV_gdt=0.\] We can let $s\rightarrow\infty$ in the left-hand side above, and applying the argument of Lemma \ref{SPC} we deduce that\[\int_0^T\tilde{\chi}^2(t)\int_0^{\rho_+(y,\theta)}\pair{(V_3-V_2)U_Aw_1,U_Aw_2}_E\circ\varphi_r(y,\theta)drdt=0.\] Since the choice of $\tilde{\chi}$ is arbitrary, we conclude that for $V=V_3-V_2$ and for each $t\in(0,T)$ and $(y,\theta)\in\d_+SM$ we have\[\int_{0}^{\rho_+(y,\theta)}\pair{VU_Aw_1,U_Aw_2}_E\circ\varphi_r(y,\theta)dr=0,\] where $\varphi_r(y,\theta)$ is the geodesic flow defined by $\varphi_r(y,\theta)=(\gamma_{y,\theta}(r),\gamma'_{y,\theta}(r))$. Then, by linearity, we conclude that\begin{equation}\label{6.4}\int_{0}^{\rho_+(y,\theta)}{U_A}^{-1}VU_A\circ\varphi_r(y,\theta)dr=0.\end{equation}

In order to finish the proof, we wish to interpret (\ref{6.4}) as an attenuated ray transform.\\

For $W\in C^\infty(SM;\C^{n\times n})$, we can define the map $BW=AW-WA=[A,W]$. Since $A(x,v)=A_j(x)v^j$, we can observe that $BW(x,v)=[A_j(x),W]v^j$.\\

Recall that the inner product $\pair{\cdot,\cdot}_E$ on the trivial bundle $E=M\times\C^n$ induces an inner product on the endomorphism bundle $\End(E)=M\times\C^{n\times n}$ via $\pair{X,Y}_{\End(E)}=\tr(X^\ast Y)$. Note that we can regard the $1$-form $B$ as a connection form on $\End(E)$. Further, by the cyclic property of trace, we observe that \[\pair{BX,Y}_{\End(E)}=\tr(-X^\ast AY+AX^\ast Y)=\tr(-X^\ast AY+X^\ast YA)=\pair{X,-BY}_{\End(E)},\] and that $B$ is, therefore, a unitary connection on the endomorphism bundle.\\

Letting $X$ once again denote the geodesic vector field, we consider $\omega$ the solution of the transport equation
\begin{equation}\label{6.5}X\omega+B\omega=-V\textrm{ in }SM,\quad \omega=0\textrm{ on }\d_-SM.\end{equation}Using (\ref{6.5}), we now observe that \[X({U_A}^{-1}\omega U_A)={U_A}^{-1}A\omega U_A+{U_A}^{-1}(-V-B\omega)U_A-{U_A}^{-1}\omega AU_A=-{U_A}^{-1}VU_A.\]

The above then implies that $\d_r[({U_A}^{-1}\omega U_A)\circ\varphi_r(y,\theta)]=-{U_A}^{-1}VU_A\circ\varphi_r(y,\theta)$, using the definitions of the geodesic vector field and geodesic flow. We can integrate this last expression to obtain\begin{equation}\label{6.6}-\int_{0}^{\rho_+(y,\theta)}{U_A}^{-1}VU_A\circ\varphi_r(y,\theta)dr=({U_A}^{-1}\omega U_A)\circ\varphi_{\rho_+(y,\theta)}(y,\theta)-({U_A}^{-1}\omega U_A)\circ\varphi_{0}(y,\theta).\end{equation}

Note that $\varphi_{\rho_+(y,\theta)}(y,\theta)\in\d_-SM$ and $\varphi_{0}(y,\theta)=(y,\theta)\in\d_+SM$. By recalling that that $U_A=\Id$ on $\d_+SM$ and $\omega=0$ on $\d_-SM$, we observe that the right-hand side of (\ref{6.6}) is just $-\omega\vert_{\d_+SM}$.\\

Therefore, (\ref{6.4}) and (\ref{6.5}) tell us that\[I_BV=\omega\vert_{\d_+SM}=\int_{0}^{\rho_+(y,\theta)}{U_A}^{-1}VU_A\circ\varphi_r(y,\theta)dr=0.\]

Hence, if $M$ is $2$-dimensional and simple, we can apply \cite[Theorem 1.3]{C3} to conclude that the above implies $V=0$. On the other hand, if $M$ is of dimension $m\geq3$ with strictly convex boundary and admits a smooth strictly convex function, we instead apply the result of \cite[Theorem 1.6]{C1} to conclude that $V=0$. Thus, it holds that $V_2=G^{-1}V_1G+iG^{-1}\d_tG$, and hence $(A_1,V_1)$ is gauge-equivalent to $(A_2,V_2)$, as required.
\end{proof}

\section*{Acknowledgments}
AT was supported by EPSRC DTP studentship EP/N509577/1.

\begin{small}