Quantum walks of two correlated photons in a 2D synthetic lattice

Quantum walks represent paradigmatic quantum evolutions, enabling powerful applications in the context of topological physics and quantum computation. They have been implemented in diverse photonic architectures, but the realization of a two-particle dynamics on a multi-dimensional lattice has hitherto been limited to continuous-time evolutions. To fully exploit the computational capabilities of quantum interference it is crucial to develop platforms handling multiple photons that propagate across multi-dimensional lattices. Here, we report a discrete-time quantum walk of two correlated photons in a two-dimensional lattice, synthetically engineered by manipulating a set of optical modes carrying quantized amounts of transverse momentum. Mode-couplings are introduced via the polarization-controlled diffractive action of thin geometric-phase optical elements. The entire platform is compact, efficient, scalable, and represents a versatile tool to simulate quantum evolutions on complex lattices. We expect that it will have a strong impact on diverse fields such as quantum state engineering, topological quantum photonics, and Boson Sampling.


I. INTRODUCTION
A quantum walk(1) (QW) is the quantum-mechanical analogue of the classical random walk, describing the evolution of a quantum particle that moves on a discrete lattice, hopping between adjacent sites. These quantum dynamics can be either continuous or discrete in time, depending on whether the couplings between neighbouring lattice positions are continuously active or can be described as sudden kicks, occurring at discrete time-steps. In the latter case, at each step the walker moves in a direction that reflects the state of a spin-like internal degree of freedom, playing the role of "quantum coin". The growing interest in QWs is due to their potential use in diverse quantum applications, as for instance quantum search algorithms (2), quantum gates for universal quantum computation (3)(4)(5)(6), quantum state engineering (7)(8)(9)(10), and quantum simulations of topological and physical phenomena (11)(12)(13)(14)(15)(16)(17)(18)(19)(20)(21)(22)(23).
Quantum walk dynamics exhibit a richer variety of phenomena when the evolution involves more than one particle. These walks indeed are characterized by multi-particle interferences (24)(25)(26)(27)(28), having no classical analogue and incorporating an inherent source of complexity, as highlighted by their central role in computational models such as Boson Sampling (29,30). The dimensionality of the lattice where the QW takes place is also a crucial ingredient. As an example, quantum search algorithms based on QWs overcome their classical counterparts solely when the spatial dimension of the lattice is equal or greater than Step 1 Step 2 Step 3 (c) 2D lattice encoding. The walker position on the 2D lattice is encoded in transverse momentum of the beam. We consider a linearly polarized input beam, without transverse momentum component. The corresponding position on the lattice is (0,0). A g-plate, oriented along y (x) with δ = π, divides a horizontally-polarized input beam in two parts, which acquire two opposite transverse momentum components along y (x), k y(x) = ±∆k ⊥ or θ y(x) = ±θ0. Then, after coin tossing, a second g-plate is oriented along the orthogonal direction. Four different beams with different angular deviations (±θ0, ±θ0) are obtained at the output, mapping the lattice positions (m, n) = (±1, ±1), respectively. The red lines on the images indicate the propagation direction of the Gaussian beams. The beam deviations have been overemphasized for the sake of visualization. The deflected beams remain spatially overlapping while they travel along the setup, except in the final imaging stage. two (31). Furthermore, two-dimensional (2D) QWs display a richer landscape of topological features, when compared to the 1D case (11,32).
Photonic platforms developed to implement QWs differ in terms of the methods to encode both walker and coin systems into optical degrees of freedom. Starting from the first experiments in linear optical interferometers composed of beamsplitters and phase shifters (33), integrated photonic technology has enabled significantly larger instances both in their continuoustime (21,26,(34)(35)(36)(37)(38) and discrete-time version (8,25,39,40). Other schemes rely on light polarization and orbital angular momentum degrees of freedom (9,41), multimode fibers (36) or fiber network loops (14,15,(42)(43)(44), where the walker position is simulated by the temporal separation between the laser pulses. Other remarkable experiments have been conducted by controlling confined waves in arrays of micro-resonators (45,46). Multiparticle regimes have been already implemented in continuous-time QWs (47)(48)(49), even in 2D lattices (38). However, the demonstration of multi-photon discrete-time QWs in more than one spatial dimension has remained elusive thus far. Here we devise and experimentally validate a compact, flexible and scalable photonic platform that achieves this goal. Specifically, we realize a three-step quantum walk dynamics with coherent light, one-photon and two-photon inputs, spanning 2, 6 and 12 modes for the first, second and third step, respectively.

A. Model and encoding
The essential elements of a discrete time QW are captured by the single-step evolution operator U , as after an evolution of t timesteps the system is described by a quantum state |ψ t = U t |ψ 0 , where |ψ 0 is the input state. The operator U typically includes a spin rotation C, acting only on the coin Hilbert space, and a spin-dependent shift. When the walker moves on a 2D square lattice [see Fig. 1 For two-photon inputs, both photons are injected in the QW implementation and they propagate along two parallel paths. A lens system enlarges their waist radius and introduces a relative angle between the optical modes. The relative inclinations of the optical modes represent different lattice positions. In this way the photons start the walk at positions (−1, 0) and (1, 0) of the lattice.In the single-particle case, one of the two photons is directly measured to act as a trigger. QW implementations. The quantum walk is performed by using waveplates and g-plates arranged in a cascade configuration. Each g-plate is controlled independently by tuning its phase retardation via a voltage controller. Measurement stage. A lens system at the output stage converts the different momentum values into a spatial grid. Then, for single photon acquisition, an array of micro-lenses is used to efficiently inject the output modes in a 2D square-lattice fiber array. Finally, each output fiber is plugged into an avalanche photo-diode detector connected to a coincidence electronic system. For the coherent light data acquisition, we inserted a beamsplitter between the last lens and the micro-lense array, and we positioned the CCD on the reflected path to perform image acquisition. details can be found in Supplementary Note 1). Our platform builds on a recent approach to the simulation of single-particle 2D QWs using coherent laser light (22). Here, the walker positions are provided by optical modes |m, n with the following spatial profile: where A(x, y, z) is a Gaussian envelope with large beam radius w 0 in the transverse xy plane, ∆k ⊥ represents a quantum of transverse momentum and the z axis is regarded as the main propagation direction. ∆k ⊥ fulfils the condition ∆k ⊥ 2π/λ, λ being the optical wavelength. Modes in Eq.1 are essentially Gaussian beams propagating along a direction that is slightly tilted with respect to the main propagation axis [see Fig. 1(a)]. Photons associated with mode |m, n carry an average transverse momentum (k x , k y ) = (m∆k ⊥ , n∆k ⊥ ). Basis states of the coin space (|↓ , |↑ ) are encoded in right-handed (|R ) and left-handed (|L ) circularly polarized states, respectively. The conditional shift is realized via a liquid-crystal polarization grating, that is a g-plate (22). These plates are made of a thin layer of liquid crystal, whose molecular orientation is arranged in a periodic pattern along one of the directions that are transverse to the propagation direction. Considering for instance a g-plate with a modulation along the x direction, with a period Λ = 2π/∆k ⊥ , its action on spatial modes defined in Eq. 1, combined with circular polarized states, has the following expression: A similar expression holds in case the modulation is along the y axis. The parameter δ is the birefringent optical retardation of the plate, that can be adjusted by applying an external voltage to the outer surfaces of the liquid crystal cell (50). Thus, in our encoding the g-plates implement the generalized shift operators T x (δ) and T y (δ), with the value of δ determining the fraction of the wavefunction that is shifted to neighbouring sites [see Fig. 1(a)]. We use optical waveplates with tunable retardation ω (based on uniformly patterned liquid-crystals) to implement adjustable coin rotations C(ω) [see Fig. 2(b)]. Therefore, the single-step operator is an ordered sequence of g-plates and waveplates. A large variety of QWs can be implemented with this platform, by tuning the parameters δ and ω, as already shown in Ref. (22). In this work, we will focus on a fully-balanced 2D-QW protocol described by the following one-step operator U = T y (π)C(π/4)T x (π)C(π/4) [see Fig. 2(c)]. x y P(x,y) x y P(x,y) x y P(x,y) x y P(x,y)

B. Experimental setup
The complete setup is shown in Fig. 3 (more details can be found in Supplementary Notes 2-4 and Supplementary Figures 1-3). A photon-pair source is employed to generate and then inject single-and two-photon inputs into the quantum walk platform. In particular, in the two-photon case, an appropriate optical system is implemented to inject the two particles in different sites of the lattice, corresponding to different optical modes (see Methods). The same apparatus can be used to inject classical laser light. The quantum walk itself is implemented via a cascade of the single-step building blocks described above. Note that, by turning on and off individual g-plates in sequence, that is by setting δ to π or 0, respectively, it is possible to measure the spatial distribution of the walker after each step of the protocol. Finally, the output of the quantum walk is sent to the detection stage.
In the focal plane of a lens, our modes can be spatially resolved as they form a grid of small spots. A suitable 2D fiber array and a micro-lens system are placed at this position, so that each spot matches the core of the corresponding fiber in the array (see Methods). With a classical light input, the output can be measured via a charge-coupled device (CCD) camera. We have employed our platform to implement a balanced quantum walk U with single-and two-photon inputs, up to three time-steps.

C. Single-photon Quantum Walk
We performed a single-particle experiment by injecting a single photon in position (1,0) with polarization |D = (|H + |V )/ √ 2, and then measuring the output distributions after each step (obtained by sequentially switching on the g-plates). The single-photon experimental distributions are shown in Fig. 4 (related analysis for coherent light input is reported in Supplementary Note 5). The slight differences between theoretical and experimental output distributions are due to experimental imperfections such as inaccuracies of g-plates tuning values and of their relative horizontal alignment. The agreement between experimental data and the expected distributions is quantified by the similarity, defined as: where P (t) (r) andP (t) (r) are the theoretical and experimental distributions of the quantum walk at the t-th step, respectively, while r is the particle position on the lattice. The similarity value of the last step S 1p = 0.9773±0.0002 shows a high agreement with the expected distribution. Similar results were observed when injecting classical light to the quantum walk platform (see Table I).

D. Two-photon Quantum Walk
Next, we realized a two-photon, 2D quantum walk by injecting two photons with polarizations |A = (|H − |V )/ √ 2 and |D in (−1, 0) and (1, 0) lattice positions, respectively. It is worth noticing that at the input of the quantum walk the two-photon state is separable. Temporal synchronization between the particles is ensured in advance by performing a direct Hong-Ou-Mandel (HOM) measurement at the output of the quantum walk, which provides a measured visibility v = 0.95 ± 0.02. Further details on the HOM measurement are reported in Supplementary Note 6 and Supplementary Figure 4. For the multi-photon case, the obtained theoretical and experimental distributions depicted in Fig. 5 show a high quantitative agreement. This is confirmed by their similarities defined as: where P (t) (r 1 , r 2 ) andP (t) (r 1 , r 2 ) are the theoretical and experimental distributions of the quantum walk at the t-th step, respectively, while r 1 and r 2 are the position on the lattice of the two particles. The theoretical distribution was computed by considering an initial state described by the density matrix ρ 0 = c 0 ρ ind + (1 − c 0 )ρ dis , where ρ ind indicates the density  (6) 0.204(2) 96 2 steps ∼ 0.996 0.9929(2) 0.9743 (7) 0.0077(5) 14 3 steps ∼ 0.988 0.9773(2) 0.914(2) 0.0084 (7) 11 TABLE I. Similarities of distributions related to classical, one particle and two particle regime for all the steps. Experimental results for the 2D quantum walk up to the 3rd step when using classical light (S cl ) one photon (S1p) and two photons (S2p). In the last columns we report the maximum value of V with the associated errors σV , and the corresponding values of V/σV for each measured configuration. matrix of two completely indistinguishable photons and ρ dis is the density matrix describing two distinguishable particles. c 0 , that corresponds to the measured visibility of the HOM test, is equal to 0.95. For the experimental distribution at the third step we obtain a similarity value S

E. Non-classical correlation witness
The presence of non-classical correlations in the output distributions is witnessed by applying the non-classicality test of Refs. (48,51), given by: where Γ (cl) m1,m2 is the classical probability that light exits from the m 1 and m 2 output ports of an interferometer. How this formula adapts to our specific case is shown in detail in the Supplementary Note 7. The obtained maximum violation of the inequality are of 96, 14 and 11 standard deviations, respectively for 1, 2 and 3 QW steps, thus unambiguously proving the quantum behavior of the reported two-photon 2D quantum walk. The indistinguishability between injected photons gives rise to quantum interferences, yielding in turn probability distributions that cannot be reproduced by classical states of light. In Fig.  6, we report the complete plots of V over the standard deviation for the second and third step, respectively. All the results of one-particle and two-particle QWs are summarized in Table I. In the last column, we also reported the maximum value of the violation V and its error, for each step. These results highlight that the proposed platform has the potential to be employed for significantly larger instances, with high degree of control on the implemented protocol.

III. DISCUSSION
We have presented and realized a platform for the implementation of two-dimensional, multiphoton discrete-time quantum walks, demonstrating experimentally single-and two-photon operations on a 2D squared lattice. The presented platform is compact, flexible and enables the implementation of a large variety of different topological quantum walks (22). Hence, it can represent a powerful tool for the investigation of rich dynamics that are experimentally unexplored. Recent works reported 2D continuous-time quantum walks of correlated photons, relying on arrays of coupled waveguides (21,38,52) , and important results have been also achieved by means of superconductive quantum processors (6). We stress that our system is based on a very different approach, exploiting a synthetic 2D lattice made of internal modes of a single optical beam, as opposed to real-space neighbouring lattice sites, and implementing discrete-time evolutions that can be actively controlled and easily reconfigured. We believe that these different approaches may have complementary advantages. In our platform, several quantum walk protocols can be dynamically realized. This can be achieved by tuning the retardation of each plate in the range [0, π], by changing their orientation and their position in the plane transverse to the photons main propagation direction. By controlling these parameters, diverse single particle QWs mimicking periodically driven Chern insulators have been reported (22). Here we implement a different split-step quantum walk, that is proved to realize the Grover search algorithm in the high step-number limit (53). The number of steps that can be currently realized is essentially limited by the optical losses, which are mainly due to photon reflections at each plate. However, these can be significantly reduced (from ∼ 15% to at least ∼ 5%) by adding a standard antireflection coating on the plate outer surfaces. Several applications can be foreseen, including quantum state engineering (7)(8)(9) or quantum algorithms based on the quantum-walk paradigm (2,3). Furthermore, given the possibility to exploit multi-photon inputs and to control the performed transformation, this approach can also represent a promising platform for the implementation of Boson Sampling and Gaussian Boson Sampling experiments in large optical lattices (30).

A. Photons preparation
The photon pairs used in our experiment are generated by a parametric down-conversion source, composed of a nonlinear beta barium borate crystal (BBO) pumped by a pulsed laser with λ = 392.5 nm. The generated photons, λ = 785 nm, are then injected in two identical single mode fibers for spatial mode selection. On each fiber, an independent polarization controller allows choosing the polarization of each input photon. Then, delay lines are used to temporally synchronize the optical paths through a Hong-Ou-Mandel interference measurement.
Once the photons are temporally indistinguishable, they are injected in the QW platform. In order to precisely control the distance between the injected photons and match the coupling conditions of the fiber array, a half-mirror on one of the two paths is used. By translating this mirror, the distance between the paths can be modified from 3 mm to 8 mm, while its tilting can change the relative orientation between the two photons. Finally, an appropriate lens system superposes the two paths and introduces a small prescribed angle between them (see Supplementary Note 4). This relative angle corresponds to into a difference in the photon transverse momentum.

B. Measurement stage
At the output of the quantum walk structure, a three-lens system is used to decrease the relative distance between adjacent beams and reduce the corresponding beam waists. Then, a micro-lens array, composed of a hundred micro-lenses with a short effective focal distance (approximately 5 mm), separated by 250 µm, is used to inject the photons into a square-lattice multimode fiber array, reaching an individual waist of 15 µm without changing the corresponding distances. The adoption of the micro-lens array enables an improvement in the coupling efficiency from 0.1 to 0.75. Overall, the fiber array sets a 8 × 8 spatial grid, where the output fibers are separated by a 250 µm pitch. Each of these multi-mode fibers is connected to a single-photon avalanche photo-diode. The output signal of each detector is directed to a coincidence apparatus, able to record single-photon counts and two-fold coincidences.

V. DATA AVAILABILITY
The data that support the findings of this study are available from the corresponding authors upon reasonable request.

VI. CODE AVAILABILITY
The codes for data processing are available from the corresponding authors upon reasonable request.

Supplementary Note 1. Theoretical background
The realization of a quantum walk on a 2D square lattice is performed by the repetitive action of an unitary single-step operator U on a quantum system. Tracing a parallel with the classical random walk, the quantum system plays the role of the walker while its internal spin-like degree of freedom represents the coin [1]. Therefore, the walker can be described by a quantum state |ψ in the bipartite Hilbert space H = H p ⊗ H c . The two dimensional internal degree of freedom space, H c , is spanned by the states {|↑ , |↓ }, and the position space H p by {|m, n }, where (m, n) label the square lattice coordinates. Thus, the discrete-time quantum walk evolution at the t-th step can be expressed as: The single-step operator U encompasses the spin rotation C(ω), which only acts on the coin space, and the spin-dependent translation operators along x and y directions, T d , with d = x, y. The general unbalanced spin rotation is given by : where parameter ω sets the balance of the coin rotation. In a balanced coin rotation, there is an equal probability to move in each direction, condition achieved if ω = π/4. The translation operator can be expressed as: where δ establishes the probability of the walker to remain in the same position at each step, while the operator S d is the shift operator along the d = x, y direction. The S x operator does not act on y coordinates and it is given by: |m − 1, n, ↑ m, n, ↓| + |m + 1, n, ↓ m, n, ↑| .
Similarly, the operator S y does not act on x coordinates and is defined as S y = m,n |m, n − 1, ↑ m, n, ↓|+|m, n + 1, ↓ m, n, ↑|.
The ordering of these operators in the single-step operator U and the chosen values for δ and ω allows to define different QWs protocols.
Here, we mainly focused on the 2D-QW protocol U = T y (δ)C(ω)T x (δ)C(ω), setting δ = π and ω = π/4. By making this choice, we selected a protocol in which the walker must change its position at each step, with a balanced coin rotation.
In the particular case of one-particle regime, we considered a family of localized initial states |ψ 0 described by: From the evolved state, |ψ t , we obtained the probability distribution as: where the vector r indicates the positions (m, n) in the lattice. In the two-particle regime, partial photon indistinguishability has to be taken into account to properly model the input state and calculate theoretical predictions. For two-photon inputs, this effect can be modeled with a density matrix defined as: where ρ ind indicates the density matrix of two completely indistinguishable photons and ρ dis is the density matrix describing two distinguishable particles. The c 0 parameter sets the degree of indistinguishability between the two photons, which is defined in the range [0, 1]. At each step, we obtain a probability distribution on the lattice position of the two walkers as: Tr(ρ t |r 1 , σ 1 , r 2 , σ 2 r 1 , σ 1 , r 2 , σ 2 |).
At each step, the elaboration of the theoretical prediction for the distributions shown in the main text was performed by using Supplementary Equations (6) and (8).

Supplementary Note 2. g-plate Technology
The key elements of our quantum walk are the g-plates [2], which are electrically controlled liquid-crystal waveplates for which the optic axis orientation is varying linearly along one of the in-plane directions. By using the same principle of the q-plates [3], an optimally tuned g-plate splits a light beam into two orthogonally polarized beams, as depicted in Fig.2 of the main text. The relative angle between the two split beams is determined by the diffraction-grating spacing of the g-plate, and it is approximately equal to α = 2λ/Λ, where λ is the wavelength of the input beam and Λ the characteristic parameter of the device, defined during the fabrication process. Thus, such device adds a fixed transverse momentum component (∆k ⊥ = 2π/Λ) to an input beam, and the sign of the "kick" depends on its polarization state. In particular, if the beam is right-handed circularly polarized (|R ), it acquires a negative component of the transverse momentum. On the other hand, a left-handed circularly polarized beam (|L ) acquires a positive fixed component of the transverse momentum. The action of an x-oriented g-plate in the polarization space is given by the following matrix, written in the circular basis (|R , |L ): where δ is the phase delay introduced by the g-plate, that defines each momentum "kick" probability, and it is tunable by an external programmable voltage source. Finally, the accumulated momentum can be visualized in position space by using a lens system that implements an optical Fourier transformation. An overall amount of 12 electrically-tuned liquid crystal devices were exploited to perform a two-photon QW up to the third step, encompassing waveplates and g-plates arranged in a cascade configuration. In Supplementary Figure 1 (a) a schematic representation of the exploited apparatus is shown, while the image of a g-plate is reported in Supplementary Figure 1 (b).

Supplementary Note 3. Micro-lens array
A micro-lens array, model M LA − S250 − f 20 of RPC Photonics, was used to efficiently inject the photons into the 2D fiber array employed at the detection stage. The micro-lens array is made of polymer-on-glass material, has a size of 250 × 250 mm 2 , index of refraction 1.56 (633 nm), and it is composed of a hundred of micro-lenses with a short effective focal length (around 5 mm), separated by 250 µm (see Supplementary Figure 2). In this way, we were able to reach an individual waist of 15 µm without changing the corresponding distances in the experiment. Consequently, the adoption of such device enabled a significant increase in the coupling efficiency into the final array of optical fibers used for detection from 0.1 to 0.75.

Supplementary Note 4. Experimental implementation
The photons used in the experiment were produced by a parametric source, consisting of a beta-barium-borate (BBO) crystal, that generates a photon pair (λ = 785 nm) by non-collinear spontaneous parametric down conversion (SPDC) when it is pumped by a pulsed laser with λ = 392.5 nm. Once the photons were generated by the SPDC source, they were injected into two identical single-mode fibers. Mechanical polarization controllers were employed to arbitrarily and independently control the polarization of each input photon. In the case of the single-particle 2D-QW, one photon was directly sent thorough a single mode fiber to an Avalanche Photodiode (APD), to be used as trigger. The other photon was sent to the platform at the fiber launch "1", shown in Supplementary Figure  3. It performed the 2D-QW by following the same path as for the 2D-QW with two-photons, presented below. On the other hand, for the case of the two-photon 2D-QW, both particles were injected in the platform. In this case, temporal synchronization of the photons was performed via independent delay lines through an Hong-Ou-Mandel experiment (see Supplementary Note 6). The single-mode fibers were then connected to the fiber launchers "1" and "2", depicted in Supplementary Figure 3. A mirror, mounted on a translation stage, and a half-mirror were used to finely adjust the parallelism between the two beams and to change their relative distance d 0 , that can vary from 3 mm to 8 mm. In the sequence, a telescope, composed of two lenses with focal distances f 1 = 300 mm and f 2 = 30 mm respectively, was used to reduce the beam waist by a factor 10. The relation between the waist and distance can be expressed as following: Subsequently, we used a lens with focal distance f 3 = 2000 mm to enlarge each waist without changing the separation length.
At the output of f 3 , each beam had a waist W 2 , respecting the following relations: where z 1 is the Rayleigh length of the beam after the telescope. When z 1 f 3 , this can be accurately approximated as:

2D Fiber Array Detection stage
Supplementary Figure 3. Scheme of the 2 dimensional quantum walk apparatus. In the single-particle regime, one photon was sent to the fiber coupler "1". For 2D-QW with two walkers, the photons entered the system through ports "1" and "2". A system composed of several mirrors and lenses was properly designed to ensure an efficient coupling into the fiber array, considering the physical constraints imposed by its pitch. After the photons propagate through the g-plates, a micro-lens array was placed immediately before the multimode fiber-array to increase the coupling efficiency.
Therefore, the two beams were highly superposed and collimated. By using the ray transfer-matrix formalism, we were able to compute the relative angle θ 2 between the two beams at the output of f 3 , assuming an effective parallelism at the output of f 2 : Hence, by changing the value of d 0 , we linearly controlled the input angle between the two beams, θ 2 . The angular lattice unity is ∆ θ = λ/Λ, where Λ is the g-plate characteristic parameter defined above. For d 0 = 3.14 cm, the two photons are prepared at (0, 0) and (1, 0) modes of the quantum walk (θ 2 = λ/Λ = 1.57 × 10 −4 rad). Thus, increasing by a factor 2 the value of d 0 , θ 2 is also magnified by the same factor. In the latter case, the two photons were prepared at (−1, 0) and (1, 0) input lattice positions. The two-dimensional quantum walk was performed by a cascade of g-plates oriented along the x and y directions, interspersed with waveplates. For n g-plates x and n g-plates y , the maximum angular difference between two beams is θ max = 2 n ∆ θ.
At the output of the quantum walk structure, we used a three lens system to decrease the relative distance between consecutive beams and reduce their waist. The f 4 = 200 mm lens was conjugated with the f 5 = 50 mm one to reduce the waist by a factor 4. Then, we used a f 6 = 400 mm lens in order to set the subsequent relative distances to d 3 and each waist to W 3 . Given two adjacent beams at the output of the g-plate-system, the angular difference ∆ θ induces a spatial shift at the focal point of f 6 : Numerically, we obtained: d 3 = 250 µm and W 3 = 80 µm. Hence, the maximum angular shift became a spatial shift equal to n × 250 µm. Finally, we used a micro-lens array (see Supplementary Note 3) to improve the coupling efficiency into the fiber array. We employed a fiber array defining a 8 × 8 grid, in which the fibers were separated by 250 µm. Each one of these multi-mode fibers has a core of 62.5 µm and its output was coupled to a single-photon avalanche detector.
The measured rate at the output of the SPDC source was approximately 60 kHz single-photon count at each arm, and 6 kHz two-fold events by using a coincidence window of 8 ns. The overall recorded two-fold events after the first, second and third steps of the quantum walk were 46926, 65611 and 90608. The corresponding coincidence-rates after 1, 2 and 3 steps were 56 Hz, 11 Hz and 7 Hz, respectively. The estimated losses for of the optical apparatus are due to the following optical elements: • tunable plate (15% for each plate); • fiber-array coupling (40%); • delay lines (40% for each photon); • bunching probabilities detection (around 50%) (See the Supplementary Note 7); As such, the overall transmission efficiency of the QW evolution (that is considering all liquid crystal devices) is 0.85 N , where N is the number of plates. This yields an efficiency equal to 14% after three steps, where N = 12. However, by putting a standard anti-reflection coating on plate surfaces, the transmission of a single cell can be increased at least to 95% and the overall efficiency to 54%.

Supplementary Note 5. Experimental results for coherent light inputs
For the coherent light data acquisition, we inserted a beamsplitter between the last lens and the micro-lense array, and we positioned the CCD on the reflected path for the image acquisition. We show the raw pictures of one-photon distributions with coherent light in Figure (4b) of the main text. The dark counts of the images are removed by averaging the intensity values of the pixels in a 10 × 10 pixel square that is far more then 15 lattice steps away from the visible spots. We find the center of each spot of the image by looking for the local maxima of each picture. Then, the waist of each spot is estimated by performing a two-dimensional Gaussian fit on the image data. We calculated the intensity of each light spot in the pictures by summing the intensity value associated to each pixel that belongs to the spot. Then, we divided the latter by the total value of the intensity. In this way, we reconstructed the probability distributions for each step.

Supplementary Note 6. Photon synchronization through Hong-Ou-Mandel interference
Temporal synchronization between the two-photon paths was performed via a Hong-Ou-Mandel [4] test. Three liquid-crystal plates were exploited as a beamsplitter, following a procedure similar to that reported in Ref. [5] with q-plates. In particular, we turned on the plates in the following order: a g-plate along x, a waveplate tuned as ω = π/4 and, finally, a g-plate along x. Two identical photons were sent in two adjacent modes (m + 1, 0) and (m − 1, 0) of the first g-plates with polarization |R and |L , respectively. We only considered the motion along the x axis and we indicated the state of the photon in the second quantization formalism as a † ix,P |0 , where i x is the x coordinate on the lattice and P is the polarization state. The action of the three plates, written in the creation-annihilation operation formalism, reads: and: Hence, the two-photon probability of having a bunching event in the initial positions x 1 = m + 1 and x 2 = m − 1 becomes: In Supplementary Figure 4 (a) we show the plots of the bunching probability as a function of the tuning parameters δ 1 and δ 2 .
For δ 1 = δ 2 = π the probability of bunching is equal to 1. For this value of the g-plates tuning, the plates mimic a beamsplitter with transitivity t = 1/ √ 2, in the subspace of g-plate modes (m + 1, R), (m − 1, L): that is the matrix representation of a 50/50 beamsplitter. We experimentally performed the Hong-Ou-Mandel experiment in this configuration. The two photon input state was prepared in lattice sites (−1, 0) and (1, 0) with |L and |R polarization, respectively. In particular, we considered one of the two possible outputs of this unitary evolution to measure the bunching configuration. Indeed, for any arbitrary unitary operation, the probability of bunching of two indistinguishable photons exiting from the same output is twice the one obtained with distinguishable ones. More specifically, by considering a generic unitary U and two photons entering in input modes n an m, the probability that the two particles exit from the same output k in the distinguishable and indistinguishable cases are respectively: independently of the unitary matrix. For this reason, the visibility of the HOM peak is independent of the g-plate unitary transformation and only depends on photon indistinguishability. To discriminate the events when two photons exit from the same output we used a 50/50 fiber beamsplitter (FBS). Then, we connected the FBS output ports to two avalanche photo-diodes (APDs). We then acquired the coincidence-rate at different positions of the translation stage. The measured interference pattern is shown in Supplementary Figure 4 (b). The visibility v of the peak was estimated by using the following formula: where the cc in and cc out are the coincidence-rates at the maximum and large time separation, respectively. The obtained value of the visibility is 0.95 ± 0.02, showing that the degree of indistinguishability of the photons is not significantly affected by the cascade of g-plates corresponding to the quantum walk platform.

Supplementary Note 7. Experimental results for one-and two-photon inputs
Three steps of the U 2D-QW protocol were performed with one and two photons. The one-particle probability distributions are measured by counting the two-fold events between the photon that performs the quantum walk and the trigger photon. Conversely, the two-photon distributions are obtained from the two-fold events between each pair of output modes. We considered only the modes whose coincidence number is much higher than accindental ones (2 modes for the first step, 6 for the second step and 12 for the third step), as analogously performed in previous experiments [6][7][8]. In order to measure also the bunching probability we use a fiber beamsplitter for each measured mode. Hence we use two avalanche photodiode detectors for each mode. The fiber beamsplitter allows one to detect bunching configurations, i.e. two photon per mode, with a probability p = 2t(1 − t), where t is the fiber beamsplitter transmissivity. We calculate the error on each distribution value by using a bootstrapping approach. Then each probability distribution P is obtained by dividing the two-fold events of each pair of modes by the total number of measured events. We remove the accidental coincidences (due to the dark counts) by considering the dark signal as a continuous stochastic signal. The accidental coincidences for each pair of modes A i,j are given by [9]: where S i and S j are the values of single-event rates of the i and j modes, T is the total time of the data acquisition and w is the coincidence windows, i.e. the time interval in which we consider two events as simultaneous. Finally, we correct each value with the coupling efficiencies. These corrections are obtained by injecting coherent laser light into each output mode and by calculating the coupling efficiency as the ratio between the output and the input light intensity.
Here, we report the details on the similarity calculation for one-and two-particle regimes. Moreover, we report also the distributions for distinguishable particles. Finally, we show the plots of the non-classicality test, presented in Supplementary Equation (24), in the indistinguishable regime for all the steps. The errors on distributions values and all the derived quantities were calculated via bootstrapping approach, by considering that the distributions values are sampled on a Poissonian probability distribution.
A. One-particle regime In the one-particle regime, the experiment was carried out by directly detecting one of the two-photons, acting as a trigger, while the other is injected in the quantum-walk platform. The probability distributions were measured at the end of each step, and the steps were performed by inserting the photon with polarization |D in (1, 0). The distributions are shown in the main text. We compared the theoretical and experimental distributions using the similarity defined as: where t is the step number, r is the particle position on the lattice and P (t) (r) andP (t) (r) are the theoretical and the experimental distributions, respectively. B. Two-particle regime In the two-particle regime, the photon paths were synchronized in the g-plate by performing a Hong-Ou-Mandel test (see Supplementary Figure 4). As explained in Supplementary Note 6, such test was performed by exploiting the g-plates as a beamsplitter. The value of the measured visibility was employed to set the value of c 0 = 0.95 in the distribution theoretical calculations. Similarly to the one-particle regime, the probability distributions were measured at the end of each step. The steps were performed by injecting the two photons with polarization |A and |D in (−1, 0) and (1, 0), respectively. The twodimensional lattice was linearized to perform a representation of the four-dimensional probabilities, thus presenting the latter quantities as two-parameter functions. The linearization of each lattice site was obtained as the following map r → m + 7n with m, n ∈ [−3, 3] . To report the experimental coincidence-rate matrix, the same linearization was performed. In the main text, we show distributions after each of the three steps for indistinguishable particles, with the corresponding value of similarity. The similarity was calculated by generalizing that corresponding to one particle in the following way: where t is the step number, r 1 and r 2 are the positions on the lattice of the particle 1 and 2 and P (t) (r 1 , r 2 ) andP (t) (r 1 , r 2 ) are the theoretical and the experimental distributions, respectively. In Supplementary Figure 5, the theoretical and experimental distributions for distinguishable particle (c 0 = 0) for 1, 2 and 3 steps are shown. The calculated similarities are S   Figure 5. Theoretical distribution and experimental reconstruction of the 2D-QW with distinguishable particles. The theoretical distribution was computed by setting c0 = 0, corresponding to the visibility of distinguishable particles. The same linearization procedure employed in the main text for indistinguishable particles was performed to obtain a representation of the experimental data. Shaded regions on top of each bar correspond to the experimental error at one standard deviation. The error bars were obtained thorough a bootstrapping approach. The bunching probabilities are highlighted by the yellow line.

C. Violation Inequality
In order to estimate the non-classicality of our experimental distribution, we employed the inequality introduced in Supplementary Ref. [10,11]. More specifically, classical light has to satisfy the following condition: where Γ (cl) m1,m2 is the classical probability that the light exits from the m 1 and m 2 output ports of an interferometer. The action of the interferometer is described by a unitary matrix U. The evolution of 2D-QW with g-plates is equivalent to the action of an interferometer with unitary matrix given by U = U t at each step t. The output modes of 2D-QW interferometer are equivalently given by m = r, σ where r is the position on the lattice and σ = ↑, ↓. By measuring the final distribution of 2D-QW we trace away the polarization state. For this reason we do not have a direct access to Γ r1,σ1,r2,σ2 . Hence, we recast Supplementary Equation (24) in terms of the measured values of Γ r1,r2 = σ1,σ2 Γ r1,σ1,r2,σ2 . In order to do that, we sum the violation inequalities for the fixed positions r 1 and r 2 on the lattice: σ1,σ2 V(r 1 , σ 1 , r 2 , σ 2 ) = 2 3 σ1,σ2 Γ (cl) r1,σ1,r1,σ1 Γ (cl) r2,σ2,r2,σ2 − Γ (cl) r1,r2 < 0.
By exploiting the inequality: