The Goldilocks Principle of Learning Unitaries by Interlacing Fixed Operators with Programmable Phase Shifters on a Photonic Chip

Programmable photonic integrated circuits represent an emerging technology that amalgamates photonics and electronics, paving the way for light-based information processing at high speeds and low power consumption. Programmable photonics provides a flexible platform that can be reconfigured to perform multiple tasks, thereby holding great promise for revolutionizing future optical networks and quantum computing systems. Over the past decade, there has been constant progress in developing several different architectures for realizing programmable photonic circuits that allow for realizing arbitrary discrete unitary operations with light. Here, we systematically investigate a general family of photonic circuits for realizing arbitrary unitaries based on a simple architecture that interlaces a fixed intervening layer with programmable phase shifter layers. We introduce a criterion for the intervening operator that guarantees the universality of this architecture for representing arbitrary $N \times N$ unitary operators with $N+1$ phase layers. We explore this criterion for different photonic components, including photonic waveguide lattices and meshes of directional couplers, which allows the identification of several families of photonic components that can serve as the intervening layers in the interlacing architecture. Our findings pave the way for efficiently designing and realizing novel families of programmable photonic integrated circuits for multipurpose analog information processing.

Levering the unique properties of light to perform computations in novel ways is a subject with a long history 1,2 .Although an all-optical processor for universal computing seems to be a far reach goal, photonics can provide exciting opportunities for unconventional computing building on analog logic and with novel information processing configurations.What makes photonics an intriguing option for unconventional computing are exotic potentials such as the intrinsically high speed and low energy consumption, the capabilities for massive parallelization, and long-range interactions.Nevertheless, a significant challenge in optical computing lies in the absence of appropriate computing paradigms, methods, and algorithms that harness the unique capabilities of this technology to develop efficient and application-specific photonic processors.In particular, matrix-by-vector multiplication is one of the most basic mathematical operations that lies at the core of various tasks ranging from optical convolution schemes 3,4 and matrix eigenvalue solver 5 to novel optical memristors 6,7 and optical artificial neural networks 8,9 .In the past decade, with rapid technological progress, there has been a resurrection in efforts devoted to developing such programmable photonic integrated circuit that performs matrix-vector multiplication [10][11][12] .The utility of such a device as an energy-efficient photonic accelerator in conjunction with electronic processors appears to be a distant possibility, considering the inherent difficulties associated with scaling and precision.However, there is no doubt that an on-chip programmable photonic matrix-vector multiplier can create exciting opportunities in classical and quantum computing through various applications that range from quantum information and quantum transport simulations [13][14][15] , to optical signal processing 16 , neuromorphic computing 17 , and optical neural networks 3,18,19 , as well as putting forward a platform for rapid prototyping of linear multiport photonic devices 20,21 .
Indeed, the optical realization of arbitrary unitary operations has been known since the seminal paper of Reck et al. 22 , which originally concerned free space optics but successfully translated to photonic integrated circuits by Miller [23][24][25] .This architecture builds on breaking down unitary matrices of any order into lower-dimension unitary matrices, which ensures the existence of an optical realization through two fundamental building blocks that are beam splitters (couplers) and phase shifters.Despite its generality, this method uses a pyramid-shaped array of Mach-Zehnder interferometers (MZI), which is impractical for larger implementations because the number of beam splitters grows quadratically with the number of ports.In turn, Clements et al. 26 introduced an alternative and symmetric rectangular-shaped array, resulting in a device with half the total optical depth and, consequently, more loss-tolerant.Such a rectangular array has been proved robust enough to create photonic realizations of Haar-random matrices 27 .Further unitary realizations related to other mesh geometries have been explored in 28,29 , as well as topological photonic lattices with hexagonal-shaped arrays of MZI 30,31 .In turn, free-space propagation setups have been devised based on plane-light conversion 32,33 and using diffractive surface layers 34 .
While the latter devices originally consisted of bulky optical components, the principle has recently been applied to on-chip structures as well.Recent studies have explored the use of particular transfer matrices (henceforth called F) alternating with phase mask layers to obtain an arbitrary unitary transformation 10,12,[35][36][37] .Pastor et al. 12 considers wave propagation in multimode slab waveguides to implement a Discrete Fourier transform (DFT) as their transformation F. They showed that an arbitrary transformation could be performed when 6N + 1 phase layers and 6N DFT elements, where N is the number of ports.Tanomura et al. 35 interleave the phase masks with multimode interference couplers connected with single-mode waveguides and use simulated annealing optimization to argue for well-approximated conversions when M ≈ N. Fully functional unitary four-, eight-, ten-, and twelve-port devices have been proposed and manufactured [38][39][40][41] .Moreover, an alternative device using polarization and multiple wavelength degrees of freedom has been considered in 42 .Markowitz and Miri 37 have explored similar structures and have found rigorous numerical evidence that interleaved phase arrays and discrete fractional Fourier transform (DFrFT) are indeed universal, whereas the use of Haar-random unitary matrices has been proved to lead to the desired universality 10 .A further waveguide array with varying propagation constants with step-like profiles has been reported 43 .
The interlacing architecture appears to exhibit interesting auto-calibrating properties, which makes it resilient to fabrication errors 44 .Furthermore, we recently have shown that the intervening structure can go beyond implementing unitaries to directly implementing arbitrary non-unitary operations when the diagonal matrices are relaxed to leave the unitary circle in the complex domain 45 .In this sense, by utilizing both amplitude and phase modulations, one can realize a fully programmable device for arbitrary matrix operations 45 .This shows an important generalization of the previous results that show by itself an advantage of the interlacing architecture over the mesh geometries.
This manuscript discusses a broad class of universal on-chip photonic architectures based on a layered configuration of phase mask layers as programmable units interlaced with a passive random matrix F. It is shown that the proposed interlaced architecture is far more flexible by showing that broad families of matrices F can serve as the fixed intervening operator.The phases are steered to reconstruct a unitary target matrix, provided that F has well-posed properties.Numerical evidence based on rigorous optimization algorithms reveals that universality is reached for dense matrices F, while a phase transition in the accuracy of reconstructed N × N targets occurs at M = N + 1, with M the total number of phase mask layers.Tests using the discrete Fourier transform (DFT) and discrete fractional Fourier transform (DFrFT) confirm the latter claim, and Haar-random matrices also show outstanding convergence.To generalize the domain of the valid intervening operators, a density criterion is derived so that matrices F can be classified according to their elements to ensure universality.To demonstrate this result, photonic lattices with uniform, nonuniform, and disordered coupling coefficients are considered as photonic realizations for the matrix F, while using the proposed density criterion, it is shown that universality is reached for specific length intervals.Furthermore, we explore waveguide coupler meshes as an alternative intervening unit F and determine the minimum number of coupler layers required to guarantee the universality of the interlacing architecture.

