\section{Introduction}

Reconfigurable photonic networks or meshes of interferometers are capable of transforming optical modes for quantum computation, energy-efficient deep learning, mode unscrambling, sensing, and beamforming \cite{Bogaerts2020ProgrammableCircuits}. In such applications, a set of $N$ waveguide modes containing data or bits are fed into a mesh network of reconfigurable photonic interferometers and phase shifters whose output represents a matrix multiplication, where the matrix is represented by the transmission matrix of the optical device in the basis of waveguide modes. The ability to perform matrix multiplication in an ultra-fast and energy efficient manner is an attractive property of reconfigurable photonic circuits, which significantly broadens the application space for such devices. This approach also promises a new way of performing self-aligning optics with no moving parts which has applications in high-speed sensing and optical phased arrays \cite{Miller2013Self-aligningCoupler, Miller2020AnalyzingNetworks}. However, a key problem with this architecture is the presence of systematic errors that build up over various stages of the photonic circuit, requiring the definition of a theory of sensitivity.

In this paper, we discuss a common class of reconfigurable photonic networks that contain a network or ``mesh'' of Mach-Zehnder interferometers (MZIs), which can be used to progressively interfere pairs of waveguide modes. These modes may be generated in a number of ways; machine learning and photonic computing might use integrated modulators on a photonic chip operated ideally at sub-nanosecond time scales using electro-optic modulation \cite{Wang2018NanophotonicModulatorsb} and using a laser source tuned to a wavelength of interest (e.g. 1550 nm). Alternatively, these modes can arrive directly from free space as outputs of photonic sensors. This idea is still relatively unexplored but is a major focus of recent work \cite{Miller2020AnalyzingNetworks} showing that multimode fields can be analyzed directly by training self-configurable networks of interferometers arranged in a ``binary tree'' configuration. Such modes convert a continuous free space intensity pattern into a discrete sequence of complex numbers with amplitude and phase characterizing light in the fundamental mode of a single-mode waveguide. Pairs of modes can be constructively and/or destructively interfered sequentially as they propagate through feedforward meshes of interferometers. This can be achieved programmatically by adjusting active phase shifting elements within MZIs placed in the photonic circuit. Such phase shifting elements, the pathway for programming analog matrix multiplication and signal processing in photonics, can operate by modulating the effective mode index thermally \cite{Harris2014EfficientSilicon}, electromechanically \cite{Errando-Herranz2020MEMSCircuits, Edinger2020CompactPhotonics}, using phase change co-integration \cite{Wuttig2017Phase-changeApplications}, or electrooptically \cite{Wang2018NanophotonicModulators}.

As shown in Fig. \ref{fig:node}(a), a triangular network of MZIs can be programmed such that the process of light propagating forward through the network is a physical (analog) implementation of any desired unitary matrix multiplication \cite{Reck1994ExperimentalOperator, Miller2013Self-configuringInvited} \footnote{The triangular unitary model was actually first discovered in 1897 by Hurwitz \cite{Hurwitz1897UberIntegration} and only popularized for physics applications over a century later.}. Here, we propose an alternate method for implementing matrix multiplication in classical photonic circuits. Our approach relies on ``nonlocal'' interactions that can be recursively defined using binary trees, a data structure commonly deployed in computer science applications. Depictions of ``unbalanced'' and``balanced'' binary trees of MZIs are shown in Fig. \ref{fig:node}(d). The recursive definition of such structures is further elaborated in Fig. \ref{fig:tree}, demonstrating the construction of arbitrary binary tree designs by connecting every MZI to an output or up to two binary subtrees.

At a high level, this binary tree ``vector unit'' construction can be programmed to route any incoming set of waveguide modes into a single waveguide \cite{Miller2013Self-aligningCoupler} (Fig. \ref{fig:node}(e, f, g)). The reverse is also true by reciprocity as light can be sent into a single mode waveguide and routed out such that any complex vector output may be generated in the device basis (Fig. \ref{fig:node}(a, b, c)). We specifically show that a vector unit can be represented by a recursive binary tree definition (Fig. \ref{fig:node}(d, h)), which can ultimately be used to explore a deeper mathematical framework for the error tolerance and dispersion of such networks. Such units have important roles for state preparation and readout of optical modes which can be invaluable for machine learning in hybrid neural networks \cite{Pai2022ExperimentallyNetworks} and cryptocurrency hash functions \cite{Pai2022ExperimentalCryptocurrency} as well as telecommunications and optical phased arrays.

Photonic networks that are universal can implement any unitary matrix, contain $N^2$ degrees of freedom, and typically include triangular and rectangular meshes \cite{Clements2016AnInterferometers}, the latter of which is more compact but not self-configurable. However, these networks typically only interact waveguide modes locally (only modes of neighboring waveguides interact), which can limit overall photonic bandwidth for which a certain matrix accuracy is achievable. This motivates exploring whether nonlocal architectures increase the bandwidth in a photonic mesh, which is explored in some detail in Ref. \citenum{Hamerly2021InfinitelyInterferometers}. While crossings for nonlocal connections may be a concern, recent advances in multilayer photonics and silicon nitride-on-silicon in commercial CMOS foundries make such designs plausible \cite{Chiles2017Multi-planarLoss}. Balanced ``binary tree''-based architectures are nonlocal architectures that are a key building block for low-depth photonic circuits such as the butterfly architecture, the Benes network (an arbitrary switching network that is two butterfly architectures back-to-back). They also have mappings to quantum architectures such as state preparation circuits \cite{Araujo2021APreparation}, quantum fast Fourier transform (FFT) \cite{Flamini2017BenchmarkingProcessing}, and the cosine-sine architecture, which is a useful architecture in quantum computation \cite{Mottonen2004QuantumGates}.

Our main contribution is to provide a new error model to explain why increased bandwidth and robustness arises from nonlocal connectivity in balanced binary trees. We generally propose ``wide'' or ``splay'' photonic architectures that provide more tolerance than ``deep'' universal architectures because light has to propagate through fewer devices on average, motivating a core design principle for photonic circuits. To prove these claims, we organize the paper as follows:
\begin{enumerate}
    \item In Section \ref{sec:node}, we define the key unit cell or node for the vector unit, specifically the self-configuration or self-error correction step, which automatically corrects for any fabrication and/or hardware error without requiring off-chip calculations \cite{Miller2013Self-configuringInvited, Bandyopadhyay2021HardwarePhotonics, Hamerly2021AccurateInterferometers}.
    \item In Section \ref{sec:network}, we define the balanced and unbalanced vector units mathematically, laying the foundation for our error analysis.
    \item In Section \ref{sec:error}, we define a second-order perturbation theory error model for general feedforward networks. More strikingly, we show that sensitivity of each phase shifter is proportional to the power going through it, which is the basis of our error model. We can leverage this property to non-invasively monitor powers in any feedforward photonic circuit using binary trees at the input and output of the device.
    \item In Section \ref{sec:results}, we bring together the results of Sec. \ref{sec:network} and Sec. \ref{sec:error} to derive a statistical power model that explains the increased robustness scaling of balanced compared to unbalanced vector units. We also perform simulations that compare various phase and coupling errors in the mesh that agree with the derived scaling properties.
    \item In Section \ref{sec:apps}, we propose a new ``binary tree cascade'' model, a generalization of the the triangular architecture \cite{Reck1994ExperimentalOperator, Miller2013Self-configuringInvited} which can be useful for error-tolerant mode conversion and principal component analysis-based signal processing. The Appendix later expands on the binary tree cascade to propose other new error-tolerant splay architectures and relates our work to existing well-known butterfly architectures.
\end{enumerate}

\section{Component model} \label{sec:node}

Before we discuss a network model, we propose a component model for the individual nodes or building blocks of the network. First we note that a node does not have to be explicitly an MZI, which has been recognized more recently \cite{Bogaerts2020ProgrammableCircuits} with the development of tunable coupling elements, so we need to describe the basic functionality of the node component that abstracts away the MZI functionality. 

\subsection{Photonic node core functionality}

In general the core functionality of a photonic node in the network is:
\begin{enumerate}
    \item The node is a $ 2 \times 2 $ component (2 left ports and two right ports).
    \item Given any input into the left ports of a node, the node arbitrarily redirects the light entirely into either one of its right ports.
    \item By reciprocity, given an input into either right port, the node arbitrarily redirects the light into its left ports in both amplitude and \textit{relative} or \textit{differential} phase.
\end{enumerate}
Here, we refer to left and right ports separately from input and output ports, following the naming convention of Ref. \citenum{Miller2020AnalyzingNetworks} because the above definition assumes that there is a phase shifter present on the left side of the device and that light can enter from either direction. The ratio of normalized powers in the output ports is known as a ``split ratio.'' We define the reflectivity $r = 1 - s$ and transmissivity $s = 1 - r$ respectively as the fractional power in the bar port and cross port from the port in which light is sent in. ``Bar state'' means $r = 1, s = 0$ (all light goes to ``same side'' port, e.g. lower input to lower output port), and ``cross state'' means $s = 1, r = 0$ (all light goes to the ``opposite side'' port, e.g. lower input to upper output port). Ultimately, the main idea of this paper is to work with photonic circuits that connect these nodes up to implement arbitrary multi-waveguide modes, as discussed in Ref. \citenum{Miller2020AnalyzingNetworks}, as well as other photonic mesh feedforward networks \cite{Pai2020ParallelNetwork}.

In this paper, our error models work for two implementations of nodes depicted in Fig. \ref{fig:node}(b): the tunable directional coupler (TDC) node and the Mach Zehnder interferometer (MZI) node \cite{Bogaerts2020ProgrammableCircuits}, the latter of which is already quite ubiquitous in commercial and academic implementations of these systems \cite{Shen2017DeepCircuits, Annoni2017UnscramblingModes, Harris2018LinearProcessors, Taballione20188x8Waveguides}. The theme here is in general that there are two elements to a node: a phase shifter and a tunable split ratio. In the MZI case, the tunable ratio consists of two directional couplers and a phase shifter. In the TDC case, there is an explicit physical mechanism to modulate the splitting ratio of a $2 \times 2$ coupler by directly perturbing the active coupling region, side-stepping potential errors in the individual passive directional couplers of the MZI node, though there are nodes that can be corrected for these errors in triangular and rectangular architectures \cite{Hamerly2021InfinitelyInterferometers}.

\begin{figure*}
    \centering
    \includegraphics[width=\linewidth]{node.pdf}
    \caption{(a) We show a universal $4 \times 4$ photonic network that arises from the well-known Reck architecture. (b) The individual nodes of the universal network are either MZIs or TDCs with errors in phase shifts and coupling. A benefit of the TDC is that it does not rely on fundamental limits based on errors in the fixed splitters $\delta$. (c) Self-configuration of a universal network of four inputs consists of three steps, where each step involves self-configuration to encode a row vector of $U$, a $4 \times 4$ unitary device operator. (d) Self-configuration of a 4-tree vector unit involves nullifying columns in parallel, similar to the proposals in Refs. \cite{Miller2020AnalyzingNetworks, Pai2020ParallelNetwork}.}
    \label{fig:node}
