Asymptotically fault-tolerant programmable photonics

Component errors limit the scaling of programmable coherent photonic circuits. These errors arise because the standard tunable photonic coupler—the Mach-Zehnder interferometer (MZI)—cannot be perfectly programmed to the cross state. Here, we introduce two modified circuit architectures that overcome this limitation: (1) a 3-splitter MZI mesh for generic errors, and (2) a broadband MZI+Crossing design for correlated errors. Because these designs allow for perfect realization of the cross state, the matrix fidelity no longer degrades with increased mesh size, allowing scaling to arbitrarily large meshes. The proposed architectures support progressive self-configuration, are more compact than previous MZI-doubling schemes, and do not require additional phase shifters. This removes a key limitation to the development of very-large-scale programmable photonic circuits.

Large-scale programmable photonic circuits are opening up radical new possibilities for optics.Of central importance in many devices is the universal multiport interferometer, which functions as an N × N programmable linear circuit (Fig. 1(a-b)).This device, usually constructed from a dense mesh of Mach-Zehnder interferometers (MZIs) [1,2], is widely employed in applications ranging from spatially multiplexed optical communications to machine learning and quantum computing [3][4][5][6][7].Sadly, component errors (Fig. 1(c)) are a critical factor limiting the size of such circuits.Since the circuit depth of MZI meshes scales as O(N ), the effect of errors grows with mesh size, meaning that, in practice, even modestly sized circuits cannot be programmed to high accuracy.Motivated by this challenge, a large body of recent work has focused on "correcting" hardware errors by global optimization [8][9][10], self-configuration [11][12][13][14][15][16][17][18], or local correction [19,20].For conventional MZI meshes, correction reduces errors by a quadratic factor [16,19]; however, the effect of errors still grows with mesh size and poses a fundamental limit to the scaling of these circuits.
To overcome this limit, various alternative mesh architectures have been proposed.Non-compact structures such as binary trees avoid the extreme splitting-ratio requirements [21,22], but suffer from large chip area and the need for many crossings.A complementary approach is to stick to conventional geometries [1,2], but insert redundant MZIs to realize the full range of splitting ratios even in imperfect hardware [23][24][25].This solves the scaling problem, but at the cost of a 1.5-2× increase in the number of splitters and phase shifters.The resulting effects on chip area (particularly on emerging high-speed platforms where phase shifters have a large footprint [26,27]), waveguide length (which affects insertion loss and latency [28]), and electronic complexity (number of pads, traces, DACs / drivers, etc.) make this option unappealing.
In this paper, we propose two mesh architectures that achieve the same perfect scaling without significant added complexity: a 3-splitter MZI that corrects all hardware errors (Fig. 1(d)) and an MZI+Crossing design that only corrects correlated errors, but has the added advantage of broader bandwidth (Fig. 1(e)).These designs take up significantly less chip area than the "perfect" redundant MZIs [23,24], and do not require additional phase shifters.Moreover, the proposed architectures support progressive self-configuration [16,17], allowing for error correction even when the hardware errors are unknown.This work will enable the development of freely scalable, broadband, and compact linear photonic circuits.This paper is structured as follows: first we introduce the formalism of error correction in MZI meshes, focusing on the self-configuration approach.Splitting ratios are visualized as points on the Riemann sphere, where hardware imperfections lead to forbidden regions around the poles (bar-and cross-state), where the probability density is at a maximum.To avoid this unfortunate coincidence, our architectures "rotate" the Riemann sphere to move the forbidden regions away from this peak, so that a larger fraction of MZIs are perfectly realized.Based on this concept, we introduce the 3-splitter MZI, which can correct arbitrary errors by rotating the forbidden regions to the equator.Using a benchmark optical neural network, we show that this modified MZI mesh is > 3× more robust to hardware errors, enabling accurate inference in a regime where standard interferometric circuits struggle.Finally, we introduce the MZI+Crossing, which flips the poles of the Riemann sphere.While this design is only robust against correlated errors, it has the added advantage of broader intrinsic bandwidth.For both architectures, we compare the matrix fidelity to the standard MZI to demonstrate the scaling advantage of both schemes.

Error Correction Formalism
To correctly configure an MZI mesh in the presence of errors, one uses a nulling method based on physical measurements [16,17].Fig. 2(a) illustrates the case of the triangular mesh [1], where the procedure is more straightforward.The transfer matrix for this system is a product of a phase screen D and a sequence of 2 × 2 unitaries W : where T mn is the n th MZI of the m th rising diagonal.We configure the mesh by building up matrix W in a sequence of steps designed to diagonalize a target matrix X = U W † .In each step, we add one crossing to W , performing the update W → T mn W , which right-multiplies the target matrix X → XT † mn (Fig. 2(b)).The phase shifts (θ, φ) are chosen to zero a particular matrix element v → 0 (green in figure), satisfying the equation (indices m, n suppressed for notational simplicity): This is illustrated in Fig. 2(c).Nulling physically corresponds to injecting w * j (the j th column of W † ) and zeroing the power at the i th output [17].If all nulling steps are performed exactly, the mesh will perfectly realize the target matrix U (see Methods and Supp.Sec.S1 for details).
Mathematically, nulling corresponds to matching the complex splitting ratio s ≡ T 11 /T 12 = −(T 22 /T 21 ) * to a target value ŝ ≡ u/v.This is not always possible, as the range of splitting ratios tan |α + β| ≤ |s| ≤ cot |α − β| is constrained by hardware imperfections, namely the splitting-angle errors α, β for the 50:50 couplers in a real MZI (Fig. 1(c)).These imperfections lead to forbidden regions (Fig. 2(d)) for small and large s, where nulling cannot be achieved perfectly.It is also instructive to view this chart on the Riemann sphere, which shows that these forbidden regions are centered around the poles (Fig. 2(b)), highlighting the well-known fact that imperfect MZIs generally have finite extinction ratio and cannot realize a perfect cross (s = 0) or bar (s = ∞) state.If in a given nulling step ŝ falls within the forbidden region, nulling is imperfect, and an off-diagonal residual prevents perfect diagonalization of the matrix, leading to an "uncorrectable" error.This residual is proportional to d(s, ŝ), the Euclidean distance on the Riemann sphere between the target ratio and the closest realizable s.The overall error is the quadrature sum of all such residuals.
For linear photonic circuits, two important fidelity figures of merit are (1) the coverage C, i.e. the probability that a matrix is realized exactly, and (2) the normalized matrix error E = ∆U rms / √ N , which is approximately equal to the average relative error for a given matrix element.
C and E depend on the error model and the distribution of target matrices.Here, consistent with prior work [16,17,19,29], we sample target matrices randomly over the Haar measure [30,31] and consider an uncorrelated Gaussian error model α rms = β rms = σ.Analytic expressions for E and C are derived in the Methods, which we summarize here.If a mesh is straightforwardly programmed without taking any account of the imperfections ("uncorrected" error), the normalized error is E 0 = √ 2N σ [16,19].The coverage C = e −N 3 σ 2 /3 (Eq.16) decreases sufficiently fast that even moderately sized meshes have vanishingly small coverage, and error correction is generally imperfect.In this case, the residual "corrected" error E c = (2/3)N σ 2 (Eq.19) is the more relevant metric.Since E c ∝ (E 0 ) 2 , self-configuration correction affords a quadratic suppression of errors, which is a significant advantage when errors are below a threshold.However, for sufficiently large meshes N 1/σ 2 , error correction will be ineffective and the mesh cannot realize most matrices at high fidelity.Thus, even with error correction, hardware imperfections set a fundamental scaling limit for standard MZI meshes.

Asymptotically Perfect Photonic Circuits
The main challenge limiting error correction here is that the forbidden regions overlap with the peak of the probability distribution, which clusters tightly around the cross state s = 0 (Fig. 2(e)) [29].This clustering happens because light must propagate all the way down a mesh's diagonals to realize generic unitaries; the forbidden regions disrupt this ballistic transport leading to clipping of off-diagonal matrix elements [10].Adding redundant components (MZI doubling) solves this problem by eliminating the forbidden regions altogether [23,24], but at the cost of added optical and electrical complexity.Here, we take the alternative approach of displacing the forbidden regions away from the cross state.This can be performed by placing a third splitter at the input of the MZI, as shown in Fig. 3(a).The extra splitter performs a Möbius transformation s → (s + i tan η)/(1 + is tan η), which for a 50:50 splitting ratio (η = π/4) maps the bar and cross states to s = ±i (Fig. 3(b)).This can be visualized as a 90 o rotation on the Riemann sphere, which pushes the forbidden regions to the equator, while the probability density is still concentrated at the poles (small errors γ in the third splitter perturb this rotation angle slightly, but this does not change the structure of the forbidden regions and has little effect on the error correction).
This "3-splitter MZI" (3-MZI) can realize the full range of (absolute value) splitting ratios |s| ∈ [0, ∞), and can thus function as a high-contrast optical switch [24,32].However, the presence of forbidden regions means that the relative phase of this splitter cannot be fully controlled; which means that errors can still occur when programming the mesh (unlike the "perfect" MZIs of Refs.[23][24][25], which cure this defect with redundant phase shifters).However, from the distributions in Fig. 2(e), for large meshes ŝ will fall into the 3-MZI's forbidden regions only rarely.The normalized matrix error, calculated in the Methods (Eq.22), takes the following form: In Fig. 3(c-d), we numerically simulate self-configuration on imperfect meshes using the Meshes package (see Methods and supplemental code); the realized E c shows good agreement with Eq. (3).For most mesh sizes, the E c is 1-2 orders of magnitude smaller for the 3-MZI design.
Remarkably, the error actually decreases with increasing mesh size, scaling as E c ∝ log(N )/N .In the asymptotic limit N → ∞, matrices can be programmed perfectly.
This non-intuitive effect arises from the fact that, under the Haar measure, only a small fraction of MZIs have significant probability density near s = ±i, where the forbidden regions are centered [29].This probability decreases exponentially with the distance from the triangle's base (see Methods for details).Therefore, although the mesh has N (N − 1)/2 MZIs, only O(N ) contribute significantly to the matrix error under self-configuration.
A naïve estimate assuming uncorrelated errors would give ∆U ∝ √ N σ 2 , which would lead to a constant E c .However, during the self-configuration process, subsequent MZIs can partially correct for errors in earlier MZIs that cannot be properly configured; the end result is to reduce the overall error of each MZI by a factor proportional to log(N )/N (see Methods), yielding the result Eq. (3).
Another advantage of the 3-splitter MZI is that the threshold for perfect error correction is higher.One obtains this threshold is found by computing the coverage C = e −16N σ 2 (see Methods, Eq. ( 20)).This is much larger than the coverage of the regular MZI mesh, and the threshold scales as σ th ∝ N −1/2 , as opposed to the N −3/2 scaling observed for the standard mesh.Consequently, perfect error correction is available under a much wider range of conditions, as shown in Fig. 3(e).