Architecture and Mathematical Foundation
Let us consider an arbitrary unitary matrix U ∈ U(N), with U(N) the group of unitary matrices N × N. Our goal is to implement a proper factorization of U in terms of another unitary matrix F to be defined and a set {P k } M k=1 composed of phase matrices n ∈ (0, 2π] for n ∈ {1 . . ., N} and k ∈ {1, . . ., M ∈ N}.The factorization proposed here is such that it intercalates F with a phase matrix P k through the relation This factorization is convenient for optical applications as the phase matrices can be implemented through phase shifter layers, which are the active optical elements of the architecture.In turn, the matrix F (passive optical element) has to be selected so that arbitrary unitary target matrices U t can be reconstructed with minimal error by adequately tuning the phase shifters φ If the latter is achieved, it is said that the universality property has been met.An arbitrary matrix U ∈ U(N) requires N 2 real parameters to be fully defined.Therefore, while performing the factorization, it is vital to consider at least the same number of parameters.For the device proposed in (1), we have MN free parameters in total and expect that M ≥ N in order to achieve the desired universality.Although there are some cases where U has a particular symmetry that reduces the number of parameters, we aim for the general case.The proposed architecture (1) and its optical implementation are sketched in Figure 1.On the one hand, the mathematical structure of the unitary F (see top panels of Figure 1) can be that of the discrete Fourier transform (DFT), discrete fractional Fourier transform, or simply a Haar-random matrix.On the other hand, the photonic implementation of F can be performed through different approaches (see bottom panels of Figure 1), such as waveguide arrays 45,46 , meshes n }, with p = 1, . . ., N + 1.The upper insets depict the modulus and argument of the potential candidates for the unitary matrix F, which have been selected as the DFT, DFrFT, and a random unitary matrix.The lower insets illustrate potential photonic implementations to perform the unitary matrix F. of directional coupler 26,28,29 , and multimode interference (MMI) 12,47 .Here, we focus on architectures based on the first two solutions, the numerical analysis and universality of which are discussed below.Layered architectures akin to Eq. 1 have been numerically validated in previous works using different optical arrays and MZI meshes, where different optimization algorithms such as gradient-descent, stochastic gradient descent, simulated annealing, and basin-hopping have been implemented.See for instance 10,40,43,48 .In order to demonstrate the universality of this device, we optimize the NM phases for a variety of randomly chosen target matrices U t generated in accordance with the Haar measure 49 .The objective function to be minimized, also called error norm, is defined by where ∥ • ∥ stands for the Frobenius definition of the norm, U t is the target matrix being tested, and U is the reconstructed matrix using the factorization (1).The Levenberg-Marquardt (LM) algorithm 50,51 is used to find the minimum of this function.
For a given target U t , the phases are randomly initialized between 0 and 2π.The optimization was performed in MATLAB.
In Figure 2(a), the norm for 100 target unitary matrices is shown for various cases, fixing the default tolerance values to 10 −10 .A phase transition occurs between the M = N and M = N + 1 layers, which is unsurprising given the system becomes over-determined by N parameters.These jumps are larger than reported by Tang et al. 52 due to their usage of a probabilistic algorithm (Simulated Annealing) rather than a gradient-based one such as LMA.A downside of gradient-based approaches is they may require many runs with different starting conditions.We can decrease the overhead by using a stopping criteria for the norm along with a maximum iteration for each run of LMA.Using a maximum iteration of 50, we find that we rarely need For systems with a lower number of ports N, the distribution skews towards lower values.To more confidently label choices of F which are not Haar-random generated as "bad" mixing layers, we set the maximum number of runs somewhat higher to 250 or 500 as found appropriate.

