Characterization of multi-mode linear optical networks

Multi-mode optical interferometers represent the most viable platforms for the successful implementation of several quantum information schemes that take advantage of optical processing. Examples range from quantum communication, sensing and computation, including optical neural networks, optical reservoir computing or simulation of complex physical systems. The realization of such routines requires high levels of control and tunability of the parameters that define the operations carried out by the device. This requirement becomes particularly crucial in light of recent technological improvements in integrated photonic technologies, which enable the implementation of progressively larger circuits embedding a considerable amount of tunable parameters. In this work, we formulate efficient procedures for the characterization of optical circuits in the presence of imperfections that typically occur in physical experiments, such as unbalanced losses and phase instabilities in the input and output collection stages. The algorithm aims at reconstructing the transfer matrix that represents the optical interferometer without making any strong assumptions about its internal structure and encoding. We show the viability of this approach in an experimentally relevant scenario, defined by a tunable integrated photonic circuit, and we demonstrate the effectiveness and robustness of our method. Our findings can find application in a wide range of optical setups, based both on bulk and integrated configurations.


Introduction
Linear optical networks are fundamental elements in several protocols for computation, communication and sensing. Recently, many schemes for computation have found their natural implementations through optical processing, such as neuromorphic and reservoir computing, 1-4 optical neural networks 5,6 and optical simulation. [7][8][9] Large-scale photonic platforms are one of the most promising candidates to implement quantum information and quantum computation protocols. 10,11 Indeed, they have been extensively employed for quantum walks routines, [12][13][14] quantum machine learning algorithms, [15][16][17] and, recently, for experiments that aim at demonstrating quantum advantage with photons. [18][19][20] All these protocols in classical and quantum optics require complex interferometric structures composed of numerous optical components. Integrated photonics is one of the best candidates to realize such optical protocols in compact devices, offering in addition the capability to reconfigure the circuit operation. 11,21 The latest examples of multi-mode optical networks 22,23 have shown a significant increase in the network complexity as well as in the number of control parameters.
The mathematical tool to model any linear optical processing in such experiments is the unitary matrix U that describes the relation between input/output complex amplitudes of the electromagnetic field. Several relevant aspects connected to optical network programming have been identified. They range from engineering the optical setup for implementation of a given U , to the identification of universal architectures that can perform any transformation. 24,25 Here we  Figure 1 Reconstruction of multi-mode optical circuits. a)The model of a multi-mode interferometer considered in this work. It is composed by the ideal optical circuit described by the unitary transformation U plus layers of mode-dependent losses at the input and at the output (represented by beam splitters in the figure), and phase instabilities (represented by sparks). Output losses take into account also possible differences in the detection efficiencies among the modes. b) Scheme for the measurement of second-order cross-correlations C hk ij with coherent light emitted by a continuous wave laser (CW). The latter is coupled in single-mode fiber and split into two beams by an in-fiber beam splitter. The two beams enter the interferometer in modes h, k. The phase modulation ϕ M performed by a liquid crystal compensates the fiber phase fluctuations ϕ = ϕ 1 − ϕ 2 to satisfy the conditions in Eq. (17). c) The 3-mode integrated chip employed to test the reconstruction algorithm. It is composed of a sequence of two tritter structures. Each tritter comprises three beam-splitters, whose reflectivity is reported in the figure, and a phase-shift equal to π/2. Between the two tritters there are three heaters {R 1 , R 2 , R 3 } that dynamically control the optical phases between the two structures via the thermo-optic effect.
focus on the characterization of linear optical circuits, which requires a systematic methodology to verify the operation of the optical circuit, or more generally, to reconstruct the unitary matrix implemented by the network. This task turns out to be a non-trivial one in certain scenarios. For example, interferometers based on a bulk optic implementation can be typically decomposed in elementary units that can be addressed and characterized separately. On the contrary, this procedure is typically not viable for integrated photonic circuits, and the network needs to be analysed as a single element.
Several techniques have been developed to achieve full characterization of integrated photonic devices by exploiting different degrees of knowledge of the network's internal structure and various measurement approaches with classical or quantum light. These algorithms can be divided into two main groups. The first one is called black box approach and exploits only the information provided by the measurements, without assumptions on the internal structure of the optical network. The second one is the white box approach, which exploits knowledge of the optical network structure together with the information obtained from a set of appropriately chosen measurements. White box approaches usually make use of an optimization algorithm to estimate the parameters of the given architecture. 26,27 It is clear that this family of reconstruction algorithms strongly depends on the knowledge of the relation between the optical component parameters and the matrix elements, and thus requires accurate modelling of the system including noise processes. On the contrary, black box approaches aim at characterizing a multi-mode linear optical interferometer without using any information about its structure. They are very useful when the structure is inaccessible, untrusted or too complex to be modelled. 23 A black-box method for linear optical circuits was proposed in Ref. 28 and subsequently refined in Ref. 29 The authors presented an analytical algorithm to reconstruct the elements of the unitary matrix U from quantum light measurement, via singleand two-photon experiments. The algorithm is advantageous since it permits to retrive the matrix U even in the presence of losses and optical phase instabilities due to, for example, fiber connections in the input and output stages. However, such a method requires quantum light input states, is a slow procedure and the accuracy is limited by noisy experimental data. In fact, the result of this reconstruction method is typically employed as a starting point for further numerical optimization to improve the stability of the solution. An alternative method exploits only classical field intensities measurements. 30 The moduli of the matrix elements are measured from the field intensities distribution in the output modes, while, the complex phases are estimated by a measurement of the interference fringes between two coherent beams. The two procedures for the moduli and the phases estimations are independent and they can be mapped directly onto the unitary matrix without the need to apply any further optimization algorithms. 26,30 With this approach, a crucial requirement for a correct phase measurement is to perform the phase scan in times much shorter than the typical timescale of phase fluctuations of the system. In addition, it cannot be used in the presence of mode-dependent losses in the collection stages. Other black-box algorithms based on coherent light measurements always require high phase stability during the scan in the input and output sections, 31-33 thus making them not viable for optical setups with in-fiber connections, which are nevertheless typical for integrated photonic devices. The last class we mention are machinelearning algorithms 34,35 that need large sets of data to learn the correct transformation. Very recently, Kuzmin et al. 36 simulated the application of a supervised-learning strategy for the calibration of a reconfigurable interferometer and experienced an unfavourable scaling of the training set size with the number of modes in the black box scenario.
In this work, we propose a new black box approach that aims at solving some open issues mentioned in the past algorithms. The goal is to provide a methodology to identify separately mode-dependent losses and matrix elements of U via coherent light measurements. In particular, we first present two algorithms to estimate the amount of loss unbalancement in the injection and collection stages of a linear optical interferometer and consequently characterize the moduli of the elements of the unitary matrix. Then, we move to the problem of measuring the phases of the unitary matrix elements. Previous algorithms 28, 29 exploited second-order quantum optical correlations such as the Hong-Ou-Mandel effect, 37 or the first-order classical correlations in Mach-Zehnder like interferometric structures. 30 Since there are mathematical analogies between classical and quantum second-order correlations, 38 we propose to replace such quantities with the second-order correlations of classical light in a Hanbury-Brown&Twiss like experiment. 39 This approach combines some advantages of both the previous methods, i.e. the simplicity in the use of classical light 30 and the independence from losses and optical phase instabilities, due to the fibers, characteristic of the methods that employ two-photon pairs as input states. 28,29 Moreover, we show that the proposed classical second-order measurements depend only on the matrix phases and not on the moduli as for the quantum correlations. This allows for completely independent estimations of the phases, input/output losses and the moduli of the unitary matrix. From an experimental point of view, this is an important improvement that reduces the propagation of the error along the characterization process. Thus, this method can be applied in a general scenario and can be effective for different experimental platforms, ranging from bulk to integrated and in-fiber optical setups.

Overview on black-box linear optical circuits reconstruction
Any ideal linear optical interferometer can be represented by a unitary matrix U , acting linearly on the annihilation (creation) operators of the electromagnetic field in the input modes a j and transforming to the annihilation (creation) operators in the output modes where U ij are the elements of the unitary matrix. The same relation holds for classical states of light, by replacing the operators a i and b i with the complex field E i in Eq. (1). In other words, the elements of the unitary matrix U ij completely characterize the field amplitudes propagation through a multi-mode optical network. In general, the set of optical modes i may represent any degree of freedom of light, such as polarization, path, time of arrival, frequency, angular and transverse momentum.
Optical losses deviate the interferometer operation from unitarity. To take into account the losses we consider biased insertion losses at the input and at collection stages, and balanced internal losses which are known to commute with the unitary matrix. 40 This assumption is not the most general since it considers the optical transformation U without any internal unbalances among the modes. Currently, most of the photonic circuits are designed in such a way that losses are practically unbiased along the evolution and can be factorized to one of the ends of the interferometer. 25 Then, our model of losses can find applications in many scenarios. We model losses as in Fig. 1a with a set of beam splitters placed on the input and output modes. We further consider the presence of unknown (and possibly unstable) phase terms on these modes. The device is thus described by a matrix T = D 1 U D 2 , where D 1 and D 2 are diagonal matrices that account for such losses and additional phases. The matrix elements of the unitary part are expressed as U ij = τ ij e iφ ij , where τ ij and φ ij are the matrix moduli and phases respectively. In what follows, we briefly analyse the two most general algorithms in the literature to reconstruct the matrix U . The first algorithm exploits quantum light 28 and the second one is based on coherent light measurements. 30 If one employs measurements with Fock states at the input of the interferometer and photonnumber-resolving detection at the output, the results will be insensitive to the (unstable) phase terms at the input and at the output. This means that the matrix U and all the matrices U in the form U = F 1 U F 2 , where F 1 and F 2 are a unitary diagonal matrix, are equivalent. Another invariance property of these measurements carried out with Fock states is that the outcomes do not change with respect to the conjugate operation U = U * . Given these equivalence relations, it is not necessary to reconstruct the actual unitary matrix implemented by the interferometer but only a representative element of its class of equivalence. This is commonly defined by choosing a matrix with real-valued elements in the first row and first column (φ 0i = 0 and φ i0 = 0) together with the condition that the element U 11 has positive phase (φ 11 > 0). In Ref., 28 the authors presented an algorithm to reconstruct the value of moduli and phases of the representative unitary matrix elements via single-photon intensity and two-photon measurements, the latter via the visibility V of Hong-Ou-Mandel (HOM) interference. 37 Labelling the input modes of the two photons as h, k, and the output ones as i, j, the visibility is defined as where (P hk ij ) I is the probability to find the two (indistinguishable) photons in output modes (i, j) when they are injected respectively in input modes (h, k), while (P hk ij ) D correspond to the same quantity with distinguishable particles. Note that the visibility V does not change in the presence of mode-dependent losses in both preparation and measurement stages, and thus its estimation gives direct access to the matrix elements. By measuring these quantities it is possible to define a system of equations and to retrieve an analytical solution to the problem, as shown in Refs. 28,29 One of the main constraints of such approach is the requirement of measurement with quantum light. Recently, a further method that exploits quantum light probes, such as two-mode squeezed states, and single-photon counting combined with heterodyne measurements, provided further refinement for reconstructing the U matrix in the presence of internal unbalanced losses. 41 The other method proposed in Ref., 30 based on classical intensity measurements, requires phase stability within the measurement time to estimate matrix phases, since these values are extracted from first-order correlation functions. This means that in such a measurement scheme the equivalence between U and U does not hold. Additionally, the correct estimation of matrix element moduli, which in this approach is performed independently from the phase estimations, is spoiled by the presence of output mode-dependent losses.

Algorithm for reconstruction of linear interferometers
In this section, we propose an alternative black-box methodology based on coherent light probes. This approach permits to estimate matrix elements moduli and losses also in the presence of loss imbalance among the modes and retrieving quantities having the same properties of V without requiring phase stability at the inputs and outputs of the interferometer. Reconstruction of moduli and losses. We start our investigation with the estimation of the matrix elements moduli τ ij . The τ 2 ij coefficients represent the probability to find a single photon in output mode j given an input mode i or, alternatively, the fraction of classical field intensity in the mode j. This observation allows to define a probability matrix P ij = τ 2 ij . The main task is to estimate such a matrix when mode-dependent losses are present in the preparation and collection stages. Let us define M as the matrix obtained from single-mode intensity measurements, estimated experimentally via single-photon input states or by injecting laser light in a single mode. Following the model presented in Fig. 1a, the actual measured matrix M can be expressed as follows where D 1 and D 2 are diagonal matrices that describe the unbalanced losses in input and output modes. The task then requires reconstructing the probability matrix P and the input and output losses matrices D 1 and D 2 starting from the measured matrix M . To this end, we now introduce below two approaches that can be used to estimate D 1 and D 2 up to a multiplicative factor and, consequently, matrix P under very few assumptions on the measured matrix M .
(1) Sinkhorn's decomposition-based algorithm. -A first approach to reconstruct matrices P , D 1 and D 2 is obtained starting from the observation that the probability matrix P is a doubly stochastic matrix, i.e. the matrix has non-negative entries and the sum of each row and each column is equal to one. It is thus possible to apply Sinkhorn's theorem and the matrix scaling algorithm on this system. 42,43 Indeed, a matrix with non-negative elements such as M admits a Sinkhorn's decomposition if it is diagonally equivalent to a doubly stochastic matrix, i.e. can be written in the form D 1 P D 2 where D 1 and D 2 are diagonal matrices and P is a doubly stochastic matrix. This decomposition exactly represents the solution to our problem (see Eq. (3)). Sinkhorn's theorem and subsequent extensions 42 guarantee that this solution exists and it is unique. The theorem gives us also an important warning since the algorithm is sensitive to the position of the zero elements of the measured matrix M . This means that an incorrect attribution of zero-valued elements in M , due to experimental errors or limited measurement sensitivity, could make the matrix impossible to decompose with Sinkhorn's theorem.
Finding Sinkhorn's decomposition for a non-negative matrix M is a special case of the matrix scaling problem that has applications in a large variety of fields. In our case, defining X = D −1 1 and Y = D −1 2 , we can write P = XM Y and using the property that P have to be doubly stochastic, we obtain the following system of equations where e is the vector with all the components equal to 1.
In the literature, different algorithms have been proposed to solve Eq. (4). Here, we present a specific choice among the possible algorithms (see Ref. 44 for a review of the possible approaches). The chosen algorithm allows to recovery all the three matrices in Eq. (3). The idea is to rearrange Eq. (4) in terms of vectors x and y which are the diagonal elements of X and Y . Formally, they can be expressed as x = X e and y = Y e. By defining the inversion of a vector as the inversion component-by-component, At this point we apply an iterative algorithm to solve the system of equation (5) and retrieve the two vectors x and y and consequently the two diagonal matrices X and Y . Then, we can recover the probability matrix P = XM Y and the two loss matrices as D 1 = X −1 and D 2 = Y −1 .
(2) Variance-minimization based algorithm. -In the derivation of the previous algorithm to solve Eq. (3) we have implicitly assumed to have measured the field intensity distribution for any combination of outputs j for any input i. In the following method, we define an alternative procedure that can be applied when only a subset of the inputs are available while requiring the capability of reconfiguring the linear optical network.
Let us then consider an interferometer in which it is possible to change the probability matrix P without affecting the input and output losses D 2 and D 1 , and to measure all the output configurations only for a subset of input modes. We first consider the scenario in which the light source, a single photon or a laser beam, is injected only in mode i. In absence of unbalanced losses, the outcome M of such intensity measurements is proportional to the vector P, (the i-th column of the matrix P ) and represents a discrete probability distribution which depends from a set of parameters θ describing the optical circuit settings. By taking into account the presence of unbalanced losses D (D 1 in the general case), the components M j of vector M are where I is the intensity attenuated by the input loss. Since we are injecting light in the same mode for all measurements, this factor can be included in D.
The separation of the probability from the losses in Eq. (6) can be done starting from the following observation. Let us define the quantity S( α, θ) as the weighted sum of the components of M : If the vector of the weights α is proportional to the element-wise inverse of the losses vector D, the quantity S does not vary with the control parameter θ. In fact, when α = β D −1 , where β is a global factor, we have: In general, S( α, θ) changes with the controls parameters θ. Then, the idea is to find the weight vector α such that S( α, θ) is constant when the controls parameters θ vary. This vector can be obtained by minimizing the variance of S( α, θ) with respect to α. The variance ξ 2 ( α) can be estimated on a sufficiently large set of settings of the parameters { θ i } as where n is the number of different settings { θ i }, and consequently represent the size of the measurement outcome vector { M ( θ i )}. This minimization is equivalent to a quadratic optimization problem. To solve such a task, let us call M ij = [ M ( θ i )] j the matrix in which each row contain the output intensities distribution for a particular configuration of the chip then we define the following positive semi-definite matrix Q as Then our problem can be rewritten, in terms of the matrix Q, as The minimization of Eq. (11) has as trivial solution α = 0 (the vector with all null components) and another one that is the eigenvector of Q corresponding to a null eigenvalue. The latter solution corresponds exactly to the element-wise inverse of the losses vector α = D −1 . In the presence of noisy experimental data, the solution is the eigenvector of Q corresponding to the lowest eigenvalue. Alternatively, the non-trivial solution can be found by an ordinary numerical minimization approach of Eq. (11). This can be fulfilled by setting a normalization constraint N · α = 1 for some normalization vector N , since losses can be estimated with this procedure up to a multiplicative factor common to all the modes.
The method can be generalized to the scenario in which one is interested in reconstructing a sub-matrix of P . Then, it is possible to reconstruct even the relative losses between the different measured input modes. Here, we suppose to have a set of output vectors { M i } k for each input k ∈ K and compute the variance function ξ 2 k ( α) and the associated matrix Q (k) . At this point we minimize the sum of all variances with the same constraints of the previous derivation: After the minimization, it is possible to recover the input losses from the value of the sum function associated with each input as follows Reconstruction of the internal phases with classical light. Here we propose a procedure to estimate the complex phases of the matrix elements φ ij . The methods reported in Ref. 28,29 employs the visibility of the Hong-Ou-Mandel (HOM) effect described in Eq. (2) for this task, by sending a two-photon input state whose indistinguishability is tuned during the experiment. In this work, we propose an analogous quantity that can be measured by intensity cross-correlation at the output of the linear network with classical light. The measurement scheme is presented in Fig. 1b. The laser source is split and sent into the network in modes (h, k). The additional phases ϕ 1 and ϕ 2 account for phase instabilities in the optical paths between the sources and the interferometer. We can define the cross-correlation σ hk ij between the output modes (i, j) when the two beams enter from modes (h, k) as where I i and I j are the field intensities in the corresponding output modes, while · is the time average. We define also the self-correlation σ hk ii of the intensity fluctuation as We make the hypotheses that (i) the external phase fluctuations ϕ = ϕ 1 − ϕ 2 have zero time average, and (ii) the input laser intensity is constant. After some calculations, reported in Supplementary Materials S.II, we can define the normalized cross-correlation C hk ij as Note that this quantity only depends on the phase of the matrix elements. Similarly to HOM visibility in two-photon experiments, the set {C hk ij } does not depend on the input and output losses. Additionally, and at variance with HOM visibility, {C hk ij } does not depend also on the moduli {τ i,j } of the matrix elements, thus permitting an independent estimation of the phases. The derivation of Eq. (16) is performed under a specific assumption on the external optical phase fluctuations at the input. More specifically, we require that In general, mechanical and thermal phase fluctuations do not satisfy these conditions. These equations can be satisfied by adding a phase modulator in one of the two input paths. In this scenario, the external phase contribution can be expressed as ϕ = ϕ M + ϕ T , where ϕ M is the modulated phase and ϕ T is the one modified by thermal and mechanical noise. Since the two contributions are uncorrelated, we can write e ıϕ = e ıϕ M e ıϕ T . Controlling the phase modulation such that e ıϕ M = 0 and e 2ıϕ M = 0, for example by adding white noise with appropriate amplitude or by a discrete set of phases, the conditions of Eq. (17) can be satisfied. Complete algorithm. Given the methods defined above, we can summarize the complete procedure to reconstruct the matrix U as follows 1. Perform field intensity measurements in the output of the circuit. Apply the Sinkhorn-based algorithm to retrieve the complete set of moduli for the matrix elements ({τ i,j }) or the variance minimization approach to estimate a specific subset.
2. Perform cross-correlation measurements in pairs of the output of the circuit. Solve the system of equations for the set {C hk ij } to extract the complex phases of the unitary matrix. Note that point 2 has some similarities to the procedure proposed in Refs. 28,29 with HOM visibilities. In particular, this approach could require a further minimization step on a larger set of {C hk ij } with respect to the minimum ensemble needed to solve the system. Nonetheless, it is worth noting that in our case we do not require (i) the use of indistinguishable single-photon states and (ii) the measurements of {C hk ij } give us directly the information about the phases without requiring the knowledge of the matrix moduli. In fact, in our algorithm points 1 and 2 are independent, as for the previous methods with coherent probe light, 30, 33 but having the additional features of permitting the estimation of losses. Furthermore, it can be applied in any scenario with phase instabilities.