\end{figure*}

\subsection{Node transmission matrix}

It is often mathematically convenient to represent the transmission matrix acting on the two element vector representing the input modes approaching from the left side $\bm{x} \equiv (x_0, x_1) = (x_\mathrm{lower}, x_\mathrm{upper})$. When nodes are connected together, this approach ultimately allows us to define a matrix $U$ to represent the transmission matrix. This may sound familiar, as this formulation is a condensed version of the $S$-matrix formulation, which considers both reflection into the input ports and transmission into the output ports. Since we assume the reflection is sufficiently small, the transmission matrix formulation is also sufficient since it considers transmission into the output ports given an input port excitation. Therefore, based on our definition above, the general ideal $2 \times 2$ node matrix is defined in terms of a phase matrix component $P$ and coupling component $X(s)$:
\begin{equation}\label{eqn:node}
    T(s, \phi) \equiv X(s) P(\phi) \equiv  \begin{bmatrix} -\sqrt{1 - s} & \sqrt{s}\\ \sqrt{s} & \sqrt{1 - s} \end{bmatrix} \begin{bmatrix} e^{i \phi} & 0\\ 0 & 1 \end{bmatrix},
\end{equation}
where any node implements some functionally equivalent form of the $XP$ representation above. In this paper, we abstract away the details of exactly how $s$ and $\phi$ are implemented, though generally it may take the form of an MZI or a tunable coupler \cite{Bogaerts2020ProgrammableCircuits}. In such cases, we typically can find a physical parameter that behaves like a phase (e.g. an arm of the MZI) or the inverse beat length (difference between the first and second mode propagation constants in an MMI or directional coupler) such that the transmissivity $s = \cos^2 \frac{\theta}{2}$ for $\theta \in [0, \pi]$. More explicitly, we will also consider the following more standard parametrization which for MZIs, assumes a ``differential mode phase shift'' \cite{Pai2020ParallelNetwork}:
\begin{equation}\label{eqn:mzinode}
    T_{\mathrm{MZI}}(\theta, \phi) \equiv X(\theta) P(\phi) \equiv  ie^{i \frac{\theta}{2}}\begin{bmatrix} -\sin \frac{\theta}{2} & \cos \frac{\theta}{2}\\ \cos \frac{\theta}{2} & \sin \frac{\theta}{2} \end{bmatrix} \begin{bmatrix} e^{i \phi} & 0\\ 0 & 1 \end{bmatrix}.
\end{equation}
A TDC node, on the other hand, has the parametrization, based on tuning a beat length of the directional or MZI:
\begin{equation}\label{eqn:tdcnode}
    T_{\mathrm{TDC}}(\theta, \phi) \equiv X(\theta) P(\phi) \equiv  ie^{i f(\theta)}\begin{bmatrix} \cos \frac{\theta}{2} & i\sin \frac{\theta}{2}\\ i\sin \frac{\theta}{2} & \cos \frac{\theta}{2} \end{bmatrix} \begin{bmatrix} e^{i \phi} & 0\\ 0 & 1 \end{bmatrix},
\end{equation}
where we include a $e^{i f(\theta)}$ overall phase term that depends on the exact TDC design and coupled mode theory \cite{Huang1994Coupled-modeOverview}. Ultimately, we will focus specifically on architectures of MZI nodes, but our analysis can be extended to uncover error models of architectures of TDC nodes as well.

\subsection{Feedforward architectures}

In this paper, we are specifically interested in how error from individual nodes affects the overall error in entire ``feedforward networks'' of nodes given by Eq. \ref{eqn:node}  \cite{Pai2020ParallelNetwork}, as opposed to networks containing no cyclic loops as in Ref. \citenum{Perez2017HexagonalInterferometers}. As previously discussed, binary tree architectures, butterfly/Benes nonlocal architectures, and rectangular or triangular universal architectures all fall under the umbrella of feedforward networks. As with the component definition above, we assume monochromatic light and light propagating from left-to-right in the network.

When a feedforward photonic circuit is operated, an $N$-dimensional input vector (or input data) $\bm{x}$ enters the left waveguide ports of the device and propagates forward (left-to-right) through the MZI network until it reaches the right side of the network where output amplitude and phase $\bm{y}$ is measured at a set of photodetectors \cite{Miller2013Self-configuringInvited}. Since ideally no light is lost in the circuit (the coupling nodes are ideally just rearranging the light), the propagation of the $N$-dimensional input through the network may be modelled as a \textit{unitary} or norm preserving transformation $U$. Thus our measured $\bm{y}$ is a result of a change in mode basis operation performed by the photonic circuit, which can be represented mathematically as a matrix product $\bm{y} = U \bm{x}$. To recover the mode basis, the column vectors of $U$ can be determined by measuring outputs given inputs into individual input waveguides of the circuit. The row vectors are recovered by sending light back in from the right and measuring amplitudes and conjugate phases on the left of the circuit. Note that there is a final set of phase shifts (a tunable phase screen) placed at the end of the feedforward network that can be generally useful to achieve unitary architectures. For simplicity, we generally do not consider the contribution of these phase shifts as they are not required for self-configuration of vector unit architectures \cite{Miller2020AnalyzingNetworks}.


As shown in Fig. \ref{fig:node}(c, d), the property that MZIs can guide any incoming mode into a single waveguide means that the universal triangular architecture in Fig. \ref{fig:node}(a) can self-configure itself to program any unitary operator in an $N \times N$ optical spatially multiplexed system \cite{Miller2013Self-configuringInvited}. The protocol consists of three steps to self-configure a $4 \times 4$ unitary $U$, where each step self-configures a vector unit of successively smaller inputs (4, 3, 2). This protocol is useful because it can also correct for any fabrication errors in the device, owing to the model-free optimizations of phase shifts at each step of the process, discussed in further detail in the Appendix. Fig. \ref{fig:node}(d) suggests that depending on the architecture, self-configuration can be completed in fewer steps (here shown for a ``balanced'' vector unit example).

As it pertains to this paper, self-configuration uses model-free feedback optimizations to automatically correct for any error while programming, a convenient property we will use later to analyze error-corrected bounds of binary tree architectures. These errors appear in the form of phase errors (wavelength tuning or environmental perturbations), coupling errors and optical losses due to fabrication process variation. In Refs. \cite{Bandyopadhyay2021HardwarePhotonics, Hamerly2021InfinitelyInterferometers}, strategies for error tolerance are proposed that enable error correction in the presence of realistic errors in the splitters, which can actually be mapped to corrected phase errors. These issues are potentially avoided by using a tunable TDC node because such nodes include couplers that are nominally tunable from cross to bar state. We will focus specifically on modelling the MZI node in this paper for simplicity, because as previously mentioned, we avoid a global phase term in the TDC definition that varies with $\theta$ and is thus less straightforward to model. However, assuming such a model is found, similar concepts can be transferred from our analysis in this paper.

\section{Binary tree networks} \label{sec:network}

\begin{figure*}
    \centering
    \includegraphics[width=\textwidth]{tree.pdf}
    \caption{(a) Generator symbol for $N = 8$. (b) Unbalanced generator network for $N = 8$. (c) Balanced generator network for $N = 8$. (d) Recursive generator network definition to generate any vector unit structure. (e) Analyzer symbol for $N = 8$, including nullified ports indicated by purple crosses. (f) Unbalanced analyzer network for $N = 8$. (g) Balanced analyzer network for $N = 8$. (h) Recursive analyzer network definition to generate any vector unit structure. Note: all red dotted lines and red labels denote the index assignments for the individual nodes (both the $\theta_n, \phi_n$ phase shifters in those nodes) compatible with the definition in our Simphox framework \cite{Pai2022Simphox:Library}.}
    \label{fig:tree}
\end{figure*}


In this section, we describe a particular implementation of $U$ that allows for self-configuration. For this, we assume lossless and ideal nodes as described in Section \ref{sec:node}.

\subsection{Embedding nodes in a larger circuit}
In an $N$-waveguide photonic circuit that is feedforward, we can consider modes propagating through a circuit along $N$ ``rails'' indexed $1$ to $N$. This framework is defined explicitly in our previous work Ref. \cite{Pai2020ParallelNetwork} and also in the original proposal for photonic networks \cite{Reck1994ExperimentalOperator}.

We now need some formalism to place $2 \times 2$ elements along these rails in a specified order for light to propagate. For this, we define an embedding matrix for $T$ of Eq. \ref{eqn:node} of the form $T_{m, n}^{[N]}$, where modes $m$ and $n$ are the waveguide indices for the modes to interfere and $N$ is the number of waveguides or ``circuit size.'' This is a unitary operator known in the mathematics literature as a \textit{Givens rotation}:

\begin{equation}\label{eqn:givensrotation}
T_{m, n}^{[N]} \,:=\,\, \begin{blockarray}{ccccccccc}
 &  & \small{m} & & \small{n} &  & N & \\
\begin{block}{[ccccccc]cl}
1   & \cdots &    0   &   0  &   0   & \cdots &    0 & & \mybigstrut\\
\vdots & \ddots & \vdots &  \vdots &  \vdots &        & \vdots & &\\
0   & \cdots &    T_{11}  & \cdots &   T_{12}   & \cdots &    0 & & \small{m} \\
0   & \ddots &    \vdots  & \ddots &   \vdots   & \ddots &    0 & &  \\
0   & \cdots &    T_{21}  & \cdots &   T_{22}  & \cdots &    0   & & \small{n}\\
\vdots &   & \vdots  & \cdots & \vdots & \ddots & \vdots & &\\
0   & \cdots &   0 & 0 &    0   & \cdots &    1 & & N\mynegbigstrut\\
\end{block}
 &  & & &  & & \\
\end{blockarray},\\
\end{equation}
where $T$ is defined as in Eq. \ref{eqn:node}.

In summary, we have a unitary matrix that is an identity matrix, with 1's along the diagonal except in row/column $m, n$ where there are also off-diagonal terms $T_{21}, T_{12}$. This deviates from the definition for embedded photonic node in Ref. \citenum{Pai2020ParallelNetwork}, but ultimately simplifies the representation for defining tree networks.

