Reconfigurable continuously-coupled 3D photonic circuit for Boson Sampling experiments

Boson Sampling is a computational paradigm representing one of the most viable and pursued approaches to demonstrate the regime of quantum advantage. Recent results have shown significant technological leaps in single-photon generation and detection, leading to progressively larger instances of Boson Sampling experiments in different photonic systems. However, a crucial requirement for a fully-fledged platform solving this problem is the capability of implementing large-scale interferometers, that must simultaneously exhibit low losses, high degree of reconfigurability and the realization of arbitrary transformations. In this work, we move a step forward in this direction by demonstrating the adoption of a compact and reconfigurable 3D-integrated platform for photonic Boson Sampling. We perform 3- and 4-photon experiments by using such platform, showing the possibility of programming the circuit to implement a large number of unitary transformations. These results show that such compact and highly-reconfigurable layout can be scaled up to experiments with larger number of photons and modes, and can provide a viable direction for hybrid computing with photonic processors.

Since the original proposal of a computational paradigm based on the rules of quantum mechanics [1], large research efforts have been devoted to identifying the optimal approach to implement a universal quantum computer [2,3].Besides the great promises provided by such computational paradigm, the implementation of a largescale universal quantum device capable to outperform or, at least, to be comparable with a classical computer is still a challenging task.In view of the very recent advances in quantum technologies, approaching the noisy intermediate-scale quantum (NISQ) era, the current target is to reach a fundamental milestone named quantum computational advantage.The goal is to achieve, unambiguously and possibly with different platforms and approaches, the scenario where a quantum device is capable of solving a specific task faster than any classical counterpart.
Among the several quantum algorithms that provide a computational speed-up, two strategies emerge as the most suitable for their experimental realizations [4].The first one requires to sample the output states produced by a random quantum circuit.This paradigm has been recently implemented in quantum processors based on superconducting qubits [5,6].The second approach is based on a different albeit related task, named Boson Sampling (BS) [7,8], which requires sampling from the distribution of non-interacting bosons scattered by a random unitary transformation.Here, the hardness of simulating such a bosonic system is strictly connected to quantum interference effects due only to particle indistinguishability [9].Classical simulation of Boson Sampling, even approximately, is computationally-hard since it requires calculation of permanents of complex-entried matrices, a #P-hard problem.A natural way to reproduce this dynamics and sampling effectively from the output distribution of such process can be obtained via a photonic quantum processor.Reaching the quantum advantage regime with this approach requires the generation of a set of n highly-indistinguishable single photons, which evolve via quantum interference in a low-loss linear optical network with m∼n 2 modes capable of implementing a random transformation according to the Haar measure (see Fig. 1a).Samples are then collected by direct measurement of the output modes.In this direction, several proof-of-principle implementations have been reported exploiting photonic platforms [10][11][12][13][14][15][16][17][18][19], including one of the latest experiments with the detection of 14 photons after propagation in an interferometer with 60 ports [20].
Starting from these results, research efforts have been dedicated to technological advances enabling to enlarge the dimensionality of photonic processors, and to theoretical investigations [21,22] aimed at a precise definition of the limits of a classical simulation of Boson Sampling with experimental imperfections.In parallel, the Boson Sampling paradigm has triggered the definition of a set of variants of the original formulation with the aim of improving the efficiency of the corresponding photonic platform while preserving the classical computational complexity of the task.Examples include (but are not limited to) Scattershot Boson Sampling (SBS) [23] and Gaussian Boson Sampling (GBS) [24].SBS exploits probabilistic parametric down-conversion sources placed at each input mode of the interferometer in a heralded configuration (see Fig. 1b), leading to an exponential growth in the samples acquisition [18,25,26].In the GBS paradigm [26][27][28], Fock states are replaced by single-mode squeezed vacuum states (see Fig. 1c).By using the GBS approach, recent experiments have reported the achievement of quantum advantage with a photonic platform [29,30].Furthermore, GBS has been recently pointed out to have potential application for hybrid classical-quantum computing, due to the connection with other problems including graph theory [31,32] or simulation of molecular vibronic spectra [33,34].In all previous experiments reporting Boson Sampling instances, including its variants, different platforms have been employed to realize linear optical transformations.However, all the requirements for a fully-developed processor, namely low-loss, reconfigurability, and the possibility to implement random transformations, are currently not fulfilled simultaneously in a single system.Indeed, low-loss systems are necessary to avoid spoiling the complexity of the process [35,36], and full-reconfigurability is a crucial requirement on two main aspects.On one side, such a feature is essential to benchmark the effective Haar-randomness of the platform, which is at the basis of complexity conjectures [7,37].On The corresponding computational problems require sampling from the output distribution using different input quantum states of light, such as Fock states in BS, two-mode squeezed vacuum states in SBS and single-mode squeezed vacuum states in GBS.The common element among the schemes is the optical random circuit, described by the unitary evolution U. d) The most widely adopted decomposition of the operator U is via a network of beamsplitters, with splitting ratios Sij, and phase-shifts φi (left); an alternative implementation exploits continuous-coupling by evanescent waves among waveguides (right) depending on the coupling coefficients cij and the propagation constants ki, where both may vary along the direction z. e) Overview of the reconfigurable 3D integrated photonic chip, realized through the femtosecond laser writing technique.The device is composed by 32 optical-modes arranged in a triangular lattice, as showed in the inset reporting the transverse section of the sample.In red we have highlighted the input modes employed in the 3-and 4-photon experiment.The transformation U is controlled by the 16 resistors fabricated on top of the glass sample.The second inset shows the top view of the electrical circuits that controls the currents Ii applied to the resistors.
the other side, reconfigurability is needed to achieve programmable processors for applications beyond the quantum advantage demonstration [28].With current-up-to-date experiments, low-loss platforms lack an active reconfiguration capability to change the transformation implemented by the optical circuit [20,29,30] and their capability of Haar-random operations has not been demonstrated yet.Conversely, integrated photonic architectures based on universal decompositions [38,39] permit full reconfiguration capabilities, but still require further technological improvements also in terms of loss-reduction to scale up the number of modes.
Here, we perform a step forward towards developing a photonic platform encompassing all the aforementioned characteristics.We report the realization of a photonic reconfigurable integrated device with a compact structure.Such device provides significant advantages in terms of losses and number of modes.In addition, it possesses a high degree of reconfigurability enabling the implementation of a large number of transformations, in contrast to previous works that implement static integrated circuits via the same technology [19,40].This architecture takes advantage from the 3D-capability of the femtosecond laser-writing technique of waveguides in glass [41,42].The optical circuit is realized via continuous-coupling of 32 waveguides arranged in a triangular lattice.The device is then controlled by changing the currents applied to the 16 heaters fabricated on top of the interferometer.We show the reconfiguration capabilities of the device by providing a thorough analysis on the set of unitary transformations which can be reached by the system.Then, we benchmark our platform in the Boson Sampling scenario by performing and validating 3-photon and 4-photon experiments with several different configurations of the optical transformation.