Experimental verification in a reconfigurable integrated circuit
We tested the complete algorithm described in the previous section in a 3-mode reconfigurable optical circuit. The waveguides of the device were fabricated in a glass sample via the femtosecond laser-writing technique 45 (See Supplementary materials S.III for more details). The structure is composed of a sequence of two tritters, a circuit that generalizes the beam-splitter over three modes 46,47 (see also Fig. 1c). Between the two tritters the presence of three resistive heaters permits to change the unitary transformation performed by the circuit via the thermo-optic effect. The measurements were performed with a continuous wave laser at the wavelength of λ = 785nm. The laser is routed in different input modes via a fiber switch connected to a fiber array that injects light into the input port of the interferometer. The field intensity distribution in the 3 output ports was recorded via a CCD camera.
Our first experimental test regards the two algorithms described above, to retrieve the losses vectors and moduli of the matrix elements. In Fig 2 a) and b) we report the results for the application of Sinkhorn's decomposition method. In particular, in panel a) we report the measured field intensity M for a particular configuration of the chip, normalized to the column total intensity. Here we inserted on-purpose additional losses by placing an attenuation filter on the output mode 0, to test the performances of the approach. In Fig 2 b) we report the probability matrix P after the applications of the Sinkhorn algorithm. We measured the moduli of a set of N different transformations U , each of them in two loss conditions, namely by inserting and removing the attenuation filter in mode 0. Defining the fidelity as F = (1/3 2 i,j=0 P ij P ij ) 2 , the average classical Similarity between the two reconstructed distributions in the two losses configurations for a given U via the Sinkhorn method isF = 0.9999 ± 0.0001 (F min = 0.9997). The average was estimated on the set of the N pairs of matrices. This confirms that the method works properly in different mode-dependent loss configurations, and retrieves as expected always the same moduli |τ i,j | 2 associated with the transformation U .
We then move to test the second algorithm based on variance minimization. In Fig. 2 c) we report the measured field intensity in the three outputs of the interferometer for different dissipated powers in one heater, thus corresponding to different evolutions U , when the signal is injected in input 0. We also report the sum (red curve) of the three intensities for each tested configuration for U . We observe that such a sum is not constant, as one would expect if the output losses were balanced. In d) we report the same curves after applying the variance minimization algorithm, showing that Normalized correlation outputs (0,1) Normalized correlation outputs (0,2) Normalized correlation outputs (1,2) a) b) c) Figure 3 Cross-correlation measurements. In this figure we report the measurement of the normalized crosscorrelations C hk ij for different pairs of outputs, entering from (h, k) = (0, 1). In particular in a) we measure the pair (0, 1), in b) the pair (0, 2) and in c) the pair (1, 2). In red we report the experimental correlations for different configurations of the dissipated electrical power in the heater R 1 . In blue the predictions according to the results of a white-box fit that makes use of the structure of the interferometer and the previous measurements of the matrix moduli. The two independent estimations are in good agreement within one standard deviation of the experimental error. the application of the algorithm makes the sum constant with respect to changes in the internal transformation. We repeat the same procedure of the previous paragraph by tuning the output mode-dependent losses, showing the capability of retrieving the correct moduli values. The average fidelity, among the same set of N internal operations, between the reconstructed distributions in the two different conditions of losses after the application of the algorithm isF = 0.9996 ± 0.0002 (F min = 0.9990). As a further confirmation, we compared the distributions, corresponding to the same optical circuit settings, retrieved by the application of the two algorithms. The mean fidelity between the reconstructed matrix with the two methods isF = 0.999 ± 0.001 (F min = 0.992), thus confirming the effectiveness of both algorithms.
Finally, we tested the measurement of the cross-correlations defined in Eq. (16) for the phase reconstruction. Since the phase fluctuations of the fibers do not fulfil the conditions of Eq. (17), we placed a liquid crystal in one arm before the first input of the photonic chip. For the phase modulation, we used a discrete set of phases instead of a continuous one. In particular we employed the set {0, 2π/3, 4π/3}. After recording the temporal fluctuations of the output field intensities, we calculated the normalized cross-correlations for various configurations of the interferometer (red dot in Fig. 3). These measurements are compared with the predictions made by an independent characterization of the device via a white-box algorithm, used as a reference to test our approach. By this alternative method, which exploits the structure of the two tritters, the moduli and the phases of the matrix are retrieved directly from the field intensity distributions of the previous paragraph. It follows that, in this white-box approach used as a reference, the cosines of the internal phases are the result of a numerical optimization between the parameters of the optical circuit, which is designed according to the structure in Fig. 1c, and the measurements of the P distributions. The blue curve in Fig. 3) represents the prediction of such an optimization. We observe that the direct measurement of cross correlation via the proposed approach with classical light, and the same quantities retrieved via the white-box method, are compatible with the experimental error, with a normalized χ 2 = 1.09. As a final comparison, we compute the fidelity between the second column of the matrix retrieved via our black-box algorithm (bb) U bb i1 and the column of the white-box approach (wb) U wb i1 . We choose the second column of the matrix since it presents non-trivial phases, namely φ i1 = 0 for i > 0. Indeed, the first row and column are imposed to be real vectors because of the equivalence between U and U = F 1 U F 2 with F 1 and F 2 unitary diagonal matrices. The fidelity defined here as i.e the overlap between pure single-photon states described by the second column of U , has been estimated for each configuration of the dissipated power in R 1 reported in Fig. 3. The average fidelity isF = 0.999 ± 0.001 (F min = 0.994).