Density Estimation Criterion
Our preliminary numerical results suggest that dense matrices F lead to a universal architecture, and it is thus desirable to find a proper measure to quantify the density for unitary matrices.Indeed, random matrix theory establishes robust criteria for studying random complex-valued matrices at the limit of large dimensions based on the singular value decomposition (SVD) of the same 53 .On the one hand, in the current setup, we focus on unitary matrices of relatively small size, as the number of ports in our architecture is not necessarily long enough to take into account the asymptotic analysis of random matrix theory.On the other hand, the singular value decomposition of unitary matrices is not particularly helpful when dealing with unitary matrices, as for a complex-valued matrix A, the SVD requires the computation of the spectral properties of AA † and A † A, which for any unitary matrix A is always equal to the identity matrix.For these reasons, the density criterion discussed below suits better for the interlaced architectures here constructed.
To better understand the importance of dense matrices in our universal architecture, let us further inspect the factorization in (1).By fixing F = diag(e iξ 1 , . . ., e iξ N ) as a diagonal unitary matrix, for ξ j ∈ (0, 2π] for j = 1, . . ., N, it is straightforward to notice that (1) reduces to a diagonal unitary matrix as well, which is far from representing a universal device.That is, diagonal F matrices can only reconstruct diagonal unitary matrices.This behavior extends to any of the N! 2 different permutations allowed to the diagonal matrix F, as the factorized matrix U acquires the same structure as that of the permuted F matrix.Now, let us consider the set of unitary matrices {V n j } k j=1 , where V n j ∈ U(n j ) and ∑ k j=1 n j = N, such that F = diag(V n 1 , . . ., V n k ) is an N-dimensional block-diagonal unitary matrix.This particular selection for F leads to a unitary operator with the same structure, which lacks the required universality property.Although block-diagonal matrices F lose their structure when randomly rearranged, anti-diagonal block matrices will always result in either block-diagonal or anti-block-diagonal matrices, neither of which are universal.
Thus, we have identified a particular class of bad-performing matrices interlacing matrices F. Still, the latter is still useful to trace a suitable definition of density.To this end, let us first remark that any N-dimensional unitary matrix can be written as V = (⃗ v 1 , . . .,⃗ v N ), where ⃗ v j ∈ C N are complex-valued column vectors that form an orthogonal set through the Euclidean inner product in C N , i.e., ⃗ From the orthogonality condition, we can simply focus on the columns (or equivalently the rows) of V , whereas the normalization imposes a constraint on elements across the columns (rows).Since the elements of V are complex numbers, we will work with the alternative matrix V which is composed of the modulus of the elements of V .Let us define v p;q as the q-th element of the p-th column vector ⃗ v p , with p, q = 1, . . ., N, so that V p;q = |v p;q |.From the unitarity of V , it follows that ∑ N q=1 |v p;q | 2 = ∑ N p=1 |v p;q | 2 = 1 for all p, q, so that we can focus on the density of either the columns or rows of V .Without loss of generality, we work with the columns.Since we are interested in how the elements are sparse across each column, we compute the corresponding variance for each column p.Following the normalization condition, it is straightforward to prove that the variance is a bounded quantity in the interval S p = 0, N−1 N 2 , where the lower bound corresponds to the case where all the elements of ⃗ v p are equal, i.e., ⃗ u p is maximally spread.The upper bound corresponds to the case where ⃗ v p is one of the canonical unit vectors (0, . . ., 1, . . ., 0) T .Thus, a given column p of V is said to be denser if its corresponding variance S p approaches (N − 1)/N 2 .The more sparse the elements of the column p, the more S p approaches 0. These are the key ideas we use henceforth to characterize density across the full matrix V .
Be S = {S p } N p=1 the set of variances associated with each column of V .We define the mean µ and standard deviation σ associated with the elements of S , so that the density of a given unitary matrix V can be characterized by defining the point Since row permutation leaves S p invariant and column permutation only permutes the index p, the quantities µ and σ , and consequently the point ⃗ R, are permutation invariant.
There are two note-worthy extremal cases, namely the maximally sparse and the diagonal cases (densest cases) unitary matrices.In the former case, V is composed of column vectors so that the variances vanish, S p = 0, for all p = 0, . . ., N. The DFT matrix of dimension N is such an example.This leads to µ = σ = 0, which we consider the ideal case.For the second case, the variances are maximal per each column, S p = (N − 1)/N 2 , for all p = 1, . . ., N, and the statistical information of the matrix reduces to µ = (N − 1)/N 2 and σ = 0. We thus have two comparison points, from which we find the bounded interval Additional reference points can be traced out if we take the block-diagonal matrices F = diag(V n 1 , . . ., V n k ), with V n j ∈ U(n j ) unitary and maximally sparse (DFT) matrices of dimension n j for j = 1, . . ., k, 1 ≤ k ≤ N, and ∑ k j=1 n j = N, so that badperforming matrices are generated (see discussion above).Particularly, let us consider the case k = 2 so that F = diag(V n 1 , V n 2 ), with n 1 + n 2 = N.One can assign the indexes n 1 = ℓ and n 2 = N − ℓ, with ℓ = 1, . . ., ⌊N/2⌋ nonequivalent ways to define the two block-diagonal matrices F that leads to the k 2 reference points Indeed, further reference points exist for k ≥ 3.However, those points are farther from the ideal (maximally sparse) case ⃗ R 0 ≡ ⃗ R k=1 than those marked with k = 2.In this form, we can just focus on the area spanned between k = 1 and k = 2. Interestingly, for k = 2 and ℓ = ⌊N/2⌋, one obtains the reference points with smaller standard deviation, which are ⃗ R k=2,ℓ= N for even and odd N, respectively.Note that in the limit N → ∞, the mean converges to the non-vanishing value 1/2.In turn, the maximum standard deviation is determined by minimizing N σ in terms of ℓ, from which one obtains the critical value and the maximum standard deviation N σ | ℓ c = 1/4.That is, the standard deviation is bounded to the interval N σ ∈ [0, 1/4], where the upper bound is independent of N and is given by max N σ | ⌊ℓ c ⌋ , N σ | ⌈ℓ c ⌉ .The latter allows us forming a polygon with vertices at ⃗ R k=1 = (0, 0) and ⃗ R k=2 , the area of which is non-null and finite even for N → ∞ (see Figure 3).We focus on unitary matrices whose vector ⃗ R lies inside the latter polygon while avoiding the vertices, as the latter are the well-known bad-performing cases (with the exception of ⃗ R 0 ).We can go a step further and make a better prediction of unitary matrices by reducing the area of the polygon and imposing a threshold to N σ and N⃗ µ.For the former, we already know that the standard deviation reaches its maximum value at N σ | k=2,ℓ σ = 1/4 for all N. We thus implement the threshold at half of the maximum allowed standard deviation, i.e., N σ th = 1/8.For even N, the reference points ⃗ R k=2,ℓ≥ℓ σ are above such a threshold for ℓ σ = N 2 − ⌈ N 4 2 − √ 3⌉.Likewise, we fix the threshold for the mean at Therefore, any unitary matrix F whose associated vector ⃗ R lies inside the interception of the polygon spanned by the set of ℓ=1 and the thresholds N µ th and N σ th has statistical properties in its matrix elements to be considered sparse enough to induce universality in our architecture.This provides a tool to quantify and classify a priori the F matrices.