We refer to a sequence of embedded nodes $\{..., T_{m, n}^{[N]}, ...\}$ that commute, or equivalently do not share any common index $m, n$ as a ``level'' or ``column'' of nodes. We showed previously in Ref. \citenum{Pai2020ParallelNetwork} that any feedforward circuit can be decomposed into efficient columns of nodes that can be tuned (or self-configured) in parallel. For example, the nodes $T_{1, 2}^{[N]}$ and $T_{2, 3}^{[N]}$ \textit{cannot} be evaluated simultaneously (do not commute) because they share a common index $2$ and thus are connected to each other differently depending on the specified order, i.e. $T_{1, 2}^{[N]}T_{2, 3}^{[N]} \neq T_{2, 3}^{[N]}T_{1, 2}^{[N]}$. On the other hand, $T_{1, 2}^{[N]}$ and $T_{3, 4}^{[N]}$ operators do commute ($T_{1, 2}^{[N]}T_{3, 4}^{[N]} = T_{3, 4}^{[N]}T_{1, 2}^{[N]}$) and equivalently can be applied simultaneously. Grouping all columns of nodes is critical in any feedforward architecture to build compact designs and representations for self-configuration or calibration, and this can be achieved via topological sorting by time-order traversal \cite{Pai2020ParallelNetwork}.

\subsection{Definition}

Vector units are architectures that encode arbitrary complex vectors up to an overall magnitude. A vector unit is also capable of implementing a single row or column vector of $U$ in a photonic circuit, the first step in generating an architecture that computes $\bm{y} = U \bm{x}$ matrix-vector products. Such vector units can be fully specified using a \textit{binary tree} data structure consisting of nodes of the form in Eq. \ref{eqn:givensrotation}. At a high level, binary tree structure ensures a single waveguide has a connection path to any $N$ waveguides, and can be used to generate either a generator or an analyzer configuration \cite{Miller2020AnalyzingNetworks}. Vector units are currently primarily deployed as arbitrary state preparation architectures as an input to a general photonic circuit. State preparation (generation) is the inverse of self-configuration (analysis). Thus, by reciprocity as in Fig. \ref{fig:tree}(a), a tree network self-configured (programmed) to an $N$-dimensional complex vector $\bm{y}_N \in \mathbb{C}^N$ can be physically flipped (input enters from root node) and used as a generator to implement $\bm{y}_N^*$. In other words, we can send light into the root node of a balanced tree circuit in generator orientation and achieve any prepared state. In this vein, we will refer to the ``orientation'' of a tree network as an ``analyzer configuration'' if the root node is at the output of the network, and as a ``generator configuration'' if the root node is at the input of the network, assuming light always goes from left to right, following the convention of Ref. \citenum{Miller2020AnalyzingNetworks}.

Balanced vector units in photonic circuits generally require $\log_2 N$ layers, and consist of either a generator binary tree \cite{Miller2020AnalyzingNetworks} or a splitting tree with programmable phase and attenuating elements to define the individual elements of the vector. A generator binary tree (as we will define in this paper) is less lossy than the attenuating element version, which can incur an additional loss of up to $-20$ dB for $N = 100$ when ``one-hot" vectors (standard basis vectors) are programmed on the device.

\subsection{Binary tree}
Analyzer-oriented and generator-oriented vector units are ``binary tree'' graph networks decorated with input and output node-edge pairs, so it is useful to define a binary tree in this context. In this paper, we define a ``binary tree'' using the standard constructive definition definition: starting from a root node, we add at most two child nodes and repeat this process on the child nodes which can also be considered binary trees (or ``subgraphs'') of the original tree. If node $A$ is a child of node $B$, then there is an edge from node $B$ to node $A$. Optionally, we may use directed edges in binary trees, which in our context, are equivalent in function to a waveguide with photons or light travelling in the direction of the edge. Any node that has no children is referred to as a ``leaf'' of the binary tree. We refer to an $N$-tree as any tree that has $N - 1$ nodes. For instance, a $2$-tree is just a single root node with no children. We will show later that an $N$-tree can fully parametrize a vector unit that parametrizes any $\bm{y}_N$.

As shown in Fig. \ref{fig:tree}(b, c), there are two structural or graph-topological extremes for defining an $N$-tree. On one extreme, each node has a single child node: this is also referred to as a ``unbalanced'' binary tree network since the resulting structure is a maximally unbalanced tree. On the other extreme, each node above the ``leaf level'' has two children, in which case we get a \textit{balanced} binary tree, which is also known as a ``divide and conquer'' scheme. Note our definition of binary tree diverges from that of Ref. \citenum{Miller2020AnalyzingNetworks}, in that linear chains are also considered to be extreme forms of binary trees and ``hybrid'' architectures are just binary trees of varying structure or graph topology. We will also use the terms ``balanced tree'' and ``unbalanced tree'' in place of ``binary tree'' and ``diagonal line.'' In other words, our use of ``binary tree'' here maps to a broader definition of a binary tree than in Ref. \citenum{Miller2020AnalyzingNetworks}, where we allow varying extremes of balance. 

Binary trees, as defined here, are the minimal representation to define any self-configuring layer. (This representation is consistent with the criteria for self-configuring layers in Ref. \citenum{Miller2020AnalyzingNetworks} (Appendix).) However, such data structures do not consider other types of nodes in the actual photonic circuit that need to be defined such as input and output nodes (representing optical interconnects such as edge or grating couplers), because these nodes are not necessary to define to arrive at a specific MZI network topology. This is because once a binary tree subgraph is defined, we ``fill'' any missing edges and assign them to input (node with outgoing edge) or output (node with incoming edge) node-edge pairs respectively. In particular, all nodes should have a $2 \times 2$ structure (two incoming and two outgoing edges), and the remaining unfilled edges must be ``filled'' by any missing input or output node-edge pairs. In the final circuit graph, all nodes with a parent have one output node-edge pair filled and the root node has two such pairs filled. All nodes with a single child node have missing outgoing edges filled by an input node-edge pair, and all nodes with no children (leaf nodes) have both missing outgoing edges filled by input node-edge pairs.

\subsection{Generator and calibration}
The elements of a generated vector can be expressed in terms of the phases in a recursive manner, where a node $D_N = T_{N_1, N}^{[N]}(\theta_{N_1}, \phi_{N_1})$ splits vector magnitudes across subtrees that have $N_1$ and $N_2$ outputs as follows:
\begin{equation} \label{eqn:generator}
    \begin{aligned}
        \bm{x}_N &:= D_N^\dagger \bm{e}_N = \begin{bmatrix}
               \cos\left(\frac{\theta_{N_1}}{2}\right) e^{i \phi_{N_1}} \bm{x}_{N_1} \\
               \sin\left(\frac{\theta_{N_1}}{2}\right) \bm{x}_{N_2}
        \end{bmatrix},
    \end{aligned}
\end{equation}
where we define the device operator $D_N^\dagger$ with input to the $N$th (``bottom'') input of the generator device, and we prune off $\theta, \phi$ for the root node of $\bm{x}$ at each recursive step. This simple representation gives a direct formula for each vector element and allows for straightforward calculation of error sensitivities with respect to phases in the specific programmed vector. This is also a convenient formula to use when computing the sensitivity model of the architecture in terms of the control voltage of the phase shifters rather than the phase shifts alone.

The generator configuration can also be used to \textit{calibrate} the $s$ split ratio or $\theta$ phase shifters of the device to generate a lookup table for each of the phase shifter settings in the network as a function of voltages or some other control parameter. Such procedures are defined in Ref. \citenum{Miller2020AnalyzingNetworks} for any binary tree network. One key point in the paper is that rearranging phase shifts into a symmetric configuration (with phase shifts in the two arms of the MZI), can allow one to effectively also calibrate the $\phi$ phase shifters via a reparametrization of the system. However, the $\phi$ phase shifters as defined in this paper cannot be calibrated this way. That can instead be achieved using an analyzer configuration and self configuration as we now discuss.

\subsection{Analyzer and self-configuration}
We now consider the self-configuring analyzer architecture, first proposed in Ref. \citenum{Miller2013Self-aligningCoupler}. We can write a more explicit formula for the unitary matrix $D_N$ first used in Eq. \ref{eqn:generator} in terms of products of matrices of the form in Eq. \ref{eqn:givensrotation}.

Shifting from a graph topology to a matrix framework, we define a formula for \textit{any} binary tree resulting in the unitary matrix $D_N$ (unitary binary tree formulation with $N$ outputs) by a recursive definition:
\begin{equation} \label{eqn:binarytree}
    \begin{aligned}
        D_{N}(\mathbf{s}_N,\bm{\phi}_N) &= T_{N_1, N}^{[N]}(s_{N_1}, \phi_{N_1})\begin{bmatrix}
                D_{N_1} & 0 \\
                0 & D_{N_2}
        \end{bmatrix} \\
        \mathbf{s}_N &= [\ldots, s_{N_1}, \ldots] \\
        \bm{\phi}_N &= [\ldots, \phi_{N_1}, \ldots]
    \end{aligned}
\end{equation}
where the only requirement is that $1 < N_1 < N$ and $N_1 + N_2 = N$ ($N_1$ decides how many outputs each subtree gets at each split, the rail index of the top waveguide, and can be considered the index of the node within a given subtree). We refer to $D_{N_1}$ as the ``top subtree'' and $D_{N_2}$ as the ``bottom subtree.'' To label the nodes further down in the tree (the ``descendants'' of node $N_1$), we apply an offset $N_1$ to any node index in the bottom subtree $D_{N_2}$ as part of the recursive step. As a part of defining the block matrix, we also let ``$0$'' denote the setting of all off-diagonal block-matrix elements to zero. In general $\mathbf{s}_N$ contains $N - 1$ values in the range $[0, 1]$ and $\bm{\phi}_N$ contains $N - 1$ values in the range $[0, 2\pi)$. This, along with the overall magnitude or power and the overall phase accumulated, accounts for all the degrees of freedom in a complex vector, thus $2N$ total degrees of freedom to define an $N$-dimensional complex vector.

Any self-configuring graph architecture may be defined by specifying $N_1$ for recursively defined subtree branches until the base case $D_1 = 1$ is reached. For a (``fully'') balanced tree, we have $N_1 = \lfloor \frac{N}{2} \rfloor$ (i.e., the ``floor'' or integer part of $\frac{N}{2}$), and for a (``fully'') unbalanced tree (chain or ``diagonal line'') we have $N_1 = N - 1$. Variations on the binary tree structure are specified by recursively specifying $N_1$ starting from the root node into repeated invocations of Eq. \ref{eqn:binarytree} until the base case is reached. We will mostly be concerned with the extreme two cases in this paper, because they exhibit key differences in how they scale, but arbitrary structures may also be defined, whose structures fall in the spectrum between balanced and (fully) unbalanced.