Discussion
In this work, we presented new algorithms to characterize the operation of multi-mode linear optical circuits. In particular, we have shown the possibility to reconstruct the moduli of the unitary matrix element by field intensity measurements, in the presence of unbalanced losses at the input and output ports. This is in contrast to the previous black-box algorithms 28,29 that can reconstruct the moduli of the matrix only after phases measurements via HOM visibility and by imposing the constraint to have a unitary matrix. In addition, our procedure can provide directly the differential losses among the waveguides at the input and output. We also proposed a new method to characterize the internal matrix phases based on intensity correlations of coherent light beams at the outputs of the linear network. These measurement methods do not require knowledge of the matrix moduli and of input/output losses. These approaches enable the possibility to characterize separately the moduli and the phases of the unitary matrix implemented by the optical network, with a reliable and independent approach. Furthermore, all presented algorithms are polynomial in the number of modes of the interferometer, both for the execution time and for the number of measurements required, sharing the same scaling of previous state-of-the-art algorithms. This property makes our algorithms suitable for adaption to large-scale interferometers. An open issue regards the inclusion of unbalanced internal losses in the model. Some recent works have proposed some solutions for taking into account such sources of noise in the U reconstruction. 41 We foresee as future perspective of our work to adopt the same strategy of independent estimations of moduli and phases to detect the effects of such imperfections. We provided experimental proof of the effectiveness of the algorithm by verifying its performance on a 3-mode reconfigurable integrated optical circuit. In this analysis, we compared the results with the predictions of an independent reconstruction that exploits the knowledge of the internal structure of the circuit, which may be in general unknown or too complex to model when taking into account noise processes.
Our findings pave the way to the successful adoption of such an effective black-box methodology to large-scale optical networks, that nowadays are approaching a high number of optical modes and components, 19,22,23 with applications ranging from quantum communication and sensing to quantum computation and simulation via photonic systems.