Random matrices and performance test
To test this measure, we randomly generate families of 6 × 6 unitary matrices and determine their location in the plane (N µ, N σ ).The random matrices are generated by using the decomposition F X = e iX , with X † = X a Hermitian matrix in C N×N to be defined.In order to preserve control over the density of the randomly generated matrices, we introduce the four non-symmetric matrices X A = (⃗ χ 1;A , ⃗ 0, ⃗ 0, ⃗ 0, ⃗ 0, ⃗ 0), X B = (⃗ χ 1;B ,⃗ χ 2;B , ⃗ 0, ⃗ 0, ⃗ 0, ⃗ 0), X C = (⃗ χ 1;C ,⃗ χ 2;c ,⃗ χ 3;C , ⃗ 0, ⃗ 0, ⃗ 0), and D = (⃗ χ 1;D ,⃗ χ 2;D ,⃗ χ 3,D ,⃗ χ 4;D , ⃗ 0, ⃗ 0), with ⃗ 0 the null-vector in C N .The column vectors ⃗ χ j;℘ are composed of zeros in the first j inputs and random numbers elsewhere, with ℘ ∈ {A, B,C, D}.We thus consider the Hermitian construction as X = X ℘ + X † ℘ , with ℘ ∈ {A, B,C, D}.Note that in each case, the number of random parameters increases in each case as additional columns are included, which implies that more parameters are available to build up the corresponding unitary matrices.We thus expect denser matrices constructed from D than A.
In this form, we establish a controlled testing scenario for our density criterion and the subsequent universality of the matrix under consideration.These results are depicted in Figure 3(a), where the reference points ⃗ R k and the points ⃗ R associated with the sets of random unitary matrices F X ℘ for ℘ ∈ {A, B,C, D}, so that we can determine the unitary matrices whose corresponding points ⃗ R lie inside the universality area.On the one hand, for the random unitary matrices F A , only one of the points, A 5 , is expected to reveal universality properties.On the other hand, we expect more matrices F D with such a universality, such as D 1 , D 5 , and D 12 , to name a few; whereas the matrices associated with D 3 , D 10 , and D 19 are expected to lack the desired property.This criterion is tested by using the latter unitary matrices, F ℓ , and testing their performance for randomly generated unitary target matrices.Here, fifty target unitary matrices are considered for each testing unitary matrix so that the relative error of the so-reconstructed targets can be analyzed for a broad number of cases (see Figure 3(b)).Such a performance test reveals a good match with our universality criterion.For instance, the error associated with F A 2 and F A 6 is relatively high for each of the target matrices tested; in turn, F A 5 performed well in each scenario, and F A 13 and F A 20 have revealed a handful of cases with a higher error.These two matrices lie far outside the universality area shown in Figure 3(a), and thus a bad performance was expected.The performance for the unitary matrices F D j reveals a better fit across all the cases, where only F D 3 , F D 10 and F D 19 contain a few cases in which the target is reconstructed with a moderate error around 10 −5 .
As the number of input channels in an architecture increases, the time required for performance testing also increases.Therefore, the density criterion allows us to classify and select in advance the required F matrices that are expected to lead to universal behavior.Another advantage is that the proposed universality area spanned in the ⃗ R-space is finite and non-null for N → ∞, making it a suitable measure for architectures with an arbitrary number of ports.Thus, in Figure 3(c), we depict a more convenient representation of our results in which we represent both N µ and N σ separately, where we determine universality if both quantities lie below their respective thresholds.This permits a better and more organized reading for each matrix F under consideration.Remark that, although our density criterion determines the boundaries where universality is likely, it may exclude some well-performing cases, i.e., this is a sufficient but not necessary condition to identify universality.

Photonic Realizations
So far, the universality of the proposed architecture has been successfully determined using Haars random F matrices possessing the required density criterion.We now proceed to analyze potential candidates for photonic realizations of such F matrices, with a primary focus on photonic lattices and meshes of directional couplers.

Waveguide Lattices
Waveguide arrays can be modeled with high precision using coupled-mode theory 54 .The latter takes into account the coupling of evanescent waves from one waveguide interacting with its nearest-neighbor while neglecting farther neighbors due to their weak coupling.In this form, the effective Hamiltonian describing an array of N waveguides is characterized by a tridiagonal and symmetric matrix Hamiltonian H of dimension N. The wave evolution through the lattice is ruled by the dynamical law , where ⃗ u(z) ∈ C N is the electric field at each waveguide at the propagation distance z.Since H ̸ = H(z), the wave evolution is determined through the unitary evolution operator F(z) = e −izH as ⃗ u(z) = F(z)⃗ u(z = 0).We can thus implement waveguide arrays in the universal architecture, provided they fulfill the desired universality.To this end, we can test the behavior of a given lattice evolution operator for specific lengths using the density criteria of Section .Particularly, we consider the photonic J x lattice 46 , the homogeneous lattice 55 , and the disordered homogeneous lattice 56 as the physical waveguide arrays under consideration described by the respective Hamiltonians H (J x ) , H (h) , and H (h,d) .The matrix elements of the latter are explicitly given by  with p, q ∈ {1, . . ., N}.Here, κ(p) = κ 0 2 (N − p)p stands for the coupling parameter between nearest waveguide neighbors in the J x lattice, whereas ∆κ p ∈ N(µ, σ ) are random numbers taken from the normal distribution N(µ, σ ) characterizing the disorder effects.The coupled waveguide implementation for each Hamiltonian is depicted in Figure 4.