Error-Resilient Optical Neural Networks
To highlight the significance of this error reduction, consider as a concrete example deep neural network (DNN) inference on coherent optical hardware.A DNN is a sequence of layers, consisting of linear synaptic connections and nonlinear neuron activations.An emerging application of photonics seeks to use optical interference to accelerate this process, encoding neuron activations in coherent optical amplitudes, while a programmable MZI mesh implements the synaptic weights and activations are performed with an all-optical or electo-optic nonlinearity [5].Scaling remains the major challenge to constructing practical optical neural networks, as large mesh sizes (N > 100) are required to achieve a significant advantages over electronic hardware, and such large meshes are especially susceptible to fabrication errors.A recent numerical study showed that even with state-of-the-art process tolerances, hardware errors can significantly degrade DNN inference accuracy [34], a difficulty that has spurred investigations into alternatives to the MZI mesh, which all have their own limitations [35][36][37][38].
Fig. 4(a) depicts a benchmark neural network.Here, 28 × 28 images from the MNIST digit dataset [39] are preprocessed by a Fourier transform and cropped to a window of size √ N × √ N , which forms the input to a twolayer unitary DNN.The DNN can be implemented optically with rectangular MZI meshes for synaptic weighting [2] and electro-optic nonlinearities for the activation (see Refs. [17,33] for details).Models with inner-layer sizes N = 64 and N = 256 are pre-trained using the Neurophox package [40], and inference accuracy is sub-sequently simulated on imperfect meshes with Gaussian splitter errors to calculate the classification accuracy.This accuracy is plotted in Fig. 4(b) for three cases: straightforwardly programming an MZI mesh without error correction, with error correction, and with the modified 3-MZI architecture.Even for small device errors σ = 1-2%, which is considered state-of-the-art for directional couplers in highly controlled fabrication processes [41], hardware errors significantly degrade the model's inference accuracy relative to its canonical value (σ = 0).For small σ, this is recovered using error correction [17,19].However, many broadband coupler designs [42][43][44][45][46][47] trade bandwidth for fabrication sensitivity and are in practice very sensitive to process variations, meaning larger splitter errors σ 5% are common.In this moderate-error regime, error correction alone is not sufficient and the network shows reduced accuracy, a problem that becomes more pronounced as the size N increases.Moving to the 3-MZI architecture overcomes this limitation, enabling effectively error-free inference (relative to the canonical model) even out to very large splitter errors σ ≈ 10-15%, far beyond what is likely to be encountered in practice.

Broadband Mesh for Correlated Errors
For generic, uncorrelated component errors the 3-splitter MZI is well-suited.However, since the correlation lengths of process variations tend to be larger than a single MZI [48], errors are correlated in practice.This is especially true for broadband couplers based on multimode interference (MMI) [42,43], subwavelength gratings [44,45], and asymmetric designs [46,47], all of which are highly dependent on the device geometry, which can vary slightly from run to run.Moreover, even with perfect 50:50 couplers, the splitting ratios are still wavelength-dependent.Oper- ating the mesh away from its design wavelength leads to correlated device errors, so sensitivity to these errors is closely tied to the operational bandwidth of the device.
Consider the case of a constant offset µ for all splitting ratios: α = β = µ.In a standard MZI, the bar-state forbidden region (around s = ∞) disappears since |α − β| = 0, while the cross-state region (around s = 0, the peak of the probability distribution) remains in place (Fig. 2).This is consistent with the common observation that the extinction ratio in an MZI is much higher in the cross port than in the bar port.The optimal error reduction strategy, illustrated in Fig. 5(a), was previously proposed in the context of broadband optical switching: place a waveguide crossing before the MZI [49].The added crossing performs the Möbius transformation s → 1/s, rotating the Riemann sphere by 180 o to move the forbidden region to the minimum of the probability distribution (Fig. 5(b-c)).
As before, we can calculate the coverage and matrix error of this "MZI+Crossing" (MZI+X) mesh by performing the nulling procedure on target unitaries, obtaining C from the probabilities that splitting ratios fall within the forbidden regions, and E c from the residuals arising from imperfect diagonalization.In this case, there is only one forbidden region, centered at s = ∞.The calculation is worked out in the Methods.For the normalized error, we find (Eq.27): This is plotted in Fig. 6.Like the 3-MZI design, this metric scales as E c ∝ log(N )/N µ 2 , in contrast to the trend E c = (4/3 3/2 )N µ 2 calculated for the standard MZI under correlated errors.The coverage also increases (Eq.( 25)), so that the threshold for perfect correction likewise scales as µ th ∝ N −1/2 , as opposed to µ th ∝ N −3/2 for the standard mesh.
Ultimately, the scalability of the MZI+X architecture is limited by differential errors |α − β| that arise from local fluctuations in waveguide dimensions.The effect of such errors is analyzed in Supp.Sec.S2.For typical photonic process variations, |α − β| µ and differential errors are insignificant for mesh sizes up to at least N = 512.
As an added bonus, the MZI+X design also reduces the effect of errors in the absence of correction.To see how, we can make an analogy to Bloch-sphere rotations.The transfer matrix of a standard MZI is (up to a phase factor) the product of four rotations: where R k (η) = e iσ k η is a Pauli rotation and σ k is a Pauli matrix.For the cross state (θ = 0), the errors µ add up constructively, while for the bar state (θ = π), they cancel out (the latter is a simple example of dynamical decoupling of spins using a pulse sequence).Most crossings in large meshes are close to the cross state, which leads to constructive addition of the errors in the standard MZI mesh.However, for the MZI+X, the input ports of each MZI are exchanged, so the physical MZIs are close to the bar state where the errors cancel out.The resulting uncorrected matrix error is (see Methods): Correlated errors (both corrected and uncorrected) are important because they are tightly connected to the operational bandwidth of the mesh, a critical design parameter for machine learning schemes that require broadband operation, e.g. for parallel processing on wavelength-multiplexed data [50][51][52][53]  wavelength-dependent splitter error, which can usually be expanded to first order µ ≈ (dµ/dλ)∆λ.Two important wavelength-dependent figures of merit are (1) the tuning range, which refers to the range of λ over which the mesh can be programmed to a given accuracy, Fig. 7(a, c), and (2) the bandwidth, which is related to the number of wavelength channels that can be (simultaneously) processed by the mesh, Fig. 7(b, d).The tuning range is limited by the corrected error E c , while the bandwidth is limited by the uncorrected error E 0 , since a mesh cannot simultaneously error-correct at two different wavelengths.Since the MZI+X design reduces both E 0 and E c , it leads to enhancements in both the bandwidth and tuning range.The enhancement factors scale as and are listed for several mesh sizes in Table 1 (see Methods for details).As Fig. 7(c-d) illustrates, the MZI+X architecture enjoys a significantly larger tuning range, in addition to modestly greater bandwidth.
Real crossings have a small amount of nonzero crosstalk, quantified by the S-matrix element S 21 ; scattering into the forward-facing port leads to a perturbation R x ( π 2 ) → R x ( π 2 + γ) in the transfer matrix, where γ = 10 −S21[dB]/20 .This does not degrade the effectiveness of self-configuration, since the additional scattering angle merely rotates the Riemann sphere Fig. 5(c) by an additional angle γ 1, and the forbidden region is still far from s = 0. In-plane crossings in silicon can achieve sub-40 dB crosstalk suppression (γ < 0.01) with insertion losses well below 0.1 dB [54][55][56][57][58]. Unlike directional couplers, crossings are inherently broadband; the insertion loss and crosstalk depend only very weakly on λ, so any crossing imperfections can be treated as (correctable) wavelength-independent errors that do not affect the bandwidth enhancements of the MZI+Crossing scheme.In addition to the forward-scattered light, a 90 o crossing will scatter light into the backward-facing port.Back-reflected light can be subsequently reflected in other crossings, leading to a spurious signal that interferes with the forward-propagating light.Provided that the phases of reflected beams are random, these add in quadrature: with amplitude γ 2 and O(N 2 ) scattering paths, we expect this to induce an O(N γ 2 ) error, which may be uncorrectable and set a limit on scaling.However, if this effect is small, gradient-based methods or iterative selfconfiguration may enable correction of these errors.