Reconfigurable integrated 3D photonic chip
Integrated photonics has demonstrated significant advances in realizing complex optical circuits for quantum information processing [42].Nowadays integrated optical circuits are the most promising photonic platforms to implement large-scale interferometer with high-level of reconfigurability.The unitary operations that describe a given optical circuit are often decomposed in elementary optical units.The scheme by Reck et al. [38] decomposes an optical circuit with m optical modes in the product of m(m−1)/2 two-mode beamsplitters and single-mode phase-shifters.A more recent algorithm by Clements et al. [39] optimizes the arrangement of the optical elements to minimize the sensitivity to fabrication imperfections and to photon losses.In this framework, any unitary transformation is determined by the set of phases φ i and splitting ratios S ij in the circuit (see Fig. 1d).In reconfigurable integrated circuit, heaters placed in correspondence to the waveguide enable to change locally the refractive index of the material thus inducing a change in the relative phases among the optical modes.Despite the aforementioned schemes have been adopted in several experiments [12,28,[43][44][45][46], scaling the circuit to large number of modes is a challenging task in terms of size of the device, related to the amount of losses and to the number of required optical elements.In Fig. 1d we depict the approach that we employ in this paper, which is different with respect to the traditional decomposition of unitaries in beamsplitters and phase-shifters.Such a scheme exploits the continuous coupling of the radiation in waveguides arrays.In this framework it is possible to associate a Hamiltonian H(z) to the system given the coupling coefficients c ij (z) between neighbouring modes and the propagation constants of the modes k i (z).Coupling coefficients and propagation constants may change along the propagation coordinate z.
The unitary transformation of the circuit can be obtained by integrating H(z) for the whole length of the interaction region, which maps the time evolution of the system.In particular, in our case we employ a three-dimensional array, in which waveguides are arranged according to a triangular-lattice cross-section.Such architecture offers practical advantages in terms of compactness, losses and circuit length in comparison to other approaches.Exemplarily, it is striking to note that in this work we have achieved reconfigurable 32-modes random unitaries in a 7.5-cm long optical chip, while a reconfigurable discretely-coupled interferometer with the same number of modes and produced with the same technology would have required a device length of about 30 cm.Such aspects are widely discussed in Supplementary Note 1.In detail, our integrated device, depicted in Fig. 1e, comprises 32 continuously-coupled single-mode waveguides in a 8×4 arrangement, fabricated by femtosecond laser direct writing in a borosilicate glass substrate [41,47].The waveguide positions are randomly modulated along the propagation coordinate z with respect to the positions of a regular triangular lattice.These modulations introduce randomness in the c ij coefficients and thus in the circuit transformation, which otherwise would be highly symmetric and not suitable for a Boson Sampling experiment.To add circuit reconfigurability, 16 resistive heaters have been patterned, by femtosecond laser ablation, on a gold film deposited on the substrate surface.The resistors are equally distributed on the two sides of the interaction region.An external power supply controls the currents applied to the resistors.Dissipated power in this process induces thermal gradients in the substrate and thus locally changes the refractive index of the waveguides [48,49].This modulates the propagation constants k i of the waveguides by the thermo-optic effect, thus allowing to change dynamically the transformation U implemented in the integrated device.
Hence, we demonstrate that the thermo-optic phase-shifting technology can also be applied to three-dimensional, continuously-coupled waveguide devices.Additional details on the circuit geometry and on the fabrication process are provided in the Methods section.

Experimental platform
Let us now illustrate the components of the experimental platform employed to benchmark the photonic device and to collect the 3-and 4-photon samples (see Fig. 2a).We exploit one-and two-pair emission in a type-II spontaneous parametric down-conversion source composed of a beta-barium borate crystal operating at 785 nm.The first stage of the apparatus includes all the optical components to generate the state resource to perform sampling with either indistinguishable and distinguishable photons.Photons spectra are filtered through a 3nm band-pass filter.Then, photons are split in four different spatial modes according to the polarization via half-wave plates and polarizing beamsplitters, and coupled into single-mode fibers.Photons are controlled in polarization and in time-of-arrival by polarization controllers and delay lines respectively, in order to tune their degree of indistinguishability.Then, they are injected into the reconfigurable integrated chip via an array of 6 single-mode fibers that has been aligned and glued to the device.A fan-in waveguide section leads the photons to selected inputs of the waveguide array (see Fig. 2b).After the evolution in the integrated device, a fan-out waveguide section leads the photons to a 8×4 rectangular multimode fiber array that matches the fan-out geometry.It is worth noting that the 2D fiber array further helps the compactness of the device by greatly reducing the length of the fan-out section (see Supplementary Note 1 and Supplementary Figures 5-6).The detection stage includes 32 single-photon avalanche photo-diodes.We have developed a custom software that simultaneously controls the delay lines, the power supply and the 32-channel time-to-digital converter module to record two-and four-fold coincidences.This implies a full control over the unitary transformation implemented in the cir- cuit, the switching between indistinguishable and distinguishable photons, and the recording and processing of the data samples.

Unitary matrix sampling and reconstruction
Control over the unitary transformations implemented in the chip is performed by tuning the currents in the resistors.To verify the classes of matrices U that can be implemented by the device due to its reconfigurability, we have reconstructed a large noumber of different evolutions, each corresponding to a different setting for the currents in the resistors.Thus, a crucial ingredient was the adoption of a fast and efficient reconstruction algorithm.Additionally, an efficient and fast reconstruction of the unitary transformations is a fundamental step also for benchmarking and validating the 3-and 4-photon experiments described in the next section.We made use of an adapted version of the method reported in [50] (see Supplementary Note 2 and Supplementary Figure 7).The latter envisages the measurements of two-photon Hong-Ou-Mandel (HOM) dips resulting from pairs of photons injected in different combinations of input ports.In our case we restricted the measurements to two of the possible input pair combinations for the reconstruction of 3×32 submatrix and to only three pairs in the case of 4×32 sub-matrix.From each input pair we analysed 496 HOM dips, namely all the possible non-redundant pairs obtained by the combination of the 32 output ports.From these measurements we have extracted the information about the moduli ρ ij and the phases θ ij of the sub-matrix elements expressed as U ij =ρ ij e iθij .In Fig. 3a we have reported an example regarding the 96 squared moduli and phases of one of the 15 different 3×32 sub-matrices reconstructed in this work.Our next step was to prove that the random unitaries, sampled by changing the currents configuration in the circuit, were drawn from a distribution as close as possible to the Haar measure.This requirement is fundamental to ensure the hardness of BS.To this aim we compared the distributions of phases and moduli of the 15 sub-matrix measured in the experiment with the one retrieved from likewise Haar-random extracted unitary matrices (Fig. 3b-c).In both cases we have obtained a good agreement between the two distributions.The slight deviations from the theoretical histograms can be explained by experimental artifacts related to losses, imperfections in the apparatus and to the algorithm for the reconstruction of the matrix.We discuss these effects in Supplementary Note 3 and Supplementary Figures 8-10.As a final benchmark of the device, we measured 200 columns of different unitaries and calculated the similarities between the distributions given by the squared moduli of each column.The measurements were performed by sending one photon in the device and measuring it at the output in coincidence with his correlated one in a two-photon experiment.The unitaries have been generated by a uniform sampling of the electrical power dissipated in the resistors.Also in this analysis we find a good agreement with expectations, signified by the overlap with the histogram of the similarities calculated from the columns of Haar-random matrices shown in Fig. 3d.The latter result represents one of the first investigations on the level of randomness that can be reached in this continuously-coupled waveguide architecture by changing only the propagation constants via the thermo-optic effect.The similarity to the Haar-random distribution could be improved by engineering the sampling strategy of the dissipated powers, which here have been extracted from a uniform distribution.Note that, in the discrete-decomposition schemes (Reck, Clements) there exist algorithms to set the optical circuit to sample from the Haar distribution.However, the phase shift values and beamsplitter reflectivitie do not display trivial distributions [51,52] which in turn require complex settings of the external control circuit.More precisely, by increasing the dimension of the matrix, the parameter distributions tend to be more and more peaked.For example, the uniform sampling of the dissipated electrical powers employed in this work is far from the correct sampling to generate random Haar matrices in discrete optical circuits.An exhaustive and conclusive answer regarding the possibility to extract matrices from a distribution closer to the Haar measure with a continuously-coupled waveguide architecture needs further studies both from a theoretical and experimental point of view.

Experimental Boson Sampling in a 3D reconfigurable circuit
After characterization of the device, we have then performed 3and 4-photon experiments in the Boson Sampling framework with our integrated system.The 4-photon state from double-pair emission generated by the source can be written as ρ in ∼ α|1111 1111|+ β|2002 2002|+γ|0220 0220| (see Supplementary Note 4 and Supplementary Figure 11), by expressing the density matrix in the occupation number of the four modes and neglecting the higher order of multi-pair emission.To inject a 3-photon |111 input state, one of the four output modes of the source is directly measured and thus acts as a trigger.To this end, we have discarded one output mode of the chip, due to the requirement of using one detector for the trigger photon.Four-fold coincidences between the trigger photon and three output modes are then recorded, providing the output samples.In the 4-photon experiment, we have sampled from the entire ρ in by directly connecting the four output modes of the source to the integrated device (see Fig. 2).In this case, we have sampled from all the output ports of the device.For all reported experiments, our measurements are restricted to the collision-free events.Such choice does not affect significantly the outcomes of the experiment since the configurations with more than one photon per mode display very low probabilities in the regime where the number of photons is much smaller compared to the number of the optical modes [53].
Let us now illustrate the analysis of the experimental samples collected in a 3-and 4-photon BS routines.In this context, the problem of data validation is pivotal to assess the correctness of the sampling process, especially in the regime in which it is not possible to reproduce the output of the experiment with classical resources.In the past years several tests have been developed to rule out classical models, such as the uniform and distinguishable particle samplers, that could reproduce some features of the BS output distribution [54][55][56][57][58][59][60][61][62][63].In Fig. 4 a-b we report one instance of the validation of a 3-and 4-photon BS against the uniform distribution test [54], assuming the input states described above for the two scenarios.We also employ the same data to validate the experiment against the distinguishable particles hypothesis [55] (see Fig. 4 c-d, Supplementary Note 5 and Supplementary Figures 12-15).In both tests each event collected in the experiment increases or decreases a counter, W for the case of uniform sampler and C for the distinguishable particle sampler, according to a likelihood ratio test.Positive slopes are the signatures of the successful validation of the data.These kinds of hypothesis tests require a good modelling of the system, including the knowledge of the input state and of the unitary transformation applied to the state.The latter has been reconstructed through two-photon measurements and exploiting the reconstruction algorithm discussed above (see also Supplementary Note 2).We have performed 10 different 3-photon BS experiments, and 3 different 4-photon with likewise configurations of the optical circuit.All the data were successfully validated against the two hypotheses.In Fig. 4 e-f we underline the sensitivity of these validation tests to the reconstruction of the matrix representing the optical circuit.The histograms report the distribution of the slopes normalized to the one of distinguishable particles when the unitary is chosen randomly and does not coincide with the actual transformation performed in the circuit.We note that, in absence of correspondence between the unitary transformation related to the data and the one related to the likelihood ratio computations, these tests assign the data to the negative hypothesis independently from the particle statistics.In the same figure we report the experimental slopes describing the set of 10 different 3-photon BS validated in this work using the unitary matrices reconstructed via our algorithm.The experimental points (stars in the Fig. 4e-f) are distant more than 3 standard deviations with respect to the average of the histograms, representing validations with random unitaries that do not correspond to the circuit from which the samples were generated as described above.Such additional result reinforce the successful validation of the performed experiments, and benchmark the reconstruction accuracy of the optical circuit showing a high degree of control of the platform.

DISCUSSION
In this article we have reported on the implementation and benchmarking of an integrated platform, with a compact 3D layout based on continuous waveguide coupling which includes a set of heaters to enable a high degree of reconfigurability.The employed architecture, achievable by exploiting the unique capabilities of the femtosecond laser micromachining technique, can be scaled up to larger number of modes.Indeed, we have successfully shown that our 32mode device can implement a relevant portion of transformations according to the Haar measure, a fundamental requirement to fulfil the randomness hypothesis at the basis of the complexity of Boson Sampling.We have then implemented and validated 3-and 4-photon experiments, showing the viability of using such platform for future large-scale Boson Sampling instances.In fact, the devices fabricated with femtosecond laser technique can be interfaced with other types of single-photon sources, such as deterministic quantum dot emitters [64].Deterministic sources are the most promising for scaling up the number of indistinguishable photons in a genuine Fock state [65,66].The same waveguide fabrication technology is able to realize integrated parametric sources [67] that could be included in our device to mitigate the current losses in coupling photons to single-mode fibers and at fiber-waveguide interfaces.Furthermore, the scheme with arrays of integrated sources could be feasible for realizing Scattershot and Gaussian Boson Sampling variants in a fully integrated platform [26,28].The capability of devising a fullyreprogrammable, large scale and compact photonic processor is fundamental also to envisage future applications beyond the original scope of the Boson Sampling.To this aim, future studies shall focus on methods that allow one to control and program the transformations that the device can implement.This means finding a model that links the parameters of the unitary transformations to the dissipated electrical powers in the heaters, which is a problem still not addressed in the literature for reconfigurable continuously-coupled waveguide circuits.We believe that the use of a black-box approach, via neural networks and optimization algorithms [68,69], is possible and promising.Further improvements in the reconfigurability could be brought from modifications of the present device architecture, e.g. by increasing the number of heaters and by changing their arrangement with respect to the waveguide positions.
Given the aforementioned potentiality and advantages reported in this work such a platform is expected to be at the basis of recent proposals of hybrid computing architectures [31-34, 70, 71].METHODS 3D photonic circuit.The circuit consists of a triangular lattice with 32 continuously-coupled waveguides characterized by an average pitch of 11µm (see Fig. 1e).Indeed, each waveguide is shifted from its standard position of a quantity randomly chosen between 0 and 2µm along a random direction.Such shifts are varied with continuity along the whole coupling region, which is overall 36 mm long.In order to allow the coupling of single photons into the circuit, 6 central waveguides of the array extend through a fan-in region, rearranged in a single row with spacing of 127µm.At the output side, a fan-out expands all the waveguides in a 8×4 rectangular lattice with a pitch of 250µm.Photons are coupled to the input of the circuit with a linear single-mode fiber array, while are collected at the output by a rectangular multimode fiber array (see Fig. 2b).Input and output fiber arrays match the geometry of the fan-in and fan-out regions, respectively.Total insertion losses are 3.5±0.1dB,depending on the input waveguide considered.
Fabrication process.The 3D photonic circuit was fabricated by femtosecond laser writing in a boro-aluminosilicate glass substrate (EAGLE XG, Corning) extending on an area of 75×12mm 2 .The laser source (PHAROS, Light Conversion), operating at a wavelength of 1030nm, was configured to produce pulses with duration of 170fs and energy equal to 290nJ at a repetition rate of 1MHz.The laser was focused with a 20× (NA = 0.5) water-immersion objective, while the substrate was translated at 20mms −1 for six consecutive scans.To ensure efficient reconfigurability, the circuit was inscribed at 30µm from the surface.After a thermal annealing step [72], single-mode waveguides operating at 785nm with a 1/e 2 mode diameter of 4.5µm were obtained.Heaters were fabricated by depositing gold on the chip surface and patterning the electrical circuit with the process reported in [73].A total number of 16 resistors (length 3mm, resistance 70±13Ω) were arranged in two parallel rows at the sides of the coupling region.To guarantee proper heat dissipation, the device was mounted on an aluminium heat sink.In a n-photon Boson Sampling experiment, the collection rate of significant output events scales unfavourably with the n-th power of the transmission of the optical setup.Hence, optical losses have to be reduced as much as possible, to enable experiments with larger and larger number of identical photons.In bulk optics, it may be natural to attribute a given amount of loss to the individual components, while propagation in free space is reasonably considered as lossless.On the contrary, in integrated optics some photon loss is intrinsic in waveguide propagation, especially in the curved parts, while components such as directional couplers may not introduce, because of their operation, specific additional loss.Therefore, in the design process of an integrated optical circuit, minimizing the propagation length is of paramount importance.
We discuss in the following different layouts for an integrated optics interferometer.We will evaluate in particular the footprint required to achieve unitary transformations useful for a Boson Sampling experiment, and we will estimate how the physical size of the circuit scales with the number of modes.We will thus bring arguments in favour of our choice for a three-dimensional array of continuously-coupled waveguides, arranged on a triangular lattice.

Planar Clements interferometer
Let us first determine the footprint of a planar circuit that performs an arbitrary unitary operation U on m optical modes, by means of a network of discrete phase shifters and beam splitters.We adopt the scheme proposed by Clements et al. [39], which is indeed more compact than the earlier scheme by Reck et al. [38], and which has an optical depth of m beam splitters.The optical depth is defined as the maximum number of such elements that a photon needs to pass through, being injected into an arbitrary input and exiting an arbitrary output port.
In particular, we consider here an integrated-optics implementation of the Clements interferometer where beam-splitters are, in practice, waveguide directional couplers (see for instance the layout of a 6-modes interferometer in Supplementary Figure 5).It is in fact possible, with the femtosecond laser writing technology, to implement static networks of directional couplers with controlled phases and splitting ratios [12].We note that reconfigurable versions of the same scheme, equipped with thermo-optic phase-shifters [74], require beam-splitters to be replaced by Mach-Zehnder interferometers.However, substitution of directional couplers with Mach-Zehnder devices further increases the length of the circuit and does not change how the length scales with the number of modes.
From Supplementary Figure 5 we observe that, excluding the terminal waveguide segments, the overall length of an integrated . Layout for universal interferometer.. Schematic of an interferometer with m = 6 modes, realized in integrated optics according to the Clements scheme.Black thick lines represent waveguide paths.Geometrical dimensions mentioned in the text are indicated.An S-bend curve, with its local coordinate system (ξ,υ) is also highlighted in blue color.
Clements interferometer equals the sum of m−1 times the length L S of a waveguide S-bend, and m times the length L C of the interaction region of the directional coupler (we may consider all couplers to have the same interaction length, as in Ref. [12]).
The technological platform chosen to implement the waveguide circuit determines the minimum curvature radius R min that can be employed, above which additional losses due to curvature are negligible or at least tolerable.Once fixed such minimum radius, the length of the S-bend depends on the transverse elongation h that it has to cover (the wider the elongation, the longer the S-bend).In fact, choosing a sinusoidal shape for the S-bend, the curve may be described by the equation: in a suitable (ξ,υ) local coordinate system.The minimum radius of curvature occurs at ξ = 0 and ξ = L S .By imposing that such minimum radius of curvature is precisely R min , the length of one S-bend in the ξ direction is then given by the formula: Directional couplers are formed by waveguides brought close one to the other at a small distance d for a length L C , and brought sufficiently far apart elsewhere, in order to quench the evanescentfield interaction, which typically decays exponentially with the distance.In particular, interaction must have substantially decayed at the distance p which is the pitch of the input and output waveguide segments.If d p, the elongation of the S-bends inside the interferometric network can be approximated just as h=p.Therefore, we may consider: In a practical case, taking the example of femtosecond-laser written circuits, d may be in the order of a few microns and p at least in the order of a few tens of microns.The length of the interaction region L C of an individual directional coupler determines its reflectivity according to: where c is the coupling coefficient between the waveguides which depends on the distance d, and φ 0 takes into account coupling that occurs in the terminal parts of the S-bent waveguides.As a first approximation, we can consider φ 0 =0, and we can assume that the interferometer is designed as in Ref. [12], hence: Overall, the estimated length of the m-mode Clements interferometer will thus be given by: It is worth noting that this length depends on the waveguide platform employed, which determines the order of magnitude for R min , p and c.On the other hand, different platforms may be characterized by different amount of propagation losses for a given circuit length.
In any case, for any waveguide platform, the circuit length given by Eq. ( 6) scales linearly with the number of modes.
For instance, in the case of femtosecond-laser written circuits it may be reasonable to take R min = 30 mm, p = 0.06 mm and a maximum exploitable c = 1 mm −1 .This makes the two terms of the sum in Eq. ( 6) comparable, and results in a circuit footprint that increases of about 4.5 mm per each added mode to the interferometer.
We could further note that, as mentioned, if we considered a fully reconfigurable circuit designed with the same architecture, each directional coupler should be replaced by a Mach-Zehnder interferometer, thus doubling the optical depth with respect to the case studied above.The precise dimensioning of this other circuit would deserve a separate and detailed discussion.As a rule of thumb, since the optical depth doubles, we could consider an increase of the circuit length of about 9 mm per each added mode, using analogous fabrication parameters.

Planar continuously-coupled interferometer
Integrated-optics platforms enable the realization of interferometers that have no analogous with bulk components, where several waveguides are placed parallel one to the other, and are continuously coupled by evanescent-field interaction.A general algorithm has still to be devised that teaches how to implement an arbitrary unitary transformation of the optical modes with this kind of devices.We will thus proceed with discussing the size constraints of these interferometers on the basis of heuristic arguments.
We can first consider an infinite array of identical optical waveguides placed parallel, on a plane, at the same distance d.Propagation of photons in the array can be described by coupled differential equations involving the destruction operators âx of the optical modes.In the approximation of nearest-neighbour coupling, these equations read: where z is waveguide direction, c is the coupling coefficient which depends on the relative distance between two neighbouring waveguides, and x is the mode index (see also Supplementary Figure 6a).We can look for a plane-wave solution of Eq. ( 7) of the kind: where β z and β x are components of the wavevector respectively parallel and transverse to the waveguide direction z.Note that in the parallel direction we measure the spatial propagation with z, that has the dimensions of a length, while for the transverse direction we use x, which is the waveguide index and is adimensional.By substituting Eq. ( 8) into Eq.( 7) we derive a sort of dispersion relation [75], that governs light diffraction in such a discretized setting.In particular, an excitation characterized by a set of transverse components centered around β x , will propagate across the waveguides of the array with a group velocity: Namely, after a distance z, the center of the wavepacket will have travelled across ∆x = v x z waveguides.Note that there is a maximum achievable group velocity |v x | max =2c.
Let us consider now a photon that enters the array localized on a single waveguide mode.The state of this photon, being spatially pointlike (though in a discretized setting) contains all possible transverse components, including the ones travelling at |v x | max .The wavepacket will spread in the array along propagation in z; however, it will not be able to reach a waveguide that is placed m positions apart in a propagation length smaller than: If we neglect boundary effects and we apply these considerations to a finite-size m-waveguide array, we can take this L m as an estimate a) x − 1  12) in the case of a Bravais lattice with primitive vectors parallel to the axis x and y.In the general case, these axis may not be orthogonal.This scheme can be applied e.g. to a triangular lattice (c): in this case the directions x and y are separated by an angle of 120 • and one can consider cA =cB =cC =c, while cD =0.
for the minimum length that such array must have, so that photons entering in any given input port can interfere together in any output port.
We note that this minimum length L m is shorter than the length of the Clements interferometer as expressed by Eq. ( 6), if the same coupling coefficient c is used in both circuits.However, a continuouslycoupled waveguide array of length L m cannot reproduce the full set of unitaries allowed by the Clements architecture; a longer array with random modulations of the optical parameters may be needed to achieve similar purposes [76,77].In addition, the length L m in Eq. ( 11) scales linearly with m exactly as L Clem in Eq. ( 6).