6.1544
The corresponding unitary evolution operators are simply given by F (J x ) (z) = e izH (Jx) and F (h) (z) = e izH (h) .Although both are functions of the lattice length z, the J x lattice (Figure 4(a)) has equidistant eigenvalues that lead to a periodic unitary evolution operator F (J x ) (z) in z, so that we can simply focus on the interval z ∈ [0, 2π).We first estimate the lengths that induce universality in our architecture, which is depicted in Figure 4(b).In the latter, we mark the particular lengths z (m) j and z (M) j that denote the local minima and maxima of the standard deviation N σ , the exact values of which have been determined numerically and presented in Table 1.In turn, the black-thick line in Figure 4(b) denotes the lengths where both N µ and N σ are below the universality threshold; i.e., the lengths where universality is expected.Without any prior performance test, one can see that the lengths z 1 , and z (M) 2 may fail in obtaining the desired universality.Recall that our estimation criterion is only a sufficient condition and may rule out positive cases.Nevertheless, all the other marked points can be considered candidates for the matrix F in our architecture, as no false positive cases will be included.This is indeed verified in the performance test portrayed in Figure 4(c), where, for each point, we have used fifty randomly generated unitary matrices as targets.The latter confirms 8/13 1 , reinforcing the fact that only two positive cases were discarded, but no false positives were included.In this form, we can confidently conclude that a universal architecture can be built using J x lattices with lengths as small as z = π/4 for N = 10.
We alternatively consider the homogeneous lattice, which contains waveguide arrays homogeneously distributed (Figure 4(d)).The eigenvalues accordingly distributed as λ (h) n = 2κ 0 cos( nπ N+1 ), and no periodic behavior is expected.We thus focus on the interval ℓ ∈ [0, 4π] for this particular lattice.The density criterion shown in Figure 4(e) reveals that lattice lengths in the interval z/κ 0 ∈ (2.4098, 5.9595) are suitable for our universal architecture.Particularly, note that the interval z/κ 0 ∈ [ℓ (m) 4 , ℓ (m) 6 ] contains lengths so that N µ and N σ remain mostly constant with minor variations.Thus, the performance test in this interval is expected to perform well.In the universality region, there is a local minimum ℓ (m) 2 that is isolated and associated with a shorter lattice length.This reference point may be useful for reducing the size of the universal structure.Figure 4(f) displays the corresponding performance test, which supports our previous statements.As expected, the performance for lattices with length ℓ (m) 1 and ℓ (M) 8 is particularly poor.However, the length ℓ (M) 1 has a generally good performance, with only two test targets displaying slightly higher errors than the other well-performing cases.
We additionally take into account the effects of disorder on the homogeneous lattice, which may be caused by impurities or imperfections during the manufacturing process, resulting in waveguides not being in their ideal positions or displaying deviations in their sizes (see Figure 4(g)).The defects here considered are such that the nearest-neighbor interactions deviate from the ideal homogeneous lattice by a factor of twenty percent; i.e., the disorder couplings in (4) take values from the normal distribution as ∆κ p ∈ N(µ = 0, σ = 0.2κ 0 ).Although the lattice structure is modified in the latter disorder, the estimation of density does not differ significantly from the ideal case, as shown in Figures 4(h)-(i).