Discussion
As photonic circuits grow larger, error tolerance becomes increasingly important.Many techniques exist to manage hardware errors, but all involve a tradeoff between accuracy and complexity.At opposite poles lie "zerochange" error correction, which has limited scalability [16,17,19,59], and "perfect" photonic circuits, which require a larger number of photonic and electronic components [23,24].This paper has introduced two designs for programmable circuits that strike a tradeoff between these extremes, as shown in Fig. 8 and Table 2, achieving performance that is almost as good as the perfect designs, but with less added complexity (see Supp.Sec.S3 for details).
The main insight from this paper is that, by adding a single passive component (either a splitter or a waveguide crossing) to the MZI, we can recover behavior that is asymptotically perfect-that is, the average normalized matrix error decreases with size.Our design choices are motivated by the elegant theory of self-configuration by matrix diagonalization [17], where splitting ratios are set to successively zero the off-diagonal elements of the target unitary.By visualizing the MZI state on the Riemann sphere, we can intuitively understand the increased error robustness of our designs in terms of "rotating" the forbidden regions away from the peak probability den-  [63], (c) 3-splitter (3-MZI) [32], (d) Portexchanged (MZI+X) [49], (e) Suzuki [24], (f) Miller [23].

Complexity
Features sity.This leads to a several-orders-of-magnitude reduction in post-correction errors compared to the standard MZI mesh.The ability to achieve near-perfect and freely scalable MZI meshes with less complexity than the MZIdoubled designs [23,24] (especially with respect to the number of active components and pads) removes a major obstacle to the realization of very-large-scale photonic circuits.
An interesting direction for future work is to explore to what extent multiport interferometers can be made robust to imperfections in the absence of error correction.For example, previous studies of 3-MZI splitters have noted a wavelength-independent coupling ratio for certain parameter choices [32].Likewise, the near-cancellation of correlated errors in the MZI+Crossing architecture explains the O( N/ log N ) reduction in the uncorrected error, and corresponding increase in bandwidth.Further design modifications based on the theory of composite pulse sequences [60][61][62] may allow this imperfect cancellation to be made exact, further improving the bandwidth (and multiplexing capabilities) of linear photonics.

Methods
Unitaries and the Riemann Sphere A generic 2 × 2 complex-valued matrix has eight degrees of freedom, and a 2 × 2 unitary has four.However, the space of 2×2 unitaries can be divided into equivalence classes based on the splitting ratio s = T11/T12.Specifically, any two unitaries are equivalent up to output phases, i.e.T = diag(e iψ 1 , e iψ 2 ) T , if and only if the splitting ratios are the same, s = ŝ.As a complex number, s can be visualized on the Riemann sphere (Fig. 2(d)), where the mapping is performed by the stereographic projection s = (x + iy)/(1 + z) (which inverts to Ordinarily, the distance between matrices is defined as the Frobenius (L2) norm ∆U = ( mn |∆Umn| 2 ) 1/2 .However, since output phases are corrected in subsequent steps, the most relevant distance metric for a 2×2 block is the Frobenius norm modulo these phase shifts, where d(s, ŝ) = 2|s − ŝ|/ (|s| 2 + 1)(|ŝ| 2 + 1) is the Euclidean distance between two points on the Riemann sphere.

Coverage and Matrix Error Derivation
The nulling method relies on successive zeroing of off-diagonal elements to diagonalize the matrix X (initialized to U ).Each nulling step zeros a single element, increasing the size of the zeroed-out off-diagonal region.Nulling steps are performed in a particular order to ensure that zeroed-out elements remain zero after all subsequent steps [1,2,17].In a given step, if nulling cannot be achieved perfectly, the "zeroed-out" region of matrix X is left with a residual of magnitude: where ŝ is the target splitting ratio, s is the closest physically realizable value, and d(s, ŝ) is the Euclidean distance on the Riemann sphere, the same metric used in Eq. ( 8).
The coverage and matrix error depend on (1) the distribution P (s) of target splitting ratios, a function of the distribution of target unitaries, and (2) the locations and sizes of the forbidden regions, a function of the specific mesh implementation (MZI, 3-MZI, MZI+X).For the Haar measure, P (s) depends on an MZI's location in the mesh; for a given Tmn it takes the following form [29]: Here, the density is defined with respect to the area measure on the Riemann sphere so that Pmn(s)dµ(s) = 1.Note that, under Eq.( 10), Pmn is uniform for the lowest row of crossings, and becomes increasingly concentrated as one approaches the triangle's apex; as a result, the overall distribution is strongly biased towards the cross state for large meshes, as shown in Fig. 2(e) (the same distribution also holds for the rectangular mesh, up to a reordering of the MZIs).
The forbidden regions F± are centered at opposite poles of the Riemann sphere and have radii R± = 2|α ± β|.In the case of small hardware errors, where P (s) ≈ P (s±) inside each F±, the probability that ŝ falls inside the region is given by πR 2 ± P (s±).The coverage C is the probability that every ŝ avoids the forbidden regions, and is well approximated by The normalized matrix error Ec = ∆U rms / √ N is approximately the quadrature sum of the residuals accumulated during nulling: Here, . . .refers to the ensemble average over both Haardistributed target unitaries U [30,31] and the distribution of hardware errors α, β.We calculate the mean residual r 2 by averaging Eq. ( 9) over the distribution P (s).This is simplified in the case of small hardware errors, because the forbidden region is correspondingly small and where we can assume P (s) is approximately constant: A detailed description of the nulling algorithm, including a comparison to the local method [19] and global optimization [8][9][10] (which has a much longer convergence time), is presented in Supp.Sec.S1.
For the MZI mesh, the coverage expression Eq. ( 13) is dominated by the s = 0 term, where Pmn(0) = n/4π.Considering only this term, we calculate: where we have replaced the discrete sum by an integral mn (. ..) which is valid in the limit of large N .Likewise, the top forbidden region dominates the matrix error, so we evaluate Eq. ( 15) including only the first term in the sum: Converting the sum to an integral and substituting R 4 + , we find: Now we redo the calculation for the 3-MZI.In this case, the forbidden regions are located at s± = ±i and contribute equally to the problem.Following Eq. ( 13), the coverage is given by: Applying Eq. ( 15), the mean residual left by crossing Tmn is: The factors of two in Eqs.(20)(21) arise because both forbidden regions contribute equally.This r 2 mn is not slowly-varying with (m, n), so we cannot convert the sums to integrals.We first perform the summation over n, which converges rapidly due to the 1/2 n+1 factor (approximating the upper bound to infinity because of the rapid convergence), followed by summation over m.We find the normalized error: where γe ≈ 0.5772 is the Euler-Mascheroni constant.

Correlated Errors: MZI & MZI+X
Under a correlated error model, α = β = µ.In this case, there is only one forbidden region, which for the MZI is centered at s+ = 0, with R+ = 4µ.The coverage and matrix error for the standard MZI can then be calculated from Eqs. (16,19) with the appropriate substitutions for R 2 + , R 4 + : Now consider the MZI+X.The additional crossing rotates the forbidden region to s+ → ∞.Only the MZIs in the bottom row of the triangle (n = 1) contribute to the sums in Eqs.(13)(14), because the probability distribution Eq.(10) vanishes at s = ∞ for the upper rows.
As before, we use the residual formula Eq. ( 15) to calculate the matrix error.In this case, there is only one forbidden region, centered at s+ = ∞, with R+ = 4µ.Only the MZIs in the bottom row contribute to the sum, because the probability distribution Eq.(10) vanishes at s = ∞ for the upper rows.The coverage is: With the mean residual given by and r 2 mn = 0 for n > 1, the matrix error evaluates to: Now we consider the uncorrected matrix error.For the standard MZI mesh, this is E0 = 2 √ N µ [17].Using the transfer matrix of the standard MZI to first order in (α, β), the norm of the matrix error is: which is maximized when the MZI is in the cross state θ = 0.For the MZI+Crossing (Fig. 5(a)), we find: Up to irrelevant output phases, the effect of the crossing is to flip the relative sign of α and β, so the component errors appear anticorrelated.As a result, ∆T MZI+X ∝ sin(θ/2)µ, which is zero for the cross state.The actual error is found by adding the ∆Tmn in quadrature and averaging over the probability distribution Pmn(θ) = n sin(θ/2) cos(θ/2) 2n−1 (equivalent to Eq. ( 10)): For a wavelength-dependent splitter error µ ≈ (dµ/dλ)∆λ, the tuning range and bandwidth can be calculated from the expressions for Ec (Eq.( 27)) and E0 (Eq.( 31)), respectively: the tuning range is the range over which Ec(λ) < Emax, while the bandwidth is the range over which E0(λ) < Emax: From these expressions, we derive the enhancement factors reported in Eqs.(7) and Table 1.

Neural Network Model
The optical neural network model is based on the architecture described in Ref. [11].Input images are first Fourier transformed, and cropped to a √ N × √ N window, where N is the DNN's inner layer size.The signal from this window (N input neurons) passes through two optical layers, with unitary connectivity realized with rectangular meshes.The activation function at the inner layer is realized electro-optically: a fraction of each output field is tapped off and sent to a detector, whose photocurrent modulates the remaining output light [28,33], implementing the activation function: where α is the power tap fraction, g is the modulator response, and φ is the phase at zero power.Here, we choose α = 0.1, g = π/20, and φ = π, so that f (E) approximates a leaky ReLU in the right power regime.Models of sizes N = 64 and N = 256 were trained using the Neurophox package [40].

Simulations and Data Analysis
All simulations were performed using the Meshes package, an open-source simulator for feedforward photonic circuits that can account for hardware imperfections [64].Figs.

Supplementary Material S1 Error Correction Methods
Correction of hardware errors is performed using the nulling method, which is based on the diagonalization of a unitary matrix using Givens rotations.This is closely related to the QR decomposition for the Reck triangle [S1], and a related decomposition for the more compact Clements rectangle [S2].The original nulling proposal was restricted to triangular (Reck) meshes and used internal tap detectors to monitor the output power of each MZI [S3, S4].Subsequently, the method was extended to generic mesh types [S5], and Ref. [S6] showed that that external detectors are sufficient for both Reck and Clements meshes.
Following Appendix A of Ref. [S6], we describe here the nulling procedure for configuring a Reck mesh.First, we write the coupling matrix for the multiport interferometer as a product of the 2×2 MZI blocks and an external phase screen: Here, the T mn represent tunable couplers (MZI, 3-MZI, MZI+X, etc.) while D is a diagonal matrix encoding the output phase shifts.The T mn are ordered along rising diagonals as shown in Fig. S1 (nulling also works on falling diagonals [S6]).We start by initializing the mesh to approximately the cross state, Fig. S2(a).We keep track of two matrices (Fig. S2(b)): W = T N −1,1 . . .T 11 is the partial product of all configured MZIs, and X = U W † , where U is the target unitary.At the beginning, none of the MZIs are configured, so W = I and X = U .At each step, with an example shown in Fig. S2(c), we configure target MZI T mn , which updates W and X by Givens rotations W → T mn W , X → XT † mn .The target T mn is chosen to zero a particular off-diagonal element X ij .Subsequent MZIs are configured in a sequence that successively zeroes off-diagonal elements of X (Fig. S2(d)).If all MZIs are configured properly, at the end of the procedure, X is diagonalized so U = DW , and the output phases (elements of D) can be read off by inspection.Nulling specifies constraints on the target Givens rotation T mn , which zeroes an element of X by right-multiplication (Fig. S3).Assuming unitarity of all matrices, the zeroing of X ij implies that: i.e. the power at the top output is zero when the fields (−v, u) are input to the crossing.This is equivalent to the splitting-ratio condition s = ŝ, where s ≡ (T mn ) 11 /(T mn ) 12 is the splitting ratio of the crossing (Eq.( 2), main text), and ŝ ≡ u/v is the target value.The difficulty in this procedure lies in the difficulty of accurately realizing Tmn in practice, since the actual transfer matrix is a function of both the control parameters (θ, φ) m mixes the elements (u, v) of X and zeroes out the rightmost one.This is equivalent to the condition Eq. (S2).
and the unknown fabrication imperfections.Therefore, for the realized Givens rotation, in general s = ŝ, which will lead to errors in the realized matrix U .
There are three distinct variants of the nulling method that accommodate hardware errors to different degrees: (1) an in-silico approach that does not correct errors [S1, S2], (2) measurement-assisted nulling, which corrects errors provided that s does not fall within a forbidden region [S3, S6], and (3) an improved measurementassisted method that partially compensates for the "uncorrectable" errors arising from unrealizable splitting ratios.

S1.1 In-Silico
Assuming ideal hardware, there is a simple relation between (θ, φ) and T .For example, for the standard MZI, T = ie iθ/2 e iφ sin(θ/2) cos(θ/2) in the absence of hardware errors.Using this formula, we can easily obtain (θ, φ) from the target splitting ratio: Following this procedure, the phase shifts of the mesh are found entirely in a computer.As a result, hardware errors are not accounted for when programming the mesh, and the realized matrix will be off by an amount called the uncorrected error: where • is the Frobenius (L 2 ) norm, U and Û are the realized and target matrices and ∆T mn = T mn − Tmn is the difference (due to hardware errors) between the realized T mn and the ideal Tmn given by Eq. ( S3).
In-silico methods were presented in Refs.[S1, S2] for the Reck and Clements meshes.The effect of hardware errors was studied in Refs.[S6, S7].Fig. S4(a) shows the flowchart for programming a mesh via in-silico nulling.

S1.2 Measurement-Assisted
In measurement assisted nulling, the first two steps are the same: find the target splitting ratio and update X and W using the corresponding Givens rotation.The principal difference is that (θ, φ) are found using an in-device measurement.For the Reck mesh, the procedure is traced out in Fig. S2, where each step attempts to zero a matrix element X ij by injecting w * j as input and adjusting the phase shifters to zero the output power at port i (see also Fig. S4(b)).This method was first proposed [S3] and demonstrated [S4] on the Reck mesh, but can be generalized to other mesh types provided that tap detectors are present after every MZI [S5].It was later shown that selfconfiguration is possible without the tap detectors [S6].Errors occur whenever a crossing cannot be programmed to reach the target splitting ratio, i.e. when ŝ lies within the forbidden region due to hardware imperfections.

S1.3 Improved Measurement-Assisted
In this paper, we have developed a refinement to the measurement-assisted nulling algorithm that allows for some of the "uncorrectable" errors to be partially corrected in subsequent nulling steps.The impetus for this refinement is the observation that, whenever uncorrectable errors occur, the s = ŝ, and the conventional algorithm as implemented in Fig. S4(b) incorrectly updates X and W . Error correction can be improved if we can accurately estimate the realized splitting ratio s; this allows the algorithm to use this information in order to partially compensate for such errors during the programming of subsequent MZIs.
The refined error correction algorithm is shown in Fig. S4(c).Here, we defer updates to X and W until the end, and after (θ, φ) have been set, we measure s through the following procedure: • If the output power is successfully nulled (P i = 0), then the coupler is configured correctly and s = ŝ.
• If P i = 0, nulling is imperfect and s = ŝ.To find s, we now perform an optimization: injecting wj (s) (the j th column of W (s) = T (s)W , which is a function of s), we vary s in the vicinity of s = ŝ (with the fixed (θ, φ) obtained in the previous step) until the output power is exactly zero.This procedure obtains the actual splitting ratio implemented in the tunable coupler.
Once s is found, we update X and W using T (s).Since the W and X updates are exact even in the presence of uncorrectable errors, the final matrix error is directly related to the residuals left by imperfect nulling of X.These residuals were calculated in the main text using

S1.4 Local Correction Method
For comparison, we also describe the local method for hardware error correction, first presented in Ref. [S8].This method is based on the principle that 2 × 2 unitary matrices are equivalent up to output phases if they share a common splitting ratio s ≡ T 11 /T 12 : ) This equivalence principle allows perfect MZIs to be substituted for imperfect MZIs columnwise, performing correction locally at each coupler (although the procedure is not strictly local: each step depends on the phases ψ i accrued from Eq. (S7) in the previous step).Errors occur only when MZI splitting ratios are unrealizable.These "uncorrectable errors" are independent of each other and add up in quadrature.Refs.In the notation of this paper, E c takes the form: Eqs. (S6) and (S9) are almost identical, differing only by the factor of q mn = |u mn | 2 + |v mn | 2 in the former.Since q mn ≤ 1 due to the unitarity of X, Eq. (S6) will always give a lower matrix error.
Table S1 lists the formulas for coverage (Eq.( 13), main text) and matrix error (Eqs.(S5-S6)) for the three mesh architectures and error models.We see that, for uncorrelated errors, only the 3-MZI is asymptotically perfect, while both the 3-MZI and MZI+X are asymptotically perfect for correlated errors.In addition, the uncorrected error can only be reduced in the correlated case, and only for the MZI+X.Finally, the examples of the 3-MZI and MZI+X highlight the superior performance of the improved self-configuration method.Under the original method, E c is independent of N , making the mesh types infinitely scalable (with respect to these errors) but not asymptotically perfect.But under the improved method, E c ∝ log(N )/N , which vanishes in the limit N → ∞.

Model
Eq. (S5)  3) and the local error correction method, Sec.S1.4.For the MZI and 3-MZI, Gaussian error models are used with σ = 0.05 and σ = 0.10, respectively.For the MZI+X, a correlated error model with µ = 0.1 is used.Dashed lines correspond to the analytic models in Table S1.
Fig. S5 plots the numerically computed accuracy on the three mesh types.For the 3-MZI and MZI+X, the difference in scaling with N is very clear.For the regular MZI, all methods give the same scaling, but self-configuration leads to an error lower by a factor of 2/3 ( 2/3N σ 2 vs. (2/3)N σ 2 ).The overall error amplitude in the figure is distorted by saturation when E ∼ 1, but the factor of 2/3 is still clearly apparent.

S1.5 Comparison to Global Optimization
Before the self-configuration and local algorithms were developed, the only way to train imperfect meshes involved global optimization [S9   S7 compare the accuracy of the self-configured solution to this global refinement.Interestingly, the improvement is fairly significant (3-4×) for Clements, but negligible for Reck.We speculate that this discrepancy may be attributed to the triangular structure of Reck, where the MZIs near the apex of the triangle are most likely to lead to uncorrectable errors.Since the upper-left corner of the matrix depends only on these MZIs, errors in this region cannot be corrected by adjustments to MZIs up-or downstream.This is in contrast to the Clements mesh, where all paths pass through an equal number of MZIs, and errors in the center of the mesh (where the probability density clusters close to the cross state) can potentially be corrected by adjustments near the edges.
In both meshes, up to a constant factor, the selfconfigured and globally-optimized solutions have the same error scaling E c ∝ N σ 2 in the MZI mesh.For large mesh sizes, the 3-MZI mesh still offers a significant improvement over the globally optimized solutions, and its E c ∝ log(N )/N scaling means that this gap grows larger with increasing mesh size.

S2 Imperfectly Correlated Errors
The splitter errors (α, β) of an MZI are best characterized by measuring the device's extinction ratio.To do so, one tunes the internal phase shifter θ and measures the con-  In most photonic platforms, splitter errors are strongly correlated so that ER cross ER bar .For example, in Fig. S8, we plot a histogram of measured MZI extinction ratios characterized for a 3-layer silicon-photonic neural network chip reported in Ref. [S8].The median bar-and cross-port extinction ratios are 23 dB and 32 dB, respectively.This correlation between splitter errors originates    to N = 512.This suggests that error correction allows the MZI+X to be asymptotically perfect on all practical mesh sizes, as scaling to meshes of size N > 512 is likely prohibitively challenging due to chip area and loss constraints.
In order to exactly cancel the differential term α − β as required for very large meshes N > 1024, one can place a heater above the directional coupler [S29].While this scheme does come with the cost of an additional active component (putting it in the same complexity category as the "perfect optics" approaches [S28, S30]), such an MZI+X with coupler trimming is unique in that it enjoys natively broad bandwidth, enhancing the WDM capacity of the system, which may prove critical to achieving competitive performance in photonic computing applications [S31].

S3 Length and Area Estimates
Table 2 of the main text provides a rough comparison of the resource costs of various mesh architectures.In all cases, the "perfect optics" designs [S28, S30] require 1.5-2× more active components, an important near-term concern as the size of existing chips is often limited by electronic packaging [S35] or power dissipation from heaters [S36].Waveguide length (which limits loss and SNR [S37] and on-chip latency [S8, S38]) and chip area are also critical parameters, but depend on the implementation.
The approximate MZI dimensions of a range of photonic mesh platforms are reported in Table S4.Most SOI devices has similar sizes, although there is a wider range of phase-shifter lengths owing to design tradeoffs (longer thermo-optic phase shifters can be more energy-efficient in certain cases [S39] and the higher heater resistance reduces the required current, but such devices suffer from increased loss and/or higher drive voltages).Non-SOI platforms such as silicon nitride and lithium niobate can support shorter optical wavelengths and offer mechanisms for faster pure-phase modulation, but suffer from reduced integration density due to the weaker phase-shift mechanisms (e.g.Pockels [S34] or piezo-optomechanical [S27]), which require much longer phase shifters.In such platforms, the length and area reduction for the 3-MZI is particularly pronounced, as these figures depend primarily on the number of phase shifters and not the number of passive components.Table S4: Waveguide (WG) length and on-chip areal dimensions (w × h) of phase shifters and beamsplitters on several published photonic platforms.The corresponding unit-cell length and area A = wh (normalized to the standard MZI) are computed from these dimensions.† MZI+X will have a size similar to 3-MZI.

Figure 2 :
Figure 2: Nulling method of self configuration.(a) Configuring MZI Tmn updates matrix W .(b) Corresponding nulling update to X = U W † , which is (c) equivalent to zeroing an output of Tmn given a fixed input.(d) Allowed range of s = T11/T12 ∈ C; regions near s = 0 and s = ∞ are forbidden due to imperfections.Contours are lines of constant (θ, φ), with α = 0.23, β = 0.07.(e) Probability density P (s) as a function of mesh size.

Figure 3 : 3 -
Figure 3: 3-splitter MZI design and simulated performance.(a) Schematic of 3-MZI.(b) Splitter Möbius transformation on s ∈ C, which pushes the forbidden regions away from s = {0, ∞}, corresponding to a Riemann sphere rotation.(c) Matrix error E0, Ec as a function of splitter variation σ (fixed N = 256), comparing the standard and 3-splitter MZI designs.(d) Scaling with mesh size N (fixed σ = 0.05).(e) Matrix error as function of both N and σ, showing the sharp onset of "perfect" error correction in regions where the coverage C is of order unity.

Figure 4 :
Figure 4: Effect of hardware errors on DNN inference.(a) Benchmark neural network consisting of FFT preprocessing, windowing, and two DNN layers, where the linear connections U1 and U2 are realized with MZI meshes [17, 33].(b) Inference accuracy as a function of MZI error.

Figure 5 :
Figure 5: MZI+Crossing architecture.(a) Schematic of MZI+X.(b) Effect of the crossing is to flip the s = 0 and s = ∞ forbidden regions.For correlated errors, the forbidden region around s = 0 disappears.(c) Riemann sphere projection.

Figure 6 :
Figure 6: Advantages of MZI+Crossing architecture for correlated component errors.(a) Matrix error as a function of µ for fixed N = 256.(b) Dependence on size N for fixed µ.

Fig. S2 traces out the nulling steps for a 4 × 4
Fig.S2traces out the nulling steps for a 4 × 4 Reck mesh.We start by initializing the mesh to approximately the cross state, Fig.S2(a).We keep track of two matrices (Fig.S2(b)): W = T N −1,1 . . .T 11 is the partial product of all configured MZIs, and X = U W † , where U is the target unitary.At the beginning, none of the MZIs are configured, so W = I and X = U .At each step, with an example shown in Fig.S2(c), we configure target MZI T mn , which updates W and X by Givens rotations W → T mn W , X → XT † mn .The target T mn is chosen to zero a particular off-diagonal element X ij .Subsequent MZIs are configured in a sequence that successively zeroes off-diagonal elements of X (Fig.S2(d)).If all MZIs are configured properly, at the end of the procedure, X is diagonalized so U = DW , and the output phases (elements of D) can be read off by inspection.

Figure S5 :
Figure S5: Comparison of the accuracy of self-configuration (SC, Sec.S1.2-S1.3)and the local error correction method, Sec.S1.4.For the MZI and 3-MZI, Gaussian error models are used with σ = 0.05 and σ = 0.10, respectively.For the MZI+X, a correlated error model with µ = 0.1 is used.Dashed lines correspond to the analytic models in TableS1.

Figure S6 :
Figure S6: Comparison of self-configuration and global optimization for N = 256 meshes of the Reck (top) and Clements (bottom) topology, with uncorrelated errors.

Figure S7 :
Figure S7: Dependence of corrected matrix error Ec on mesh size.Model: uncorrelated splitter errors with σ = 0.05.

Figure S9 :
Figure S9: Effect of finite cross-port extinction ratio on corrected error for MZI+X; compare Fig. 6 (main text).(a) Dependence of matrix error on µ for a Reck mesh of fixed size N = 256.(b) Dependence on N for fixed µ = 0.1.

Table 2 :
Characteristics of the major tunable crossing types.

Table S2 :
Reported bar-and cross-port MZI extinction ratios.Extinction ratios are reported in dB.