S.I Quantum correlations
In this section, we recall the proof of Eq. (2) of the main text. First, we calculate the states at the outputs i and j when two photons, respectively in the state |φ e |ψ , are injected in the input modes h and k. We first choose an orthonormal basis to express the states on the two input modes: Then, the two-photon states is expanded as: The state after the unitary evolution is: Since we are interested in the probability to measure two photons in two different output ports i = j, we calculate the probability that the output state is in the non-collisional term In this expression, we observe that the indistinguishable photons scenario is obtained for | φ|ψ | 2 = 1, while the distinguishable particle case corresponds to | φ|ψ | 2 = 0. Using the definition of the HOM visibility, we find Eq. (2) of the main text. When residual distinguishability is present the measured visibility follows the following equation: If this effect is not properly taken into account, this can lead to errors in unitary reconstruction algorithms based on such quantity. Finally, we observe that is possible to obtain a quantity analogous to the normalized correlation of Eq. (S15), starting from visibility measurements. This procedure requires the measurements of the bunching probabilities, that is, the probability that two photons exit from the same output mode. Indeed, using this information it is possible to can reconstruct the following quantity

S.II Classical correlations
In this section, we show that second-order classical correlation measurements are described by an expression which is equivalent to two-photon Hong-Ou-Mandel visibilities. Starting from Figure S1, we consider two laser beams at the input of the interferometer. The input fields are described as:Ẽ where φ 1 (t) and φ 2 (t) are the phases introduced via propagation in optical fibers, and by all the other possible optical delays in the apparatus. After propagation through the interferometer, the electric fields on the output modes read: S1. Conceptual scheme of the apparatus for the reconstruction of matrix phases via second-order correlations with classical light.
Then, the output intensity of light is measured via two photodiodes: If we suppose that the intensity of the input lasers I 1 and I 2 are constant and that Ẽ 1Ẽ * 2 = 0 (that in our scenario is equivalent to e ı(ϕ1−ϕ1 = 0), we can calculate the residual as: At this point we can define the cross-correlation σ hk ij between the output modes (i, j) when the two beams enter from modes (h, k), and the self-correlation σ hk ii of the intensity fluctuation a: σ hk ij = (I i − I i )(I j − I j ) σ hk ii = (I i − I i ) 2 (S13) Let us know define a parameter γ = Ẽ 1Ẽ * 2Ẽ * 1Ẽ2 /(I 1 I 2 ), related to the first order correlation functions of the two beams that it is in turn related to visibility of the interference fringes. In addition, we assume that the fields satisfy the following hypothesis (Ẽ 1Ẽ * 2 ) 2 = 0, which is equivalent to the condition e ı2(ϕ1−ϕ1 = 0. Under these assumptions, the cross-correlations can be calculated as: Finally, by defining the normalized cross-correlation C hk ij we obtain Eq. (19) of the main text: We note that these quantities do not depend on the visibility γ nor from the moduli of the unitary matrix. We obtained something equivalent in the quantum correlation section but with the need for single-photon number resolving detection to evaluate the bunching probability.

S.III Chip fabrication
The photonic chip was fabricated in a borosilicate glass substrate (EagleXG, from Corning Inc., USA) by femtosecond laser direct writing. In detail, irradiation with focused femtosecond laser pulses causes a permanent refractive index increase in the glass, localized in the focal region (see Fig. S2); translation of the substrate with respect to the laser beam allows for the direct inscription of waveguides along the desired paths. A Yb-based cavity-dumped oscillator was adopted as a laser source, emitting ultrashort pulses with 1030 nm wavelength and 300 fs duration at the repetition rate of 1 MHz. For waveguide inscription, laser pulses of 250 nJ energy were focused 30 µm below the substrate surface by a 0.6 NA microscope objective, while the substrate was translated at a constant speed of 30 mm/s. These parameters enabled the fabrication of single-mode waveguides operating at the wavelength of 785 nm, with a 1/e 2 mode size of 7.2µm × 8.4µm and propagation losses (for vertical polarization) lower than 0.8 dB/cm. The depth of the waveguides was chosen to allow the control of the optical phase via thermo-optic phase shifters. The thermal shifters consist of metallic resistors, manufactured by the same femtosecond laser through the ablation of a thin and uniform gold layer (∼ 60-nm thickness) deposited on the chip surface. Each resistor has a width of 100 µm and a length between 5 and 7 mm (parallel to the waveguide direction), which gives resistance values in the range 60-100 Ω. Electrical connections are provided via standard pins, directly glued on the electrical circuit terminations.
FIG. S2. Typical cross-section of an optical waveguide inscribed with a single irradiation step of femtosecond laser in the borosilicate glass substrate, imaged with an optical microscope. The parameters of the laser beam and the focusing conditions are described in the text. The femtosecond laser is impinging from the top in the figure, and the scale bar corresponds to 10 µm. The refractive index modification assumes a drop-like shape with complex features; however, the waveguide supports a single optical mode at the wavelength of interest.