A proof of self-configuration can also be done straightforwardly using recursion and the methods of Eqs. \ref{eqn:phisc} and \ref{eqn:ssc}. Consider a lab setting where we are given some complex vector $\bm{x}_N \in \mathbb{C}^N$ and we do not know anything about the settings of the network implementing $D_N$. Through successive minimizations of powers in detectors, we can experimentally self-configure the device so that all the light exits the bottom waveguide indexed at $N$. In mathematical terms $D_N \bm{x}_N = \widetilde{x}_{N} \bm{e}_N$, where $\bm{e}_N$ is the $N$th standard basis vector in $\mathbb{C}^N$ and $\widetilde{x}_{N}$ represents some  phase and amplitude output by the self-configuration with $|\widetilde{x}_{N}|^2 = \|\bm{x}\|^2 \equiv P_N$ (power conservation). Along with these equations, we can use induction to prove that Eq. \ref{eqn:binarytree} implements self-configuration. Specifically, we first split $\bm{x} = [\bm{x}_{N_1}, \bm{x}_{N_2}]$ and apply our inductive hypothesis that $D_{N_1} \bm{x}_{N_1} = \widetilde{x}_{N_1}\bm{e}_{N_1}$ and $D_{N_2} \bm{x}_{N_2} = \widetilde{x}_{N_2}\bm{e}_{N_2}$. Note that power conservation ensures $P_N = P_{N_1} + P_{N_2}$, i.e. the total power in the first and second subtrees of the overall binary tree is always additive (no power is lost or gained).

To complete the proof, we note that we can apply the procedure in the methods of Eqs. \ref{eqn:phisc} and \ref{eqn:ssc} to program the matrix $T_{N_1, N}^{[N]}$ to achieve the result:
\begin{equation} \label{eqn:selfconfiguringproof}
    \begin{aligned}
        D_N \bm{x}_N &= T_{N_1, N}^{[N]}(D_{N_1}, \phi_{N_1})\begin{bmatrix}
                D_{N_1} & 0 \\
                0 & D_{N_2}
        \end{bmatrix} \bm{x}_N\\
        &= T_{N_1, N}^{[N]}(s_{N_1}, \phi_{N_1}) (\widetilde{x}_{N_1}\bm{e}_{N_1} + \widetilde{x}_{N_2}\bm{e}_{N}) \\
        &\equiv  \widetilde{x}_{N} \bm{e}_{N}.
    \end{aligned}
\end{equation}

Note that the base case for the inductive proof is trivial: $D_1 \bm{x}_1 = \bm{x}_1 = \widetilde{x}_{1} \bm{e}_{1}$ since $\bm{e}_{1} = U_1 = 1$. We have now proven that using our recursive definition, any binary tree (not just ``fully" balanced or ``fully" unbalanced) can be self-configured according to the physical process outlined in Eqs. \ref{eqn:phisc} and \ref{eqn:ssc}.

Now that we have proven self configuration is always possible for a (general) binary tree network, such a network can be parametrized in one of two ways: in terms of node settings (which we have defined above) or in terms of the inputs into the layer $\bm{u}_N$ satisfying the condition $D_N(\mathbf{s}_N, \bm{\phi}_N) \bm{u}_N = \bm{e}_{N}$. To define the latter, we define a new matrix function $R_N(\bm{u}_N) = D_N(\mathbf{s}_N, \bm{\phi}_N)$, where $\bm{u}_N \in \mathbb{C}$ satisfies $\|\bm{u}_N\|^2 = 1$ and $\mathrm{arg}(u_1) = 0$ so that the remaining $2N - 2$ degrees of freedom match the number provided by $N - 1$ nodes ($s, \phi$ for each node).

Finally, the time to self-configure relies on the number of columns in the vector unit since MZIs in a column can be self-configured simultaneously \cite{Pai2020ParallelNetwork}. The self-configuration time would need to be sufficiently small to adjust to any incoming training signal, e.g. for a sensor reading modulated optical modes. The number of columns in a binary tree is given by $\log_2 N$, whereas the number of columns in a fully unbalanced tree is given by $N$, so therefore self-configuration is faster by an order of $N / \log_2 N$, which grows quickly with $N$. Hereafter we will use $\log N$ to refer to $\log_2 N$ for ease of notation.

\section{Error model} \label{sec:error}

In this section, we set up the core results of this paper by deriving phase error sensitivity properties for feedforward photonic networks rigorously defined using the ``Hessian'' of a least squares error function describing the intended and measured behavior of the devices. This more general result is needed to ultimately compare the overall performance of various vector unit architectures, including both phase and coupling errors. Unlike in previous sections, here we assume that the feedforward network is already programmed to some desired setting of phase shifts, so the actual implemented unitary operator on the device is already intended to be some $U$. Given this ideal state, we would like to analyze some perturbed ``error'' state $\widehat{U}$ where we vary one or two of the phases.

\subsection{Error function and Hessian}

Here, we proceed to define architecture-dependent ``error sensitivities'' of vector units. We define the error sensitivity of a component to be a ratio between the overall mean square error (due to that component alone) and component error. Here, we seek to describe how various errors (other than loss) in various elements of the circuit affect the \textit{overall} circuit error. Such errors (e.g., fabrication and environmental errors) ultimately can be modelled as coupling errors and phase errors; since coupling errors can be effectively reduced to phase errors in an ideal node \cite{Bandyopadhyay2021HardwarePhotonics}, the phase error can effectively be used to describe either.

We begin by defining the mean square circuit error: 
\begin{equation} \label{eqn:vectorerror}
    \epsilon^2 = \|\bm{y} - \hat{\bm{y}}\|^2= 2 - 2\mathcal{R}(\bm{y}^\dagger \hat{\bm{y}}),
\end{equation}
where $\hat{\bm{y}}$ is the measured output vector, $\bm{y}$ is the predicted output vector, and $\mathcal{R}$ denotes taking the real part. Note again that we are now dealing with the output of a general feedforward mesh.

The first step in calculating the sensitivity of any device is to realize the gradient of the error function in Eq. \ref{eqn:vectorerror} is zero when $\bm{y} = \hat{\bm{y}}$, hence the use of a second-order term called a ``Hessian'' to describe the errors in the circuit. Given that we are aiming for generality, note that this holds regardless of whether the intended input/output behavior $\bm{x}, \bm{y}$ is known a priori. We could in principle run an experiment using an ideal analyzer and generator to determine the full matrix $U$ implemented by the device without any knowledge of the internal settings of the device (assuming Hermitian operator). However, the internal settings of the device, whatever they are, will still have some sensitivity based on the \textit{current} implemented settings.

Assume that the vectors $\bm{\theta}, \bm{\phi}$ are the true phases and the vectors $\widetilde{\bm{\theta}}, \widetilde{\bm{\phi}}$ are the phases with error. We define the error vector $\bm{\Delta} := [\bm{\Delta}_{\theta}, \bm{\Delta}_{\phi}] := [\bm{\theta} - \widetilde{\bm{\theta}}, \bm{\phi} - \widetilde{\bm{\phi}}] := \bm{\eta} - \bm{\eta}'$. This gives the following expression for the error in the network phases:
\begin{equation}\label{eqn:hessian}
    \begin{aligned}
        \epsilon^2(\bm{\Delta}) &= \cancelto{0}{\epsilon^2(\bm{0})} + \bm{\Delta}^T \cancelto{0}{\frac{\partial \epsilon^2}{\partial \bm{\Delta}}} + \frac{1}{2}\bm{\Delta}^T \mathcal{H}_{\epsilon^2} \bm{\Delta} + \cdots \\
        &= \frac{1}{2}\bm{\Delta}^T \mathcal{H}_{\epsilon^2} \bm{\Delta} \\
        &:= \frac{1}{2}[\bm{\Delta}^T_\theta, \bm{\Delta}^T_\phi]\begin{bmatrix}
               \mathcal{H}_{\theta \theta} & \mathcal{H}_{\theta \phi} \\
               \mathcal{H}_{\phi \theta} & \mathcal{H}_{\phi \phi}
        \end{bmatrix} \begin{bmatrix}
               \bm{\Delta}_{\theta} \\
               \bm{\Delta}_{\phi}
        \end{bmatrix},
    \end{aligned}
\end{equation}
where as we have just claimed, the first order gradient term evaluates to zero leaving us with a Hessian $\mathcal{H}_{\epsilon^2}$ where the first rows and columns specify all $\theta$ phase shifts and the final rows and columns specify all $\phi$ phase shifts.

We now can determine the matrix elements of $\mathcal{H}_{\epsilon^2}$, which describes how both the individual phase shifter sensitivities and correlations among phase shifters contribute to the overall error. Note that based on the properties of the second-order derivative terms, we have $\mathcal{H}_{\phi\theta} = \mathcal{H}_{\theta\phi}^T$.

For phase shifters indexed at $i, j$ (corresponding to the concatenated $\bm{\phi}, \bm{\theta}$ vectors $\bm{\eta}$), we can write the formula for each element of the Hessian. We will also need to compare to simulation. To evaluate the Hessian given errors $\delta_i = \eta_i - \hat{\eta_i} \ll \eta_i$, the Hessian can be written in terms of central finite difference (which is required to compute the diagonal terms correctly):
\begin{equation}\label{eqn:hessianfd}
    \begin{aligned}
        \mathcal{H}_{ij} &:= \frac{\partial^2\epsilon^2}{\partial\delta_i \partial\delta_j} \approx \frac{\epsilon^2(\delta_i\bm{e}_i + \delta_j\bm{e}_j) - \epsilon^2(\delta_i\bm{e}_i - \delta_j\bm{e}_j)}{2\delta_i\delta_j}.
    \end{aligned}
\end{equation}
Note that to compute the Hessian here we are subtracting error contributions where perturbations go in the same direction from those where they go in the opposite direction. Note that if $i = j$, the first term is nonzero and the second term is zero.

However, we can also use the generator formula of Eq. \ref{eqn:generator} to derive an exact formula for the Hessian to avoid needing to perform highly computationally intensive tasks as Hessian finite differences in Eq. \ref{eqn:hessianfd}. We will now explicitly evaluate the Hessian matrix elements, and in the process derive the sensitivity and correlations across many phase elements in any given feedforward mesh network and expected input/output mode pair.

\begin{figure*}
    \centering
    \includegraphics[width=\linewidth]{binarytree_hessianmeasurement.pdf}
    \caption{We show how we might experimentally implement a phase sensitivity or power monitor in a feedforward mesh (green) using an error detection circuit (blue and red). The error detection circuit assumes a perfectly ideal input generator (red) and output (blue) analyzer set to input/output mode pair $\bm{x}, \bm{y}$ respectively. As shown in the ``real part'' inset in (a), we also use an extra reference path to interferometrically determine the real part of the output signal corresponding to our mean square fidelity (or $1 - \epsilon^2$) assuming equal amplitudes at the input between reference and mesh signal (This can be adjusted using attenuation on the reference path.). (a) To measure power or sensitivity at any phase shifter, we simply perturb a phase shifter $\eta$ by some fixed amount $\delta \eta$ in the middle of the mesh and measure the corresponding change in the output power. (b) Measuring the Hessian in any feedforward mesh requires perturbing a descendant $n'$ and ancestor $n$ phase shift by the same amount and subtracting the  response resulting from going in the same and opposite directions. If there are no paths between the phase shifts, the Hessian contribution is zero, and if they are the same phase shift or $n = n'$, reduce to the case in (a). (c) The binary tree mesh problem that we consider involves sending just a single mode into a generator binary tree and evaluating the Hessian on the output vector.}
    \label{fig:sensitivityhessian}