Directional Coupler Mesh
Alternatively, the matrix F can be optically realized through proper mesh arrays of directional couplers.Particularly, we consider a construction based on two-port passive elements, which act as a power divider (50:50 MZI) equivalent, up to a global phase, to the unitary matrix U(2) ∋ T 0 = 1 √ 2 (σ 0 + iσ 1 ), with σ j the conventional Pauli matrices for j ∈ {1, 2, 3} and σ 0 the 2 × 2 identity matrix.The latter can be used as a building block to construct other U(2) matrices 57 as well as higher dimensional unitary matrices U(N) through appropriate Kronecker products 22 .In this section, we consider the symmetric 10-port array portrayed in Figure 5(a), composed of power dividers T 0 interconnected through different layers L 1 and L 2 .
Here, each layer is described by the U(10) matrices L 1 = I 5 ⊗ T 0 and L 2 = I 1 ⊕ (I 4 ⊗ T 0 ) ⊕ I 1 , with ⊗ and ⊕ the Kronecker product and Kronecker sum 1 , respectively, and I n the n × n identity matrix.The p-layered unitary matrix describing the power 1 Such operations are also known as direct product and direct sum.

9/13
divider array is thus given by where we have truncated the maximum number of layers to ten.It is not mandatory to truncate these layers, and the procedure can involve additional layers if needed.However, for practical physical implementations and limitations, we aim for compact devices and need to limit the number of layers as much as possible.To this end, we estimate the number of layers p for which universality is achieved.The preliminary estimation displayed in Figure 5(b) ensures the goodness of F for p = 7, 8, 10, being p = 10 the case in which both N µ and N σ minimize.Based on the performance test shown in Figure 5(c), it was found that the p = 7, 8, 10 layers are capable of achieving a F capable of achieving the required universality in the architecture (1).Notably, optimization results indicate that our device could perform well with only p = 6 layers, which was originally deemed inadequate from the density criterion.That is, we obtain a false negative for F. Still, no false positives were detected during the analysis, which is strictly necessary to avoid faulty designs of the final architecture.