3D continuously-coupled interferometer
The femtosecond-laser-writing technology has the unique capability to inscribe waveguides at different depths below the glass surface, and thus to realize three-dimensional waveguide arrays with arbitrary cross-section [78][79][80].Light propagation in these arrays is studied analytically by extending the coupled-equations formalism discussed above for the planar devices.In particular, Ref. [81] reports general solutions for the coupling configuration schematized in Supplementary Figure 6b, which is suitable to describe waveguides arranged according to a two-dimensional Bravais lattice with (possibly nonorthogonal) axes x and y.Thus, Eq. ( 7) is generalized as follows: where x and y are waveguide indices on the two directions of the Bravais lattice, c A and c B are coupling coefficients between neighbouring sites on these two separate directions, c C and c D take into account also diagonal coupling.Plane-wave solutions of Eq. ( 12) take the form: âx,y =Ae i(βzz+βxx+βyy) (13) which accounts for distinct β x and β y on the two Bravais directions.The case of an infinite square lattice, with negligible diagonal coupling, corresponds to c A =c B =c and c C =c D =0.In this case, the following dispersion relation is retrieved: The group velocity of a wavepacket propagating only along the x or y direction takes an expression analogous to Eq. ( 10), with the same maximum value Therefore, if we compare a planar array and a three-dimensional array with a square-lattice waveguide arrangement, having the same coupling strength, we observe that light can propagate transversally across the waveguides of the former with the same maximum velocity experienced in the latter, along the x or y direction taken separately.However, in the case of point-like excitation of the square lattice, light travels simultaneously in both directions thus spreading to a quadratically larger number of waveguides.
We can transfer these considerations to a finite-size array, again neglecting boundary effects for simplicity.In a finite-size square lattice with a total of m waveguides, the maximum waveguide indices along the two directions x and y are proportional to √ m.Therefore, in order to spread on the full area of the array, a wavepacket will have to propagate along a transversal direction across a number of waveguides proportional to √ m.It is thus reasonable to assume that the minimum array length, required for allowing a photon injected into an arbitrary input waveguide to reach an arbitrary output waveguide, is: where B is some constant.We note that here L m scales more favourably with respect to the case of the planar array, being proportional to √ m and not to m.An even faster light spreading can be reached in triangular arrays (see Supplementary Figure 6c), as the ones we adopted in the experiments described in the Main Text.In this case c A = c B = c C =c while c D =0, and the dispersion relation reads [79,81]: Note that, due to the definition of the Bravais vectors, here x a y are non-orthogonal directions, oriented with a relative angle of 120 • .The group velocity along the x or y direction is given by: with a maximum value |v x | max =|v y | max =4c, which is twice that retrieved for linear or square lattices.Scaling properties of this minimum length L m , namely the dependence on √ m, are analogous to the ones discussed for the squarelattice case and resulting in Eq. (15).However, in the triangularlattice case L m should roughly halve its value with respect to the former case.
All considerations made up to now regard homogeneous arrays, in which coupling coefficients are uniform across the array and along the waveguide direction.In addition, the evaluated lengths L m are only related to the spreading of the light to all output waveguides, with no constraints on the unitary transformation provided by the circuit.As a matter of fact, homogeneous arrays produce highly symmetric transformations, which may not preserve the complexity of the Boson Sampling problem [76].
In our work, we are indeed introducing static randomness in the continuously-coupled interferometer by taking slightly different coupling coefficients c between different waveguides, namely by varying the relative interwaveguide distances, and by further modulating them along z.In addition, we use an array length that is more than twice the estimated L m with our experimental parameters, in order to promote further mixing of the light across the waveguide array.In fact, in the case of our device, with m=32 and an average coupling coefficient c∼0.2mm −1 , L m is in the order of 15 mm while the waveguide array reaches the length of 36 mm.The considerations about fast light spreading in the arrays, presented in this Section, have guided us in the choice of the triangular array geometry with respect to other possible interferometer layouts.
We may note that previous theoretical works have shown that a random time-modulation of the site energies, in a quantum walk on a one-dimensional chain of sites, is able to provide Haar-random transformations of the input states in the long-time scale [76,77].In a photonic setting, such a quantum walk can be implemented by a photon propagating in planar waveguide arrays, where the propagation constant of the waveguides is modulated randomly along the z direction, and in a different way in each waveguide.A random transformation would then be achieved with a sufficiently long device.In our experiments, fluctuations of the waveguide positions in the cross-section of the array provide a static random modulation of the coupling coefficients, while non-uniform modulations in the propagation constants are provided in a reconfigurable way through the thermo-optic effect, by employing the resistive micro-heaters patterned on the chip surface.The achieved randomness has been investigated experimentally, as described in the Main Text.