\end{figure*}

\begin{figure*}
    \centering
    \includegraphics[width=\linewidth]{hessian_with_examples.pdf}
    \caption{Assuming an MZI node as in Fig. \ref{eqn:node}, we calculate relative magnitudes of balanced (a) and unbalanced (b) Hessian terms based on Eqs. \ref{eqn:hessianbt}, \ref{eqn:hessianpower}. These results clearly depict larger correlations in unbalanced architectures and also larger sensitivities for the individual nodes (along the diagonal). Next to each plot, we include example diagrams of the correlation links between nodes lined up with the quadrants of the Hessian $\mathcal{H} = \begin{bmatrix} \mathcal{H}_{\theta \theta} & \mathcal{H}_{\theta \phi} \\ \mathcal{H}_{\phi \theta} & \mathcal{H}_{\phi \phi} \end{bmatrix}$ so it is possible to reason out why certain Hessian terms are zero (nodes not connected) or why some are larger than others (more expected power in the descendant node).}
    \label{fig:hessian}
\end{figure*}

\subsection{Sensitivity in feedforward networks}

We want to prove the following claim: ``Given some phase shifter $\eta$, a photonic feedforward mesh implementing $U$, and the input/output mode pair $\bm{x}, \bm{y}$, the sensitivity of any individual phase shifter $\eta$ in the device is equal to the power going through that phase shifter $p_\eta$ given a linear square error function as in Eq. \ref{eqn:vectorerror}.'' Note that we are given a single input and output mode pair because Eq. \ref{eqn:vectorerror} considers a single vector error, but we will consider a more general case later. Additionally, note that in previous scenarios, we considered the specific case where the desired $\bm{y} = \bm{e}_N$ ($N$th standard basis vector) for a vector unit successfully programmed to ideal $\bm{x}$. More generally, given any input $\bm{x}$ we attain some ideal outcome $\bm{y} = U \bm{x}$.

We define operators such that $U = R_\eta P_\eta L_\eta$, where $L_\eta, R_\eta$ represent operators before and after (to the left and right) of the phase shifter $\eta$ in a given device and $P_\eta$ is the operator for some applied $\eta$ phase shift, i.e., a diagonal unitary matrix where a $e^{i \eta}$ phase shift is applied to any single waveguide mode of the system as shown in Fig. \ref{fig:sensitivityhessian}(a).

Substituting into Eq. \ref{eqn:vectorerror}, we calculate the error of that phase shifter for some $\hat{\eta} = \eta + \delta \eta$ error:
\begin{equation} \label{eqn:powererror}
    \begin{aligned}
        \epsilon^2(\delta \eta) &= 2 - 2 \mathcal{R}(\bm{y}^\dagger \hat{\bm{y}})
        \\ &= 2 - 2 \mathcal{R}(\bm{x}^\dagger L_\eta^\dagger P_{-\eta} \cancelto{1}{R_\eta^\dagger R_\eta} P_{\hat{\eta}} L_\eta \bm{x}) \\
        &= 2 - 2 \mathcal{R}(\bm{x}^\dagger L_\eta^\dagger P_{\delta \eta} L_\eta \bm{x}) \\
        &= 2 - 2 \mathcal{R}(\bm{y}_\eta^\dagger P_{\delta \eta} \bm{y}_\eta) \\
        &:= 2 (1 - p_\eta \cos\delta \eta) \approx p_\eta \delta \eta^2,
    \end{aligned}
\end{equation}
where $\delta \eta = \hat{\eta} - \eta$, $p_\eta$ is the relative power in the phase shifter $\eta$ and $\bm{y}_\eta = L_\eta \bm{x}$ is the vector preceding that phase shifter as light propagates through the feedforward device. Again, note that this Eq; \ref{eqn:powererror} holds for any feedforward programmable optical device given any individual input/output mode pair $\bm{x}$ and $\bm{y} = U \bm{x}$. It is important to note that while Eq. \ref{eqn:powererror} holds specifically for $\bm{x}, \bm{y}$, but all phase shifters contribute equally to error in the overall matrix $U$ \cite{Bandyopadhyay2021HardwarePhotonics} as we later address. However, since the emphasis in this paper is that we care about only a subset of modes and not the full Hilbert space spanned by the rows of $U$, the expression of Eq. \ref{eqn:powererror} is of increased importance.

To help with understanding this concept, we show that it is possible to apply the results of Eq. \ref{eqn:powererror} to perform a direct measurement of the sensitivity (and thus monitor intermediate powers) in an arbitrary feedforward network in a direct way experimentally. This is shown diagrammatically in Fig. \ref{fig:sensitivityhessian}(a). We dump half of the light into a reference path between the original generator for $\bm{x}$ and the output of the final analyzer for $\bm{y}$ measured at the output of our feedforward mesh. We can define the ``real'' part of a signal as the contribution in the top waveguide of the final 50/50 coupler. This trivially follows from the definition of a 50/50 beamsplitter matrix when the output mode of the analyzer and the reference path light are out-of-phase by $\pi / 4$.

\subsection{Hessian sensitivity for vector units}

As shown in Fig. \ref{fig:sensitivityhessian}(b) and (c), it is possible to specify the various Hessian matrix elements by directly evaluating correlations betweden perturbations across phase shifts as they affect the mean square error. The sensitivities given by $\epsilon^2(\delta \eta)$ of Eq. \ref{eqn:powererror} are equivalent to diagonal terms $\mathcal{H}_{\eta \eta} / 2$ for the Hessian for any feedforward network. For binary tree vector unit architectures, we now consider the more general case of Hessian sensitivities that involve correlations across different phase shifters rather than just the individual phase shifters, i.e., where $\eta \neq \eta'$.

The Hessian off-diagonal terms relate to correlated errors while the on-diagonal terms relate to uncorrelated sensitivities as in Eq. \ref{eqn:powererror}. More explicitly, we can calculate the total error due to error in the individual phase shifters as shown previously in Fig. \ref{fig:node}(b) as:

\begin{equation} \label{eqn:errorcorr}
    \begin{aligned}
        \epsilon^2(\bm{\Delta}) &\approx \frac{1}{2}\bm{\Delta}^T \mathcal{H}_{\epsilon^2} \bm{\Delta} = \frac{1}{2}\sum_{\eta, \eta'}\mathcal{H}_{\eta \eta'} \delta\eta \delta\eta' \\&= \underbrace{\frac{1}{2}\sum_{\eta}\mathcal{H}_{\eta \eta} \delta\eta^2}_{\mathrm{uncorrelated}} + \underbrace{\frac{1}{2}\sum_{\eta \neq \eta'}\mathcal{H}_{\eta \eta'} \delta\eta \delta\eta'}_{\mathrm{correlated / bias}}
    \end{aligned}
\end{equation}

The statistics of the individual phase errors may be characterized using a covariance matrix, which give the distributions of $\delta\eta^2$ and $\delta\eta \delta \eta'$ and may be measured experimentally. If we assume Gaussian error or noise $\delta \eta \sim \mathcal{N}(0, \sigma_\eta)$, then $\mathbb{E}[\delta \eta^2] = \sigma_\eta^2$, where $\mathbb{E}$ represents the average or expected value. The correlation of a pair of phase errors $\delta \eta, \delta \eta'$ is given by $\mathbb{E}[\delta \eta \delta \eta']$, which is 0 only if the errors are both completely uncorrelated and centered at zero. When using an MZI node, bias may be removed such that the errors are zero by adjusting the input wavelength until the bias is zero (i.e. the expected MZI coupling is 50/50), ensuring $\mathbb{E}[\delta \eta \delta \eta'] = 0$.

In practical linear photonic network implementations, the reason the entire Hessian needs to be considered rather than just the on-diagonal uncorrelated errors is that crosstalk and wavelength errors in phase shifters or couplers might be correlated. Measuring the bandwidth of a photonic network, for example, relies on the measurement of correlated error across phase shifters in the photonic circuit. Additionally, lithography does in many cases introduce spatially correlated errors (e.g. errors in waveguide widths and coupling gaps that are spatially closer on the photonic circuit will have more correlation than those far away).

The Hessian off-diagonal terms require two pieces of information: the ancestor and the descendant phase shifts $\eta', \eta$. If the ancestor and descendant phase shifts are not connected, then there is no Hessian contribution because those phase shifts do not affect each other (This can be also deduced from explicit evaluation of the second derivative of the error with respect to the two phase shifts.). If they \textit{are} connected it should be evident that error in one phase shift will influence the error in the other phase shift.

The computation of the off-diagonal Hessian terms is more straightforward for binary tree vector units compared to other feedforward architectures that do not obey the tree property. Examples of ancestor and descendant ``error flow'' for balanced and unbalanced trees, are provided in our evaluation of the Hessian terms for 8-tree vector units (i.e., $N = 8$) in Fig. \ref{fig:hessian}(a) and (b). In the example 8-tree diagrams, we denote connected phase shifters of types $\theta \to \theta$, $\theta \to \phi$, $\phi \to \theta$, and $\phi \to \phi$ each of which account for a matrix element in each quadrant of the Hessian matrix. Following the same theme as our individual phase shifter sensitivities, the power in the descendant node is the only quantity needed for the relevant Hessian elements.

For any binary tree, the Hessian terms of Eq. \ref{eqn:hessian} given any phase shifter $\eta$ have different expressions based on whether the ancestor is an internal phase shift $\theta$ or external phase shift $\phi$:
\begin{equation} \label{eqn:hessianbt}
    \begin{aligned}
        \mathcal{H}_{\eta \eta} &= 2 p_\eta \\
        \mathcal{H}_{\theta \eta} &= \mathcal{H}_{\eta \theta} = p_\eta \\
        \mathcal{H}_{\phi \eta} &= \mathcal{H}_{\eta \phi} = 2 p_\eta
    \end{aligned}
\end{equation}
In words, ``each element of the Hessian matrix is nonzero if the phase shifters are connected, is proportional to the power going through the descendant phase shifter closer to the leaves of the tree with scaling of 1 if the ancestor is internal phase shifter and 2 if ancestor is an external phase shifter.'' This Hessian definition gives a complete picture of error flow and is ultimately a major step towards efficient and thorough error modeling of binary tree photonic architectures. Further details are provided in the Appendix.