Conclusions
We have introduced the design for a lossless universal photonic architecture based on a layered scheme of interlaced active phase shifter layers and passive random matrices.Numerical results obtained from the LMA optimization revealed that generating Haar random matrices F leads, in a vast majority of cases, to the desired universal architecture.It is observed that well-behaved matrices F show a phase transition on the error norm L of the reconstructed target at M = N + 1, with M the total number of phase shifter layers.In such a layer number, the error drops significantly to numerical noise values.While this is not proof that the factorization is exact, the error involved in the reconstruction process lies in the numerical error regime, and it is thus low enough to ensure that any unitary matrix is reconstructed with the desired accuracy.
Despite the accuracy of the LMA optimization, the computational time required for testing the universality of the random matrices F scales with the total number of ports, which becomes impractical for particularly large architectures.Numerical evidence shows that denser matrices perform better than sparse ones, usually involving relatively large errors.Therefore, a density criterion has been devised and introduced to classify the candidates for the matrix F used in the architecture.This criterion is built on preliminary knowledge of bad-performing matrices, such as diagonal and block diagonal matrices, which are analytically known to fail but serve as reference points to look for good-performing matrices, such as the DFT case.In this form, instead of performing a long optimization routine on the candidate for F, we simply analyze the standard deviation of the modulus of its columns or rows, which provides information about its density.This allows defining a mapping ⃗ R : U(N) → R 2 , which renders a vector that estimates whether F is suitable for the architecture.We thus possess a tool to preselect matrices F beforehand, making the design process more practical than generating and testing several random matrices.
Our tests using randomly generated unitary matrices showed that matrices within the threshold marked by the density criterion led to the required universality.Thus, universality is not limited to a specific realization of F; as shown in the results, infinitely many unitary matrices can meet our requirements.This paves the way for more efficient construction and optimization of compact devices that are simultaneously resilient to random defects.Particularly, the photonic Jx lattice was found suitable for this task at lengths different than the previously reported critical value ℓ = π/2 37 .This defines intervals in the lattice length for which the architecture is universal, leading to more flexibility in the manufacturing process so that one can allow for deviations in lattice length.This fact is further supported in the context of homogenous lattices, which are also suitable for our architecture and robust against disorder effects due to waveguide impurities or mismatching sizes.The latter was tested by introducing deviations of up to 20% into the homogeneous lattice, from which the density estimation showed no significant difference in the universality performance for the lattice lengths considered.Further constructions for the F matrices are indeed allowed, and an alternative construction based on a layered array of power dividers was shown to be efficient for our purposes, the analysis of which allowed us to determine the optimal number of passive elements required for the architecture.
The presented results, in particular the density criterion, are not limited to the waveguide and direction coupler realization discussed in the manuscript.Indeed, any optical element described by a unitary matrix, such as multimode interference couplers, can be further tested by the density criterion and classified before being used in the layered architecture (1).The latter helps facilitate the design of the passive layer to reduce the overall architecture size, account for potential manufacturing errors, and diminish the device footprint.This is particularly handy when deploying more complex optical circuits.

Figure 1 .
Figure 1.Universal architecture scheme.The proposed architecture involving alternating layers of random unitary matrices F and diagonal phase shifts layers (PL) {φ(p)

Figure 2 .
Figure 2. Numerical universality test.(a) Architecture depiction (left column) and optimization objective function (right column) for 100 target matrices at various values of M and N. Black boxes denote any possible realization for the F matrix.(b) Multiple trials for N = 8 and M = 9 using 250 random F matrices were considered; 250 targets were used for each matrix F. Shown is the distribution of the number of LMA runs to achieve a norm lower than the stopping norm of 10 −10 , with a maximum of 50 iterations per run.(c) Norm (log 10 ) in terms of the number of iterations for the run with the best norm.Using 100 random matrices F, each with a single target matrix.

Figure 3 .
Figure 3. Density estimation and performance test.(a) Points ⃗ R associated with density estimation for the set of unitary matrices {e A j , e D j } 50 j=1 (left column) and {e B j , e C j } 50 j=1 (right column).The blue heat maps denote the absolute value, U , for some particular choices of unitary matrices.(b) Error norm (log 10 ) L in (2) for each unitary matrix under consideration with fifty testing targets per matrix.(c) Mean and standard deviation N µ and N σ , respectively, related to the density estimation for each unitary matrix in (a).The horizontal blue and red lines denote the universality threshold for N µ and N σ , respectively.

Figure 4 .
Figure 4. Photonic lattice realization and their universality.Sketch for the waveguide array associated with the J x lattice (a), homogeneous lattice (b), and homogenous lattice with disorder effects (c).Density criterion as a function of the lattice length ℓ (d-f) and the corresponding numerical performance test at the reference lengths ℓ (m) j and ℓ (M) j (g-i) for the J x lattice (left column), homogeneous lattice (central column), and disordered lattice (right column).
criterion N σ for the J x and homogeneous lattice.

Figure 5 .
Figure 5. Geometric array using power dividers.(a) Power divider array composed of p layers as defined in (5).Light-shaded and dark-shaded layers denote L 1 = I 5 ⊗ T 0 and L 2 = I 1 ⊕ (I 4 ⊗ T 0 ) ⊕ I 1 , respectively.Density criterion (b) and error norm (log 10 L) (c) of the mesh architecture in (a) as a function of the number of layers p.