Fan-in and fan-out sections
A waveguide interferometer reasonably needs to be coupled, at the input and output ports, with optical fibers.To connect several fibers in parallel to the same facet of the optical chip, commercial fiber blocks or fiber arrays are a convenient choice.These components typically contain the desired number of fibers, arranged on a line at a fixed pitch p F , whose precise alignment is guaranteed by means of V-grooves fabricated lithographically.
In the case of the Clements interferometer, one may choose to design the pitch p ≡ p F , so that fibers can be coupled in a direct fashion.However, in the general case, and especially in continuously-coupled waveguide arrays, the pitch of the waveguides is different, and likely quite smaller, than the pitch of the fiber array.In addition, in case of three-dimensional arrays, waveguides are not even arranged along a single line.In such circumstances, fan-in and fan-out sections need to be added to the circuit, which are composed of S-bends that bring the waveguides at the correct spacing and with the correct arrangement to be interfaced with the fibers.
It is important to study how the length scales also for these parts of the circuit.In fact, we should at least check that advantages in compactness gained in the interferometer section are not vanished by long fan-in and fan-out sections.
As previously discussed, the length of sinusoidal S-bends is governed by Eq. ( 2).If the fan-in or fan-out is judiciously designed, it is reasonable to assume that the maximum lateral elongation will be smaller than half the width of the fiber array, namely h max < p F •(m−1)/2.This constrains the length of the fan-in or fan-out sections to: We note that this bound for L F scales with √ m as m grows larger.Therefore, as the number of modes increases, the length of these sections tends to be less relevant than the length of a Clements interferometer, which scales as m.On the other hand, it may share a similar scaling law as a three-dimensional continuously-coupled interferometer.
An improved compactness may be gained if the m fibers to be coupled are arranged on a two-dimensional grid, as in the experiment described in the Main Text.In this case, the maximum size of the fiber array cross-section, in either dimension, scales as √ m.This maximum size determines also the maximum elongation of the Sbend.Applying Eq. (2), it follows that the length of the fan-in or fanout section here scales as 4  √ m.This further contributes to limiting the device footprint and, hence, its optical insertion losses.