\begin{figure*}
    \centering
    \includegraphics[width=\linewidth]{treestats.pdf}
    \caption{(a) Average distribution of power in a balanced binary tree. (b) Gamma-beta flow is a diagram labelling the statistics of $p_n$ and $s_n$ the power entering node $n$ and the transmissivity of node $n$. The statistics of the programmed analyzer settings are beta distributions ($\beta_{N_1}^{N_2}$) given gamma-distributed input powers ($\gamma_{N'}$). (c) Relative power distributions enforce the constraint that the total power in the system is 1, so the measured powers in each of the waveguides follow beta distributions ($\beta_{N'}^{N - N'}$) rather than gamma distributions ($\gamma_{N'}$). (d) Gamma distributions given an input power $p$. (e) Beta distributions given an input transmissivity or relative power $s$. (f) Phase-beta distributions (for the internal phase shifts).}
    \label{fig:treestats}
\end{figure*}


\begin{figure*}
    \centering
    \includegraphics[width=\textwidth]{error.pdf}
    \caption{Here, we analyze correlated and uncorrelated error, in the presence and lack of error correction for balanced and unbalanced photonic vector units. Legends for each column of this figure are provided in the final row, specifying the various dimensions of our analysis. (a) Random coupling error for unbalanced and balanced trees from $N = 64 \to 4096$. (b) Constant coupling error for unbalanced and balanced trees from $N = 64 \to 4096$. (a) Random phase error for unbalanced and balanced trees from $N = 64 \to 4096$. (b) Constant phase error for unbalanced and balanced trees from $N = 64 \to 4096$. (e) After error correction, the unbalanced architecture now has error proportional to roughly $N \sigma^2$ for large $N$ and $N \sigma^3$ for small $N$. (f) After error correction, the balanced architecture now has error proportional to roughly $\log N \sigma^2$ for large $N$ and $\log N \sigma^3$ for small $N$.}
    \label{fig:error}
\end{figure*}

\section{Results} \label{sec:results}

We now have a prescription for determining the various sensitivities of any feedforward vector unit regardless of structure, but we still need to perform a fair comparison of robustness among various architecture choices. In this section, we define first a new statistical model of tree vector units which ultimately allows us to perform a thorough comparison of the robustness of varying photonic architectures. All calculations are performed in our simulation framework Simphox \cite{Pai2022Simphox:Library}.

\subsection{Statistical model of tree vector units}

We begin by defining how assumed random input distributions transform into powers in the binary tree architectures used to compute the Hessian as in Eq. \ref{eqn:hessian} and resulting in the average behavior of Fig. \ref{fig:hessian}(a) versus (b) (balanced versus unbalanced).

Recently, error models have been proposed for rectangular and triangular architectures under the assumption that inputs to such networks are complex normal vectors \cite{Russell2017DirectMatrices, Bandyopadhyay2021HardwarePhotonics}, which accounts for most realistic scenarios for these devices. We approach the problem from a similar angle, but with more of a statistical focus as we attempt to write a similar framework for the class of binary tree vector units, a larger class of architectures. We ultimately show that maximally balanced trees are the most robust to coupling matrix errors, both analytically and via simulation. Note we use the term ``maximally balanced'' to account for cases a binary tree is never entirely balanced, i.e. where $N \neq 2^K$ for nonnegative integer $K$.

Keeping things simple to start, consider the distribution of equal powers in what we may call the ``average'' case as shown in Fig. \ref{fig:treestats}. The average input has all equal power magnitudes entering the network. Depending on the structure of the network, for an analyzer to route the equally distributed powers to a single output, the appropriate splitting ratios must be defined at each node $n$. We need to come up with a theory to specify what these splitting ratios are as well as what happens when we no longer obey this average simple case. For this reason, we need to make some assumptions about the input distribution into the network, and this will allow us to fairly compare the error tolerances of balanced and unbalanced tree vector units.

The (standard) assumption that inputs to the network are complex normal random numbers means that the inputs are of the form $\bm{x} = \bm{a} + i \bm{b} \sim \mathcal{C}\mathcal{N}(0, 1)$, where $\mathcal{C}\mathcal{N}$ represents a complex normal distribution defined where $a_n, b_n \sim \mathcal{N}(0, 1 / 2)$ are independently distributed complex normal distributions centered at 0 with variance $1 / 2$:

\begin{equation}
    \mathcal{P}(x_n) = \frac{e^{-a_n^2 - b_n^2}}{\pi / 2}
\end{equation}

In physics, we typically work with the phasor representation $x_n = \sqrt{y_n} e^{i \varphi}$, where $y_n = a_n^2 + b_n^2$, so the above probability distribution becomes much simpler:
\begin{equation} \label{eq:exponential}
    \mathcal{P}(y_n, \varphi) = e^{-y_n}
\end{equation}
This proves that that $\varphi \sim \mathcal{U}(0, 2\pi)$ is a uniform distribution, i.e. all values between $0$ to $2\pi$ are equally likely. Additionally, the powers are exponentially distributed which by definition obeys the Gamma distribution defined as $y_n \sim \mathrm{Gamma}(1)$. More generally, the formula for a gamma-distributed random variable distributed as $\mathrm{Gamma}(N)$ is
\begin{equation} \label{eq:gamma}
    \mathcal{P}_\Gamma(y; N) := \frac{y^{N - 1}e^{-y}}{\Gamma(N)},
\end{equation}
where substituting $N = 1$ gives the exponential distribution in Eq. \ref{eq:exponential}. 

We can extend the definition in Eq. \ref{eq:gamma} even further. Specifically, we can use the convenient property that Gamma distributions, like powers in self-configured nodes, are additive. In particular, the average and variance of the distribution for $\mathrm{Gamma}(N)$ is both $N$. The sum of gamma-distributed powers add up both in average and variance to a total gamma distributed power, an interesting statistical analog of energy conservation. 

More formally in our scenario, the total power in a self-configured branch of the binary tree with $M$ inputs is distributed as $P_M \sim \mathrm{Gamma}(M)$ because it is the sum of identically distributed inputs distributed as $\mathrm{Gamma}(1)$. In Fig. \ref{fig:treestats}(b), we label edges (waveguides) with the distribution $\mathrm{Gamma}(N)$ as $\gamma_N$ edges, which gives a shorthand notation to denote the distribution of powers expected in that waveguide given a self-configured vector. The inputs are labelled $\gamma_1$ as expected by the random complex vector condition.

The final step is to find the distribution for $s, \phi$ in the binary tree nodes. A reparametrization is needed using the recursive update in Eq. \ref{eqn:selfconfiguringproof} from $(\widetilde{x}_{N_1}, \widetilde{x}_{N_2})$ to $(s, \phi, \widetilde{x}_{N})$. As we previously showed in Eq. \ref{eqn:ssc}, we can parametrize $s$ as:
\begin{equation}
    s = \frac{|\widetilde{x}_{N_1}|^2}{|\widetilde{x}_{N_1}|^2 + |\widetilde{x}_{N_2}|^2} = \frac{P_{N_1}}{P_{N_1} + P_{N_2}},
\end{equation}
where $P_{N_1} \sim \mathrm{Gamma}(N_1)$ and $P_{N_2} \sim \mathrm{Gamma}(N_2)$. Recall that $s \in [0, 1]$, so we also have $1 - s \in [0, 1]$. As before, we also have $\phi \sim \mathcal{U}(0, 2\pi)$.

The distribution for $s$ is determined using a change-of-basis, and is known as a \textit{beta distribution}, written as $s \sim \mathrm{Beta}(N_1, N_2)$. In statistics more generally, beta distribution can be thought of as a way to measure fairness of a coin given $N_1$ head trials and $N_2$ tail trials. The smaller $N_1, N_2$ are, the less certain we are of the fairness and the larger the variance of the corresponding beta distribution. Here, we propose a new analogy to optical power statistics; we apply the same concept to our physical platform where coin flip probabilities are instead represented as fractions of powers in various segments of the circuit. The beta distribution (also shown in Fig. \ref{fig:treestats}(e)) is defined as:
\begin{equation} \label{eqn:betadist}
\begin{aligned}
    \mathcal{P}_\mathrm{B}(s; N_1, N_2) &= \frac{s^{N_1 - 1} (1 - s)^{N_2 - 1}}{\mathrm{B}(N_1, N_2)} \\
    \mathrm{B}(N_1, N_2) &= \frac{\Gamma(N_1 + N_2)}{\Gamma(N_1) \Gamma(N_2)},
\end{aligned}
\end{equation}
where $\mathrm{B}(N_1, N_2)$ is just a normalization function for the beta distribution (similar to $\Gamma(N)$) that depends on the parameters of the beta distribution. Analogous to coin fairness in our above example, the beta distribution parameters $N_1, N_2$ tell us the average or  ``expected'' fraction of power expected for random variables, which is generally $\langle s\rangle = N_1 / (N_1 + N_2)$. We can further the analogy by relating coin fairness to the expected power fraction. Just as more trials increase coin fairness confidence, an increase in $N_1$ and $N_2$ (the number of inputs leading into the first and second subtrees of a node) corresponds to increased confidence (decreased variance) in allocating power to each subtree, with variance given by $\frac{N_1N_2}{(N_1 + N_2)^2 (N_1 + N_2 + 1)}$. If $N_1 = N_2 = N / 2$ (maximally balanced case), the variance is $\frac{1}{4(N+1)}$. If $N_1 = N - 1, N_2 = 1$ (maximally unbalanced case), then the variance is  $\frac{N - 1}{4N^2(N+1)}$, which is roughly a factor of $N$ less than the variance for the maximally balanced case. The increased confidence in the fraction to attribute to each subtree also decreases the error tolerance for $s$. Each node may be labelled with the notation $\beta_{N_2}^{N_1}$, which indicates that there were $N_1$ inputs that went into the top subtree and $N_2$ inputs that went into the bottom subtree, and that the results of these inputs are now being funneled or combined into this node. In summary, given input powers into a self-configured node labelled $\gamma_{N_1}, \gamma_{N_2}$ for the top and bottom, the node will be labelled $\beta_{N_2}^{N_1}$ and the output edge labelled $\gamma_{N_1 + N_2} = \gamma_N$. This process is applied recursively for $N = 8$ in Fig \ref{fig:treestats}(b, c). We also consider the change-of-variable $s \to \theta$ in Fig. \ref{fig:treestats}(f) for internal phase shifts $\theta$ based on the formula we derived earlier $s = \cos^2 \frac{\theta}{2}$:
\begin{equation} \label{eqn:phasebeta}
    \mathcal{P}_{\mathrm{B}, \theta}(\theta; N_1, N_2) = \frac{\left(\sin \frac{\theta}{2}\right)^{2N_1 - 1} \left(\cos \frac{\theta}{2}\right)^{2N_2 - 1}}{\pi \mathrm{B}(N_1, N_2)},