II. RECONSTRUCTION ALGORITHM
A fundamental point for the characterization of reconfigurable photonic chip and for Boson Sampling itself is the reconstruction of the associated unitary matrix.In the case of continuous-coupling devices it is not trivial to find an analytical model linking the elements of the unitary matrix to the internal architecture of the device.For this reason, we model the sub-matrix elements in the most general way, i.e as complex numbers U lm = ρ lm e iφ lm identified by the modulus ρ lm and the phase φ lm .We start by considering two quantities related to twofold experiments.The first is the probability to detect two distinguishable particles in the output i,j, when the particles are injected from the input h,k We observe that the values of a hk ij correspond to a direct measurement of the unitary matrix moduli.The second quantity is the visibility of the Hong-Ou-Mandel (HOM) dip, defined as the difference between a hk ij and the probability of detecting two indistinguishable particles divided by a hk ij : The visibility V hk ij is sensitive to the difference among the phases of matrix elements.The 496 values of the a hk ij and V hk ij for a fixed input pair h,k can be retrieved from the experimental data via different procedures.The a hk ij can be calculated from the intensity distribution of classical light or from the distribution of heralded single-photon injected in h and k.The visibilities can be derived from direct measurements of the two photons probabilities in two different positions of the delay lines.Alternatively, the interpolation of the data resulting from a complete scan of the position x of the delay lines in the HOM dips with a Gaussian function provides an estimation of the quantities under investigation.We find the moduli of the matrix elements by minimizing the χ (2) quantity between the measured a hk ij and the right side of Eq. ( 19).For the phases we exploited an analytical algorithm proposed in Ref. [50] to solve the equation (20).One limitation of this algorithm is that it is not robust with respect to experimental imperfections.For this reason, we use such solution only as a starting point for a further minimization to find the phases φ lm by fixing the moduli to the values retrieved from the first minimization.The quantity to optimize is where hk ij are the experimental error associated to the quantity a hk ij 1+V hk i,j .Furthermore, this last minimization enables us to reconstruct the 3 × 32 and 4 × 32 sub-matrices from a subset of the possible input pairs necessary to evaluate the phases values.In fact, using such pool of data, the analytical algorithm provides two or more solutions that differ only in the sign of the phases.The minimization discriminates which solution provides a better agreement with the experimental data.
In Supplementary Figure 7 we report some examples of 4×32 sub-matrix reconstructed using the method described in this section.