\end{equation}
which, in the case of $N_2 = 1$, reduces to the findings of Ref. \citenum{Russell2017DirectMatrices} for locally interacting photonic networks. We will call this the ``phase-beta distribution.''

The key difference between balanced tree and unbalanced tree architectures is in the beta distribution for leaf nodes (connecting directly to the inputs into the device) as shown in Fig. \ref{fig:treestats}(a-c).

Note that relative magnitudes in the mesh also behave like a beta distribution, which is useful for our phase and coupling sensitivity analysis. To see this, consider node $n \leq N - 1$ leading to $N'$ outputs with relative output power $p_n$:
\begin{equation} \label{eqn:waveguidestats}
    \begin{aligned}
        p_n &= \frac{P_{N'}}{P_{N'} + P_{N - N'}} \\
        p_n &\sim \mathrm{Beta}(N', N - N'),
    \end{aligned}
\end{equation}
which shows that \textit{relative} powers in the mesh follow beta distributions as well as the nodes of the mesh. Note that the powers in the mesh are not independent because of the tree structure, and therefore, these relative powers can only be used to evaluate sensitivities of the circuit to individual nodes, but not how the overall circuit responds to a perturbation of all nodes at once. We will address this point specifically later.

For specifying the power going through phase shifts $\eta$ in the device more specifically, we define the following based on previously defined $p_n, s_n$ at the various nodes:
\begin{equation} \label{eqn:phaseshiftpowers}
    \begin{aligned}
        p_{\eta} &= \begin{cases}
            p_n / 2 & \eta = \theta_n\\
            p_n s_n & \eta = \phi_n,
        \end{cases}
    \end{aligned}
\end{equation}
where the factor of $2$ comes from the fact that half the power entering the node goes through the top $\theta$ phase shift.

Interestingly, the variance of $\mathcal{P}_{\mathrm{B}, \theta}(\theta; N_1, N_2)$ scales such that it is near-constant for a given $N$, i.e. balance no longer affects the error model for $\theta$. However, the variance does depend strongly on $N$, and binary trees have a scale invariant property that half the nodes in the network always have $N = 2$, or by our notation are $\mathrm{B}_1^1$ nodes, which explains why balanced trees are so much more error tolerant and broadband as compared to unbalanced trees. More generally, balanced trees have the property that $N / 2^\ell$ nodes in column $\ell$ are $\beta_{2^\ell}^{2^\ell}$ nodes as shown in the labelled waveguides in Fig. \ref{fig:treestats}(c).

The number of possible binary tree architectures for a given $N$ is given by the Catalan number $C_{N - 1} = \frac{(2N)!}{(N + 1)! N!}$. Application of Stirling's approximation suggests that the number of possible vector unit designs scales as roughly $4^N$, and though we only consider the two extremes of these designs, our theory may be applied to any of these designs by considering an arbitrary choice of $N_1$ at each recursive step for Eqs. \ref{eqn:binarytree} and \ref{eqn:selfconfiguringproof}, and then applying the statistics of Eq. \ref{eqn:phasebeta}.



\subsection{Balanced trees are robust to phase error}

We now compare the Hessian and the error scaling in balanced and unbalanced binary tree networks, the core theoretical result of this paper. Our goal is to (1) analytically show that balanced trees are robust to phase error compared to unbalanced trees and (2) perform simulation analysis to numerically verify our scaling arguments given both coupling and phase errors.

Comparing these large class of networks requires incorporating the statistical analysis of flow of power in the network in Fig. \ref{fig:treestats}. Specifically, the Hessian elements for the 8-tree shown in Fig. \ref{fig:hessian} arise from the mean of the beta distributions representing powers in the various waveguide segments as labelled in Fig. \ref{fig:treestats}(c). As a result, we find that the Hessian matrix elements for the balanced tree are much smaller than those of the unbalanced tree. Ultimately, we attain expressions for the power in different segments of the network $p_n$ which is entirely sufficient to determine the necessary Hessian terms and even the distributions of those Hessian terms given random inputs.

Under a simplified assumption of uncorrelated error, it is possible to directly compare the performance of balanced and unbalanced trees by simply tracking the amount of light in various branches of the photonic network. In particular, we can assume that the uncorrelated error is the dominant contributing term in Eq. \ref{eqn:errorcorr} where $\delta \sim \mathcal{N}(0, \sigma^2)$. The overall sensitivity can therefore be thought of as simply the sum of optical powers in various branches of the network assuming a complex random input distribution from which $p_n$ statistics arise. 

Let us now assess the scaling relations according to this now simplified framework for uncorrelated error. We apply Eqs. \ref{eqn:waveguidestats} and \ref{eqn:phaseshiftpowers} to compare expected (average) errors $\mathbb{E}[\varepsilon^2(\bm{\Delta})]$ in balanced networks and unbalanced networks using linearity of expectation:
\begin{equation} \label{eqn:hessianpower}
    \begin{aligned}
        \mathbb{E}[\varepsilon^2(\bm{\Delta})] &= \sum_\eta \mathbb{E}[p_\eta] \\
        \mathbb{E}[\varepsilon^2_{\mathrm{bal}}(\bm{\Delta})]& = \sum_{k = 1}^{\log N} \sum_{k' = 1}^{N / 2^k} \frac{2^k}{N} \sigma^2 = \log N \sigma^2 \\
        \mathbb{E}[\varepsilon^2_{\mathrm{unbal}}(\bm{\Delta})] &= 2 \sum_{n = 1}^N \frac{n}{N} \sigma^2 \propto N \sigma^2 \\
    \end{aligned}
\end{equation}
where $N_\eta$ is the number of inputs spanned by phase shifter $\eta$'s subtree, and $\sigma$ represents uncorrelated phase shift error yielding $\log N$ scaling for balanced trees and $N$ scalings for unbalanced trees. This scaling argument is the key argument of this paper and applies not only for phases but for coupling too as we will see in our numerical analysis. This is not surprising due to the relationship between coupling and phase error uncovered by Ref. \citenum{Bandyopadhyay2021HardwarePhotonics}.

If $N$ is not a power of 2, the structure of a fully balanced architecture is slightly different but the Hessian scaling relations are similar. Additionally, we empirically find that the scaling laws for the phase shifter and for the coupler errors are roughly the same (shown in Fig. \ref{fig:error}), which as shown earlier in Eq. \ref{eqn:nodeerror} can be attributed to the fact that coupling and phase error map to each other \cite{Bandyopadhyay2021HardwarePhotonics}. One notable difference though is the constant phase error for unbalanced architectures in Fig. \ref{fig:error}(d) has an error proportional to $N$ rather than $\sqrt{N}$, and this is explained by the Hessian off-diagonal contributions as shown in Fig. \ref{fig:hessian}(b).

We have now shown how to evaluate the distribution of phase errors in any node within any vector unit implementing a random vector as well as various correlations across elements of the circuit. At a high level, the amount of light present in the ``average device'' (device given average input) is larger for an unbalanced tree as compared to a balanced tree. As previously indicated in Eq. \ref{eqn:powererror}, the power in each waveguide segment of an ``average device'' is also a measure of how sensitive that part of the circuit is and how it effectively couples its error to other parts of the circuit. This is not surprising qualitatively, but it is useful to also quantify the sensitivity of optical circuits based on powers expected in the circuit, and our Hessian formalism accomplishes exactly that. The proportion of inputs that are connected to a waveguide segment gives the statistics of the phase sensitivities of an individual node, and the unbalanced tree tends to ``load'' all of its nodes with the maximum possible inputs. later on, specifically Fig. \ref{fig:hessianstats}, we specifically plot the Hessian statistics for balanced and unbalanced trees and observe they obey the appropriate beta statistics predicted by Eq. \ref{eqn:powererror}.

\section{Mode decomposition networks} \label{sec:apps}

We now construct the unitary device operators $U$ to perform matrix-vector products $\bm{y} = U \bm{x}$ by cascading and/or interleaving vector units in various forms, which establishes the many applications of programmable photonic circuits discussed in this work. There are many useful implementations of $U$, some universal and some not universal (i.e., not all unitary matrices can be programmed).

\subsection{Binary tree cascade}

\begin{figure*}
    \centering
    \includegraphics[width=\linewidth]{cascade.pdf}
    \caption{(a) Universal binary tree cascade consists of a sequence of analyzers that are programmed (or equivalently, mathematically computed) via self-configuration to implement rows of $U$ in order $8, 7, 6 \ldots 2$. (b) The SVD architecture implements any arbitrary unitary operator by connecting universal architectures on either end of a set of MZI attenuators with a free input and output. (c) When the error is a constant factor, the balanced architecture wins for for small $N$, but the margin decreases until it is roughly the same at $M = 3N / 4$. (d) When the error is random, the balanced architecture wins for small $N$, but the margin decreases until it is roughly the same at $M = N$. For (c, d), we compare the component error to the overall error and find that constant error is significantly larger than random error due to higher sensitivity to biased error (e.g. circuit bandwidth) versus random error (e.g. fabrication error). (e) Full singular-value decomposition consists of two universal cascades facing each other wih MZI attenuators in between. (f) Low-rank singular-value decomposition consists of two low-rank cascades facing each other with MZI attenuators connecting the relevant dimensions whose bases can be defined arbitrarily.}
    \label{fig:network}
\end{figure*}


For this section, we define a new ``multimode" error function that generalizes the ``single mode" error of Eq. \ref{eqn:vectorerror} to $M \leq N$ orthogonal basis vectors here denoted as $\bm{x}_m$, i.e. not necessarily $N$ ``full rank'' basis vectors. Specifically, we define $\epsilon_{N, M}$ as follows:
\begin{equation}
    \epsilon_{N, M}^2 = \sum_{k = 1}^M\frac{\|\bm{x}_m - \hat{\bm{x}}_m\|^2}{M} = \frac{\mathcal{R}(2 - 2\mathrm{tr}(U_M^\dagger \hat{U}_M))}{M},
\end{equation}
assuming $U_M$ is a $N \times M$ matrix with $M$ orthogonal basis vectors (so $U_M^\dagger U_M$ is $M \times M$). If $M = N$, $U$ is a square matrix and is unitary. The case $M = 1$ degenerates to our previous definition for single vectors in Eq. \ref{eqn:vectorerror}.

In a ``binary tree cascade,'' we use $M$ vector units to construct a set of $M \leq N$ normalized and mutually orthogonal basis vectors that form $U_M$ (which becomes square and unitary when $M = N$). Self-configuring networks in the default orientation defined in Eq. \ref{eqn:binarytree} can be cascaded to form universal unitary architectures, namely architectures that can be used to define any unitary matrix. This follows from the proofs in Refs. \citenum{Reck1994ExperimentalOperator, Miller2013Self-configuringInvited, Miller2013Self-aligningCoupler}, but we will describe a simplified proof here. A \textit{cascade} of self-configuring networks can be defined as:
\begin{equation} \label{eqn:cascade}
    \begin{aligned}
    U_{N, M} &\equiv \prod_{n = 1}^M R_n^{[N]}(\overline{\bm{u}}_n), \\
    R_n^{[N]} &= \begin{cases}
         R_N & n = N \\
         \begin{bmatrix}
                R_{n} & O \\
                O & I_{N - n}
        \end{bmatrix} & \mathrm{otherwise}
        \end{cases}
    \end{aligned}
\end{equation}
where $\bm{u}_n$ is the $n$th column of $U$, $\overline{\bm{u}}_n$ is the $n$th column of $U$ after passing through the first $n$ layers, and $I_{N - n}$ is an identity matrix of size $N - n$. Finally, $R_n^{[N]}$ means, in words, a self-configuring layer matrix representation over the first $n$ inputs of an $N$-waveguide rail system. More explicitly we define $\overline{\bm{u}}_m$ using the recurrence relation:

\begin{equation}
    \overline{\bm{u}}_m = \left(\prod_{n = 1}^{m - 1} R_n^{[N]}(\overline{\bm{u}}_n) \right) \bm{u}_m  \\
\end{equation}

Due to orthogonality, we always have that the last $m$ elements of $\overline{\bm{u}}_m$ are zero, since $\bm{u}_m$, a row in $U_N$ is orthogonal to all other rows in $U$, including the preceding $\bm{u}_1, \bm{u}_2, \ldots \bm{u}_{m - 1}$.

Now, we rewrite Eq. \ref{eqn:cascade} in terms of the actual parameters of the physical system that must be programmed:
\begin{equation}
    U_{N, M}(\mathbf{s}, \bm{\phi}) = \prod_{n = 1}^K D_n^{[N]}(\mathbf{s}_n, \bm{\phi}_n), \\
\end{equation}
where we have applied the self configuration process in Eqs. \ref{eqn:phisc}, \ref{eqn:ssc}, and \ref{eqn:selfconfiguringproof} for each layer in order from left-to-right starting from $R_N \to D_N$ to $R_2 \to D_2$.

While there exist architectures for any binary tree, some vector units can be more compactly cascaded than others. For instance, the unbalanced tree (diagonal line), can be cascaded to form a triangular architecture that is $2N - 3$ nodes deep. However, a balanced binary tree cannot be packed compactly and requires up to $N \log N$ photonic layers to implement. Therefore, for larger $N$, balanced binary tree cascades should ideally be used when $M \ll N$ is sufficient to solve some problem, with the key benefit being that $N$ can now be much larger than what would typically be used for a rectangular or triangular locally interacting network. In such a case, we would have a total of just $M \log N$ layers; an interesting case that warrants further investigation would be a binary tree with $M = N / \log N$ units, which would be $N$ layers deep, equal to the optical depth of a rectangular universal network with a tradeoff of fewer degrees of freedom \cite{Clements2016AnInterferometers}.

\subsection{Singular value decomposition networks}

A singular value decomposition (SVD) network, as first proposed in Ref. \cite{Miller2013Self-configuringInvited} and shown in Fig. \ref{fig:network}(e), is capable of performing any arbitrary complex linear operation (i.e., not just unitary), and is actually a specific case of the cosine-sine (CS) photonic mesh (discussed in the Appendix). The SVD architecture provides the necessary $2N^2$ degrees of freedom by decomposing a matrix in the form $A = U \Sigma V^\dagger$ where $A$ is $N \times M$ matrix, $U, V$ are $N \times N$ and $M \times M$ unitary matrices and $\Sigma$ is $N \times M$ matrix with singular values along the diagonal, i.e. $\Sigma_{nn} = \sigma_n \in \mathbb{C}$ and $\Sigma_{nm} = 0, n \neq m$. The $\Sigma$ matrix is represented as an array of $\min (M, N)$ coupling matrices $X(\theta)$ placed between the universal architectures implementing $U$ and $V^\dagger$, where each element of the coupling matrix act as an effective attenuator (tunable loss element) (with an additional phase shifter if require to implement a complex coupling element).

This feature of these networks is particularly important as we consider using wide or low-rank rather than square matrix multiplication techniques for matrix acceleration architectures. In particular, if $ M \ll N$, then a cascade of $M$ balanced tree architectures for inputs of size $N$ and a small unitary $M \times M$ architecture may be a prudent strategy for scaling up such architectures and reducing the optical depth (number of devices light has to pass) for a lower-loss and lower-systematic error device. This architecture is shown in Fig. \ref{fig:network}(f).

In signal processing applications, there are many cases where a small number of principal components are necessary to characterize some system, which can benefit from the architecture of Fig. \ref{fig:network}(f). This is certainly the case in mode conversion and telecommunications applications \cite{Miller2013EstablishingAutomatically, Annoni2017UnscramblingModes}. However, principal components analysis (PCA) is also a commonly employed tactic for reducing the dimensionality of a problem, first popularized in a 1991 paper by Turk and Pentland \cite{Turk1991EigenfacesRecognition}. Our scalability arguments suggest that future work in assessing the value of low-rank photonic architectures in principal components analysis could very possibly lead to a new paradigm for photonic network-based computing and analog signal processing that is significantly more robust and scalable to large numbers of input modes $N$.

As a final note, many of the error tolerance scaling properties for low-rank matrices in Figs. \ref{fig:error} and \ref{fig:network} apply to SVD architectures. Because random binary tree cascades also encode random unitary matrices, it is straightforward to assign $\mathbf{s}, \bm{\phi}$ to the appropriate beta distribution that depends on the number of inputs in the  cascade meshes corresponding to $V^\dagger, U$ respectively.

\section{Discussion and Conclusion}

A major theme of this work is that wide, balanced architectures are more robust than deep, unbalanced architectures by a factor of roughly $N / \log N$. This should be evident from the fact that light must pass through more components in deeper architectures, resulting in a larger sensitivity. This general trend can be explained by the theory of error tolerance provided in Sec. \ref{sec:error} which explores the effects of both uncorrelated and correlated errors on the overall photonic device performance. Our analysis in Sec. \ref{sec:error} indicates that much of the error is concentrated in waveguides that propagate more optical power and in Sec. \ref{sec:results}, it is revealed that in a balanced tree, those are the waveguides routed closer to the root of the tree. Defining shorter (i.e., less phase sensitive) and more robust components near the root of the tree is therefore paramount in scaling photonic technologies.

Equipped with this background, it is now important to discuss some specific applications where our theory may be considered, such as analog computing and sensing when the rank (number of supported optical modes) is low, i.e. $M \ll N$. The focus in this paper has been in matrix rank (Fig. \ref{fig:network}) and error tolerance (Figs. \ref{fig:error} and \ref{fig:network}), which are considered along with other application-specific criteria such as footprint, power consumption, loss, and speed.

In transceiver applications such as sensing, LIDAR, and telecommunications \cite{Bogaerts2020ProgrammableCircuits, Annoni2017UnscramblingModes}, the main factor is the speed of phase modulators as well as the scalability in the mesh needed for optical phased array communications. While speed can be a crucial factor in wireless data transfer nodes, high resolution optical phased arrays require large numbers of emitters (say, a million outputs or a $1024 \times 1024$ array), and our error analysis in this paper is highly relevant to reaching devices of such scale in photonics. In particular, as proven throughout this paper, our analysis can help elucidate the optical sensitivities of various parts of a optical phased array circuit. Such circuits are typically a balanced tree design (which may follow an H-tree fractal) with errors (sensitivities) in various waveguide segments in the passive splitter circuit leading up to the all-important final output phase shift array. Thankfully we have shown that the overall sensitivity of a balanced photonic network scales with $\log N$, making the prospects for scaling up to a million phase-controlled emitters more feasible given $\log N = 20$ layers. 

Separately, we have mentioned in the Introduction how photonic computing (e.g., matrix-vector multiplication) tasks can suffer from low error tolerance at sufficiently high circuit sizes, e.g. for machine learning \cite{Shen2017DeepCircuits, Pai2022ExperimentallyNetworks} and cryptography or blockchain \cite{Pai2022ExperimentalCryptocurrency}. Computing tasks are evaluated in terms of OPS (operations per second), which are limited by detection limits like integration time for sufficiently low signal-to-noise ratio and modulation switching speed limits for setting up inputs into the feedforward mesh. It is predicted that we might be able to achieve petaops level efficiencies using photonic mesh circuits of sufficient circuit size $N > 64$ \cite{Pai2022ExperimentalCryptocurrency, Nahmias2020PhotonicNetworks}. 

An alternative approach for photonic computing inspired by our theory is to embrace the ``wide'' over ``deep'' architectures. If we engineer an architecture to minimize the number of layers light has to propagate through, we can reduce both the overall photonic loss as well as the overall accumulated error. This concept is introduced as part of our ``splay architecture'' framework in the Appendix, which has a larger footprint but is significantly more error tolerant and possibly a less lossy scheme. This same theme is obeyed by Hadamard and FFT-like photonic mesh networks, also discussed in the Appendix, are nested binary trees that implement a subset of unitary space but use nonlocal connections to efficiently couple large numbers of modes ($N$) with fewer degrees of freedom ($N \log N$ vs $N^2$). Such efficient unitary representations have already been suggested for machine learning applications and also implement FFT and permutation operations, both of which have wide-ranging applications.

Another key advantage of wide, balanced architectures is speed. As discussed in the Appendix, fast calibration (self-configuration) is also important for any applications (including computing and possibly sensing) where interrupting device execution to re-calibrate the system is important. Self-configuration for instance is much faster in balanced trees and balanced tree cascades versus unbalanced trees and unbalanced tree cascades. In conjunction with higher tolerance to error, this means balanced trees need to be calibrated less often, and when they do need to be calibrated, the calibration is likely faster as well.

In conclusion, we find that low-rank ``wide'' photonic computing offers high error tolerance and bandwidth as well as more efficient device operation and calibration which can serve a large range of signal processing applications. Our introduction of binary tree cascades and splay networks (and analysis using our theoretical framework of binary tree architectures) is a key step to realizing higher-scale photonic circuits.

\section*{Acknowledgements}

We would like to thank Nathnael Abebe, Annie Kroo, Dirk Englund, Ryan Hamerly, Srikrishna Vadlamani, and Saumil Bandyopadhyay for helpful discussions. Additionally, we would also like to acknowledge funding from Air Force Office of Scientific Research (AFOSR) grants FA9550-17-1-0002 in collaboration with UT Austin and FA9550-18-1-0186 through which we share a close collaboration with UC Davis under Dr. Ben Yoo.

\section*{Data and software}
All software for running the simulations and Hessian calculations are available via Simphox \cite{Pai2022Simphox:Library}.