III. UNITARY TRANSFORMATIONS
In this section we describe the preliminary tests performed with the photonic chip regarding the ability to cover a large amount of unitary transformations and to reproduce a given distribution in the outputs over time.In order to check the last feature of the device, we injected a single heralded photon in a fixed input port of the chip, and we collected the output intensities of the 32 modes.Thus, we obtained the squared moduli of a single column of the implemented unitary.We repeated the measurement three times in different days, preserving the same setting for the applied currents.In Supplementary Figure 8 we report the measured output distributions.We quantified the reproducibility through the similarity between the three patterns, whose average is 99.8%.
For what concerns the reconfigurability of the chip, namely the capability to implement different transformations on the input state, we performed the same measurement of the stability test.In this case we varied the applied currents in the circuit through a uniform sampling of the dissipated electrical power.The collected output distributions are included in Supplementary Figure 9.
In Fig. 3b-c we have reported the distributions of the squared moduli and the phases of the various experimental matrices and compared them with likewise matrices extracted according to the Haar distribution.There are slight discrepancies between the experimental and theoretical distributions.In this section we investigate a simplified model to identify the experimental imperfections that generate such deviations from the expected distributions.The histogram of the squared moduli in Fig. 3b is more peaked towards zero than the theoretical one.This is likely due to errors in the estimation of input and output losses.To support this hypothesis, we performed a numerical simulation in which an error of at most 10% on the estimated values of the losses was inserted in the Haar random matrices.The results are shown in Supplementary Figure 10a.The green distribution that displays losses is in good agreement with the experimental one reported in the main text.For what concerns Fig. 3c , we observe a slight peak around zero in the experimental phase distribution.In this case, such discrepancy can be attributed to residual correlations between the phases of the neighboring waveguides and to the reconstruction method of unitary matrices.To study this aspect we have performed another numerical simulation that exploits a simplified model of our device.We consider the case in which the first neighbouring couplings are static along z and the whole surface is uniformly heated.This produces a linear gradient of the temperature at different depths in the sample, which in turn generates correlated changes in the propagation constants and thus in the matrix phases.Note that this condition is quite far from the experiment in which each heater has been controlled independently.The red curve in Supplementary Figure 10b corresponds to the resulting ideal sub-matrices from Haar-random-extracted unitary transformations (red) and the same sub-matrices with randomly modulated insertion losses in the input and output stages (green).The lossy case reproduces the experimental distribution reported in Fig. 3b in the main text.(b) Numerical simulation for the sub-matrix phases distributions in the presence of correlations due to uniform heating of device surface.In red the resulting phases distribution and in green the same phases expressed by using the phases of the first column and row as reference.In this second case the distribution is not flat.
phase distribution.Interestingly, this distribution is completely flat as for the case of Haar random matrices.This highlights the fact that a flat distribution of the unitary matrix phases alone is not a sufficient proof for the sampling from the Haar distribution.The green histogram, on the other hand, corresponds to the same ensemble of matrices, multiplied both at left and at right by diagonal matrices of unit-valued complex elements (equivalent to phase shifters placed at the inputs and outputs of the device), in such a way that the phases of the first columns and rows are set to zero.The latter zero-valued phases were not included in the histogram distribution.Such scenario reproduces what we did in the reconstruction algorithm of the experimental matrix; in fact, the measured HOM visibilities are not sensitive to the input and output individual phases and this allows us to set them to 0 in the corresponding first row and column of the matrix.This procedure generates the slight concentration of the phases in zero, as shown in Supplementary Figure 10b.

IV. SINGLE-PHOTON SOURCE
In this section we illustrate the adopted model used for the four-photon states generated by the spontaneous parametric downconversion in a BBO source.The double-pair emission by a nonlinear crystal can be considered as two independent emission processes by likewise sources.Then, the state produced by each source is , where g ij is the nonlinear gain for the source emitting in the modes (ij) = {(1,2), (3,4)} (see Supplementary Figure 11).For a complete description of the state, the losses, labelled by η i in Supplementary Fig-ure 11, were taken into account.In our experimental setup the corresponding losses are associated to the coupling in single-mode fibers (particularly in delay-lines).Then, the input state resulting from the product ρ in ∼|ψ 1,2 •|ψ 3,4 , includes the contributions with different number of photons, weighted by coefficients that depend from g (ij) and η i .The state, post-selected by measuring four-fold coincidence, has the following form when expressed in the occupation numbers |n 4 ,n 1 ,n 2 ,n 3 , i.e through the number of photons in the corresponding input mode with α=g 12 g 34 η 1 η 2 η 3 η 4 (24) The above coefficients can be retrieved from direct measurements of two-fold coincidences which are related to the contributions |1001 and |0110 in ρ in .In fact, the ratio between these counts is R = g34η3η4 g12η1η2 .Then, the three coefficients have a straightforward expression in R, namely γ =1, α=R and β =R 2 .
In the 3-photon Boson Sampling experiment we post-selected the contribution |1111 by detecting one photon in mode 4 and the other three photons in the chip's outputs.In this case, the input state does not depend from the relative weights among the various contributions in Eq. (23).The estimation of α, β and γ is pivotal in the 4-photon experiments, where the whole state is injected in the optical circuit.We have characterized all coefficients by measuring the R parameter from the two-fold coincidences before the chip.This preliminary characterization of the source was necessary for the validation of the 4-photon samples.
We conclude this section by providing some information about the generation rate of the source.The rate of twofold generation was ∼ 18 kHz in the input modes 1,2 and ∼ 16 kHz in 3,4.For what concerns the Boson Sampling experiment, we had an average rate of ∼300 events per hour for the threefold case and ∼60 for the fourfold events after the chip.

V. BOSON SAMPLING VALIDATIONS
In the experiment we have benchmarked the integrated device by performing several Boson Sampling experiments with 3-and 4-photon states.We repeated the measurements for 10 different configurations of the optical circuits for the 3-photon case, and further three configurations for the 4-photon state in Eq. (23).We adopted likelihood ratio tests to assign the data to a given hypothesis.First, the data were validated against the uniform sampler [54].This algorithm requires the estimation of the quantifier P = i j |U ij | 2 where the index i labels the modes in which photons are detected, the index j the input modes and U the unitary matrix representing the circuit.The counter W initialized to zero is updated after the measurement k according to the following rule where n and m are the number of photons and modes in the optical circuit.The intuition behind this method is that the quantifier P reflects somehow the probability to observe the outcome k.If such quantity is greater than the uniform probability is plausible that the event was sampled from a nontrivial distribution.The second test applied to the Boson Sampling data regards the validation against the distinguishable particles [55].In this case the quantifier is the ratio between the probability q =|PerU (ij) | 2 to detect indistinguishable particles, and the probability d=Per|U (ij) | 2 to detect distinguishable particles in the set of output modes j given the input modes labelled by i. U (ij) stands for the sub-matrix identified by the input labels i and output labels j and Per is the matrix permanent.By defining L = q d , the counter C is updated after each outcome k from the Boson Sampling as Note that the expressions of p and q are related to the collisionfree subspace accessible in the reported experiment.Furthermore, both quantifiers P and L depend from the the element of the matrix U representing the interferometer.In the main text we have shown how a incorrect reconstruction of the unitary matrix affects the outcome of the test.
In Supplementary Figures 12-13 we report the validations tests performed for the two 4-photon experiments not included in the main text.In Supplementary Figures 14-15 we have reported the nine 3-photon Boson Sampling validations whose slope values are reported in Fig. 4f of the main text.We report the validations against the distinguishable particle hypothesis for the state in Eq. ( 23).These further 4-photon experiments were not reported in the main text.We report the validations against the uniform hypothesis for the state in Eq. ( 23).These further 4photon experiments were not reported in the main text.

Figure 1 .
Figure 1.Boson Sampling in a 3D continuous-coupling integrated device.a) Boson Sampling (BS) and the most recent variants, b) Scattershot Boson Sampling (SBS) andc) Gaussian Boson Sampling GBS).The corresponding computational problems require sampling from the output distribution using different input quantum states of light, such as Fock states in BS, two-mode squeezed vacuum states in SBS and single-mode squeezed vacuum states in GBS.The common element among the schemes is the optical random circuit, described by the unitary evolution U. d) The most widely adopted decomposition of the operator U is via a network of beamsplitters, with splitting ratios Sij, and phase-shifts φi (left); an alternative implementation exploits continuous-coupling by evanescent waves among waveguides (right) depending on the coupling coefficients cij and the propagation constants ki, where both may vary along the direction z. e) Overview of the reconfigurable 3D integrated photonic chip, realized through the femtosecond laser writing technique.The device is composed by 32 optical-modes arranged in a triangular lattice, as showed in the inset reporting the transverse section of the sample.In red we have highlighted the input modes employed in the 3-and 4-photon experiment.The transformation U is controlled by the 16 resistors fabricated on top of the glass sample.The second inset shows the top view of the electrical circuits that controls the currents Ii applied to the resistors.

Figure 2 .
Figure 2. Scheme of the experimental apparatus.a) A parametric down-conversion process in a Beta-Barium Borate crystal generates one-and two-pair photon states.Generation in the single-pair regime is employed for the unitary reconstruction procedures, while 3-and 4-photon states are employed for experiments in the Boson Sampling framework.Photons are prepared in their polarization and temporal degrees of freedom before coupling in the input single-mode fiber array.After evolution, photons are finally detected via a set of 32 single-photon avalanche photodiodes connected to a 32-channel time-todigital converter for the reconstruction of the coincidence pattern.b) Schematic of the in-and out-coupling of the single photons with the 3D photonic circuit: one-and two-dimensional fiber arrays connect to the fan-in and fan-out sections of the circuit.Legend: BBO -beta-barium borate crystal, BPF -band-pass filter, HWP -half-wave plate, PBS -polarizing beamsplitter, PC -polarization controller, DL -delay line, PS -power supply, TDC -time-to-digital converter, 1D-/2D-FA -one-/two-dimensional fiber array, FI -fan-in, FO -fan-out, CCWL -continuously-coupled waveguide lattice, RS -resistors.

Figure 3 .
Figure 3. Experimental reconstruction of the [3×32] sub-matrix and comparison with the Haar-random matrices.(a) Experimental reconstruction of the squared moduli ρ 2 ij (red) and phases θij (blue) for three input ports of the reconfigurable photonic chip, highlighted by the labels in the figure.The chip is set on a random configuration of currents.Each input port represents a row of the unitary transformation applied to the input state.(b)-(c) Comparison of the phases and the squared moduli frequency distribution respectively, between 15 experimentally reconstructed sub-matrix (blue) and 15 [3×32] sub-matrix sampled from the Haar-random unitaries (red).(d) Distribution of the similarity between the squared moduli of two columns of different sub-matrix.We repeated the measurement for ∼ 200 different configurations of the currents in the chip.In red it is shown the theoretical similarity distribution obtained by sampling columns distributed according to the Haar measure.The overlap between the two histograms is the 62.4% of the total area.

Figure 4 .
Figure 4. Boson Sampling data.Validation of the 3-photon (blue points in (a)) and 4-photon (blue points in (b)) Boson Sampling experiments against the uniform sampler hypothesis (green point in a-b).(c-d)Validation against the distinguishable photons sampler of 3-and 4-photon events.The test is applied using the same reconstructed sub-matrix in (a) and (b).The distinguishable samples (red points) were collected by adjusting the relative time delay between the input photons.(e) Histograms of simulated 3-photon slopes for counter W (normalized to the slopes of distinguishable particles data) in the case of validation against the uniform sampler using a random extracted U transformation, i.e.U matrices that do not match the actual operation implemented in the circuit.The obtained slopes of the 10 different 3-photon BS experiments are highlighted by the arrows (and by the corresponding stars).The experimental points were validated with the reconstructed sub-matrix retrieved from the two-photon reconstruction.All the experimental slopes values are in the positive range and far from the histograms average, thus showing the correct validation of the performed experiment.In (f), we performed the same analysis for the validation against the distinguishable photons hypothesis.

F
.H. and S.P. contributed equally to this work.F.H., T.G., G.C., N.S., F.S., S.P., A.Cr., and R.O. conceived the experiment.S.P., Z.-N.T., F.C., A.Cr., and R.O. fabricated and characterized the integrated device using classical optics.F.H., T.G., M.I., C.E., A.Ca., G.C., N.S., and F.S. carried out the quantum experiments and performed the data analysis.All the authors discussed the results and contributed to the writing of the paper.SUPPLEMENTARY INFORMATION: RECONFIGURABLE CONTINUOUSLY-COUPLED 3D PHOTONIC CIRCUIT FOR BOSON SAMPLING EXPERIMENTS I. FOOTPRINT OF INTEGRATED-OPTICS INTERFEROMETERS

Figure 6 .
Figure 6.Geometries for continuously-coupled waveguide arrays.(a) Scheme of the cross section of a planar array of identical waveguides.Optical modes (indexed with x) are coupled with a coupling coefficient c.(b) Concept scheme of the couplings considered in Eq. (12) in the case of a Bravais lattice with primitive vectors parallel to the axis x and y.In the general case, these axis may not be orthogonal.This scheme can be applied e.g. to a triangular lattice (c): in this case the directions x and y are separated by an angle of 120 • and one can consider cA =cB =cC =c, while cD =0.

Figure 7 .
Figure 7. Reconstruction of the unitary matrix.Examples of the squared moduli and phases retrieved through the algorithm.Each case corresponds to different configuration of the circuit.The rows represent the four input of the photons and the 32 columns are the output modes of the device.

Figure 8 .
Figure 8. Reproducibility of the implemented transformations.We report the intensity patterns of the output modes for a fixed circuit configuration when a single heralded photon is injected in one input port.The average similarity between two of the three maps is 99.8%.The enumeration on the left map denotes the respective output mode's intensity.The other two maps follow the same enumeration.

Figure 9 .
Figure 9. Reconfigurability of the photonic chip.Here we report the intensity patterns of the 32 output modes (which are illustrated accordingly to the labelling in the top left pattern) for 3 sets of circuit configurations and for 3 different input ports.Thus, we can observe the capability of the chip to implement different unitary transformations while varying the applied currents.

Figure 10 .
Figure 10.Discrepancy from the Haar random matrices distributions.(a) Comparison between the squared-moduli distributions of 15 [3×32]ideal sub-matrices from Haar-random-extracted unitary transformations (red) and the same sub-matrices with randomly modulated insertion losses in the input and output stages (green).The lossy case reproduces the experimental distribution reported in Fig.3bin the main text.(b) Numerical simulation for the sub-matrix phases distributions in the presence of correlations due to uniform heating of device surface.In red the resulting phases distribution and in green the same phases expressed by using the phases of the first column and row as reference.In this second case the distribution is not flat.

Figure 11 .
Figure 11.Scheme for the double-pair emission process.The source is split in two independent double pairs emissions characterized by two nonlinear gain gij.Each source generates photon in pair of optical modes.The latter are identified by propagation losses ηi.Legend: BB0 -Beta-Barium-Borate crystal; BPF -Band Pass Filter; HWP -Half Wave-plate; PBS -Polarizing Beam-splitter.

Figure 12 . 4 -
Figure 12. 4-photon experiments validation.We report the validations against the distinguishable particle hypothesis for the state in Eq. (23).These further 4-photon experiments were not reported in the main text.

Figure 13 . 4 -
Figure 13.4-photon experiments validation.We report the validations against the uniform hypothesis for the state in Eq. (23).These further 4photon experiments were not reported in the main text.