Introduction

A particularly striking departure from classical physics in quantum mechanics is the existence of two fundamental particle classes characterised by their behaviour under pairwise exchange (). For example an arbitrary number of bosons (obeying Bose-Einstein statistics: ‘+’) can occupy a single quantum state, leading to bunching1 and Bose-Einstein condensation; whereas fermions (obeying Fermi-Dirac statistics: ‘−’) are strictly forbidden to occupy the same state via the Pauli-exclusion principle2. Behaviour intermediate between bosons and fermions can also be considered; particles confined to two dimensions known as anyons3 are predicted to exhibit intermediate behaviour with fractional exchange statistics (, 0 < ϕ < π)—a concept applied in superconductivity4 and used to explain the fractional quantum Hall effect5.

Here we show that by tuning one phase parameter, entanglement and a single experimental setup can exactly simulate the measurement statistics of an arbitrary mode transformation A (Fig. 1 (a)) of N identical fermions, bosons or particles with fractional statistics. It is well known that the symmetric and anti-symmetric states |ψ±〉 = |H〉1 |V〉2 ± |V〉1 |H〉2 of two polarised photons input into a 50% reflectivity beam splitter exhibits fermion-like anti-bunching and boson-like bunching, respectively6. Furthermore, intermediate behaviour for a beamsplitter has been observed by controlling a phase to tune between the two Bell states |ψ+〉 and |ψ−〉7,8. It is less widely appreciated that the resulting detection statistics can also be exactly mimicked with two photon path-number entangled NOON states9; however this behaviour does not generalise to processes other than beamsplitters, the simplest example being the identity operation for which bunching always occurs regardless of the phase of the NOON state. The approach we present shares N-partite, N-level entanglement, parameterised by phase ϕ, across N copies of A (Fig. 1 (b)).

Figure 1
figure 1

Simulating quantum interference of quantum particles.

(a) Quantum interference of N identical particles with arbitrary statistics launched into a network described by matrix A leads to quantum correlated detection events at the output. (b) Experimental simulation of quantum interference can be directly achieved by sharing N-partite, N-level entanglement across N copies of A, where the choice of simulated particle statistics is controlled by the phase of initial entanglement.

Full size image

After presenting the scheme in the case of two particles in an arbitrary linear evolution, we experimentally demonstrate the simulation of arbitrary particle statistics with a complex photonic network—a continuous time quantum walk10,11. We achieve the required entangled state with post-selected polarisation entanglement shared across the TE and TM modes of an evanescently coupled waveguide array. This enables us to directly observe the evolution of symmetric and anti-symmetric two particle wave functions in a coupled oscillator Hamiltonian by switching the phase of the entangled state between 0 and π. Mixed-symmetric quantum walk behaviour is also observed by tuning the entangled state's phase between 0 and π, corresponding to behaviour intermediate between bosons and fermions. We note the fermion simulation directly maps to the dynamics of two excitations in a chain consisting of ten sites: Using the Jordan-Wigner transformation, excitations correspond to fermion creation operators and the vacuum to the ground state of each spin12. We finally provide a mathematical demonstration of our scheme in the case of an arbitrary number of particles evolving through any linear Hamiltonian. It is important to note that neither the two-particle case, nor the N-particle generalization, require any post-selection and are, in princlple, scalable.

Results

Exchange statistics refer to a wave-function acquiring phase ϕ after interchanging two indistinguishable quantum particles. Occupation of labeled modes j are modelled by annihilation operators aj and creation operators that obey the commutation relations

for j < k. In three dimensions particles are either bosons (ϕ = 0) or fermions (ϕ = π) and successive swapping is represented by the permutation group while the correlated detection statistics are associated respectively to the permanent and determinant of a given matrix13. In one- and two-dimensions anyons3 can exist where swapping two identical anyons ω times (represented by the braid group) results in the combined wave-function acquiring phase ϕ = ωθ where 0 < θ < π and the positive (or negative) winding number ω keeps track of the number of times two anyons braid in trajectory space with clockwise (or anti-clockwise) orientation4.

The quantum correlation function for two identical non-interacting particles (input into j and k of A and detected at outputs r and q) is given by

for bosons14 (ϕ = 0) and fermions13 (ϕ = π). Here, the process A is represented by a matrix whose element Ar,j is the complex transition amplitude from j to r. Non-classical behaviour is evident in the constructive and destructive interference between typically complex terms inside the |·|2, dictated by ϕ = ωθ. By defining fractional exchange algebra (1–3) for 0 < θ < π and some integer winding number ω, we can observe the effect of fractional exchange statistics on quantum interference between two indistinguishable anyons undergoing braiding operations dependent upon A; Eq. (4) represents the resulting quantum interference between two abelian anyons (0 < θ < π) in superposition of braiding ω times and not braiding at all. To see this we define without loss of generality the input labels by j < k and the output labels by r < q and consider two identical anyons created at modes j and k, transformed according to A and detected at modes r and q. Correlated detection at r and q reveals two indistinguishable possibilities that correspond to the two terms in Eq. (4): A maps the anyon at j to mode r and the anyon at k to mode q with amplitude Ar,j Aq,k; or A maps the anyon at j to mode q and the anyon at k to mode r via ω clockwise braids with fractional phase ϕ and amplitude Aq,jAr,k. From this definition, we use relations (1–3) to derive expression (4) for 0 < θ < π and ϕ = ωθ (see Supplementary Information).

To physically simulate quantum interference given by Eq. (4), we use two copies of A, denoted Aa and Ab, together with bi-partite entanglement of the form where creation operators and are associated to identical networks Aa and Ab respectively; creation operator sub-indices denote the modes on which the creation operators act. Once |ψ(ϕ)〉 is generated, the class of input entangled particles becomes irrelevant; only one particle is present in each of Aa and Ab, therefore no further quantum interference between particles will occur. This implies application with other quantum platforms, using for example entangled fermions to simulate Bose-Einstein quantum interference. Evolution of the two particle state through the tensor product of transformations is therefore modelled by

It follows that correlated detection probability across outputs r of Aa and q of Ab simulates Eq. (4) according to

where the 1/2 pre-factor is due to normalisation of the initial entanglement. Since the phase of the entangled state can be defined as ϕ = ωθ, this method simulates the interference statistics of two anyons of fractional phase θ that have not braided, with two anyons also characterised by θ that have braided ω times (equivalently anyons that are in a superposition of braiding ω1 and ω2 times for ϕ = (ω1 − ω2)θ—the ideal statistics are the same).

We experimentally simulated Eq. (4) for arbitrary ϕ in a coupled oscillator Hamiltonian equivalent to sharing entanglement across two continuous time quantum walks, realised with arrayed waveguides10 (see Supplementary Information). We measured the correlations of the polarisation entangled state shared across two copies of a process A: Aa and Ab, realised as the TE and TM modes of a waveguide array fabricated in silicon oxynitride (SiOxNy) (Fig. 2 (b), see also Methods). The waveguides and butt-coupled polarisation maintaining fibres (PMFs) are birefringent to ensure negligible cross-talk between horizontal and vertical polarisations. Detection statistics of the output photons in the {|H〉, |V〉} basis were analysed using polarising beam splitters (PBS) and pairs of multi-mode fibre coupled avalanche photo-diode single photon counting modules (SPCM) (Fig. 2 (c)). By inputting the state (|H〉j + |V〉j)(|H〉k + e−iϕ|V〉k)/2 and detecting one |H〉 photon and one |V〉 photon at any pair of outputs of Aa and Ab we post-selected the input state , with a measurable phase off-set (from birefringence of PMF and the waveguide circuit). This phase was compensated using a polarisation phase rotation (Fig. 2 (a)). Note that because of this compensation, the approach described is immune to birefringence of the 2 m PMF and input waveguides.

Figure 2
figure 2

Experimental setup.

(a) The parametric down-conversion based photon pair source with quarter-waveplates (QWP), half-waveplates (HWP) and a polarisation beamsplitter cube (PBS) used as a polariser, to control the initial state before input into the chip. Down conversion in the nonlinear bismuth borate (BiBO) crystal yields degenerate photon pairs (808 nm, filtered with 2 nm full width, half maximum interference filters, IF) collected into polarisation maintaining fibres (PMF). See Methods for further details. (b) Two copies of the quantum walk unitary are realised by accessing the TE and TM modes of the polarisation preserving waveguide array. Three input modes of the array are used in the experiment (−1,0,1) accessed via waveguide bends that fan to a pitch of 250 μm and butt-coupled to an array of PMF. Output is accessed via a fan of waveguide from the coupling region to a pitch of 125 μm. See Methods for details. (c) The detection scheme consists of five PBS and single photon counting module (SPCM) pairs. PMF coupled to the chip launch the output photons onto each PBS which separates TE and TM modes guided in the waveguide. The two outputs of each PBS are collected into multi-mode fibre (MMF) coupled to SPCMs for correlated detection via counting electronics.

Full size image

We measured two sets of correlation matrices for two different inputs: immediate neighbouring waveguides in the centre of the array (j, k) = (−1, 0) (Fig. 3); and next but one neighbouring waveguides (j, k) = (−1, 1), with vacuum input in waveguide 0 (Fig. 4). Each entry is the probability to detect a horizontal and vertical polarised photon at outputs r and q, respectively. Five cases are shown for ϕ = 0, π/4, π/2, 3π/4, π (left: (a), (c), (e), (g), (i)). The corresponding ideal correlation matrices given by in Eq. (4) are also plotted (right: (b), (d), (f), (h), (j)). For the (j, k) = (−1, 0) input (Fig. 3) we observe the expected bunching behaviour of bosons for ϕ = 0 which is contrasted with the observation of the Pauli-exclusion in the diagonal elements of the matrices for ϕ = π. We note that a form of exclusion on the diagonal entries would also occur for experiments with hard core bosons, but the remaining correlation matrix would still be dependent upon Bose-Einstein exchange statistics. The spectrum of phase 0 < ϕ < π exhibits a variable Pauli-exclusion and a continuous transition between simulation of Bose-Einstein and Fermi-Dirac quantum interference. In the bosonic case, we observe a less pronounced bunching behaviour for the (j, k) = (−1, 1) input (Fig. 4 (a)) than for the (j, k) = (−1, 0) input (Fig. 3 (a)), as expected14. This behaviour is inverted in the fermionic case, since the particels tend to propagate to the same peak of the walk (while still not occupying the same waveguide) for the (j, k) = (−1, 1) input (Fig. 4 (i)) and tend to propagate in different peaks for the (j, k) = (−1, 0) input (Fig. 3 (i)). This is the result of the different phases in the Hamiltonian that describes the quantum walk evolution.

Figure 3
figure 3

Experimental simulation of particle statistics correlation functions for input (−1, 0).

Experimental correlation matrices (left: a,c,e,g,i) for simulating quantum interference of non-interacting particle pairs with five different particle statistics. Ideal correlation matrices for this experiment are plotted (right: b,d,f,h,j) using Eq. (4). The data ranges from corresponding to Bose-Einstein statistics (ϕ = 0, (a,b)) through fractional statistics (ϕ = π/4, (c,d); ϕ = π/2, (e,f); ϕ = 3π/4, (g,h)) to Fermi-Dirac statistics (ϕ = π, (i,j)). (k) The correlation matrices are displayed as normalised probability distributions with white for probability and black for . Solid coloured squares represent measured correlations, hatched squares correspond to coincidences that cannot be measured due to mis-match between fibre array separation and waveguide pitch (see Methods). Coincidence rates are corrected for relative detection efficiency.

Full size image
Figure 4
figure 4

Experimental simulation of particle statistics correlation functions for input (−1, 1) with data laid out as for Fig. 3.

Full size image

Agreement between the ideal implementation of our scheme, with perfect realisations of Aa and Ab and the data in Figs. 3,4 is quantified in Table I using the similarity S between the ideal correlation matrix Γϕ and the measured distribution of detection statistics Pϕ: , a generalisation of the average (classical) fidelity between two probability distributions that ranges from 0% (anti-correlated distributions) to 100% for complete agreement. We also calculated S when we deliberately made the photons distinguishable by temporally delaying one with respect to the other. This enabled physical simulation of the classical correlations ΓC = |Ar,jAq,k|2 + |Ar,kAq,j|2 that are independent of particle statistics. For these measurements, we found agreement with the ideal model of S = 97.4 ± 0.7% and S = 97.5 ± 0.6% for inputs (−1, 0) and (−1, 1) respectively. These data indicate that Aa and Ab are near identical while keeping crosstalk between horizontal and vertical polarisation to a minimum. While non-unit similarity (S = 100%) in Table I is attributed to non-perfect entanglement generation, the high level of agreement between the ideal and measured cases indicates accurate quantum simulation of arbitrary quantum interference in a uniform coupled oscillator Hamiltonian.

Table 1 Similarity S for each experimental simulation of particle quantum interference with parameter ϕ, for input (−1, 0) given in Fig. 3 and input (−1, 1) given in Fig. 4
Full size table

Using multi-partite entanglement allows simulation of quantum interference of an arbitrary number N of identical non-interacting particles. N indistinguishable particles with exchange statistics ϕ, subjected to mode transformation A at inputs , are correlated due to quantum interference across outputs according to

where σν is an element of the permutation group SN acting on and τ (σν) is the minimum number of neighbouring swaps that maps σν to ; denotes the element j of σν. For ϕ = 0 and ϕ = π, Eq. (7) is the modulus-square of the permanent and determinant respectively of an N × N subset of A. For 0 < ϕ < π Eq. (7)—which is the modulus-square of the imminent of the N × N subset of A—is derived in the Supplementary Information. The multi-particle correlation function Eq. (7) is simulated by subjecting the generalised N-partite, N-level entangled state

to N copies of A such that one particle is launched into each copy A(j). For example, for inputs j, k, l and N = 3, the required entangled input state is equivalent to three entangled qutrits, shared across three copies of A denoted Aa, Ab and Ac:

N-fold coincidence detection across outputs μj of A(j) equals . Note as for the case with two particles, the 1/N! pre-factor can be canceled by considering all N! permutations of the outputs, therefore directly simulating quantum interference of N bosons (ϕ = 0) or N fermions (ϕ = π). For example, in the 2-particle case, the 1/2 pre-factor in Eq. (6) is canceled by summing both the detection outputs of a particle at r on Aa, q on Ab and at r on Ab, q on Aa. All other values of ϕ allow continuous transition between Fermi and Bose quantum interference. Note that the approach of initialising the simulation with a pre-symmetrized (or pre-antisymmetrized) state is in the spirit of first quantised quantum simulation15, for which quantum algorithms have been designed to efficiently generate the required states16,17 (see Supplementary Information). In a photonic realization of the proposed sheme, the required N-particle entangled state can be generated using additional photons and feedforward to implement quantum gates following the KLM scheme18.

Discussion

Correlations resulting from quantum interference of multiple indistinguishable non-interacting bosons in large linear optical circuits cannot be efficiently simulated with classical computers19,20,21. Large scale repeated quantum interference is achievable with photons in waveguides11,22,23,24, which is likely to provide the platform for the realisation of such a boson scattering device, to realise an engineered system that can outperform the limits of classical computation. A major issue will be verifying boson scattering correlations that exceed the capabilities of conventional computers. In contrast, the analogous situation for non-interacting fermions is efficient with classical computers (currently it is not as straightforward to perform quantum interference of fermions2,25 as it is with photons in integrated optics). The ability to tune continuously between tractable and intractable phenomena with a single setup could be essential for verification of bosonic scattering.

Other possible applications for simulating particle statistics with photonics include preparing appropriately symmetrized states in quantum chemistry simulations15 and combinatorics, where exchange statistics have been shown to alter mathematical graph structure26—three-dimensional waveguide architectures27 may be important here. Our simulation of arbitrary fractional exchange statistics also provides a novel method to study fractional quantum interference phenomena28; whether our approach can be used in the context of topological fault tolerant quantum computation with photons29,30 is an open question.

The method we report can be directly applied to any quantum process realisable with linear optics31 and is not reliant on using polarisation—the entangled state could be encoded in time-bin or interferometrically stable path with reconfigurable waveguide circuits32. In order to introduce interactions between simulated particles, it may be possible to use effective non-linearities using auxiliary photons18,33. In the long term, the method we have demonstrated here could be applied to other quantum platforms (consisting of fermions or bosons) where nonlinear behaviour may be more readily exploited.

Methods

The photon pairs are generated in a 2 mm thick, nonlinear bismuth borate BiB3O6 (BiBO) crystal phase matched for non-colinear type-I spontaneous parametric downconversion, pumped by a 60 mW 404 nm continuous wave diode laser. The photons are spectrally filtered using 2 nm full width at half maximum interference filters centred on 808 nm and collected with aspheric lenses focusing on polarisation maintaining fibres (PMFs). The typical two-photon count rate directly from the source with two 60% efficient single photon counting modules (SPCMs) is 50 kHz. For a coherent input state, arrival time of the two photons at the input of the waveguide array is made identical via characterisation by performing a Hong-Ou-Mandel experiment1. The polarisation state is then prepared for both photons using a half wave-plate with optic axis set to π/8 radians. The polarisation in the left arm is matched and purified to that of the right using a rotated polarisation beam splitter in transmission. The two photon state launched into the waveguide array is then . The collecting PMF and the birefringence of the waveguide array preserve horizontal and vertical polarisation components for both photons and prohibit crosstalk between horizontal and vertical light; PMF-coupled polarisation beam splitters at the output of the waveguide array allow two-photon correlated detections across the two orthogonal polarisations to project the initial input state to be the entangled state ; the extra phase is compensated using a Z(θ) rotation on the right arm of the photon source comprised of two quarter waveplates (optic axis at π/4 radians) with a half waveplate between (optic axis at radians, see Fig. 1 of the main text).

The waveguide array was fabricated from silicon oxynitride (SiON) thin films were grown by plasma enhanced chemical vapour deposition (PECVD) onto thermally oxidised <100> oriented 100-mm silicon wafers applying an Oxford 80 Plasmalab Systems. A deposition process based on 2% SiH4/N2 and N2O was performed at 60 W power, 187.5 kHz frequency, 650 mTorr pressure and 300°C substrate temperature. The SiON films were annealed for 3 hours at 1150°C in a nitrogen atmosphere. By optical inspection at 632 nm wavelength a film thickness of d = 603 (±3) nm and a refractive index (for TE) of nTE = 1.5232 (±1 × 10−3) were measured. The channel waveguides were obtained from standard lithography followed by reactive ion etching in CHF3/O2 chemistry at 350 W plasma power. The etch depth (d = 613 ± 5 nm) was confirmed by profilometry. A 5-μm thick upper cladding layer was deposited in two steps applying a Low Pressure CVD TEOS-process and standard PECVD SiO2 deposition. The upper cladding was annealed at 1150°C for 3 hours. For a core width of 1.8-μm the total birefringence of the waveguide, including material and form birefringence34, is of the order Δn = nTM − nTE = 2 * 10−4. This value provides polarization preserving propagation for the TE and TM photons throughout the device. Moreover, the birefringence is small enough to ensure a nearly identical mode overlap and coupling ratio, between different waveguides in the coupled array for TM and TE modes. This implies that the processes Aa = ATE and Ab = A™ are matched, as required and verified in the reported experiment. The coupled waveguide arrays were accessed with lithographically patterned bends (minimum bend radius of 600 μm) that gradually mapped the 250 μm pitched input (for 250 μm pitched PMF array access) to the 2.8 μm pitch in the evanescent coupling region. The coupling region is 350 μm long, before bending gradually to an output pitch of 125 μm which was accessed by 250 μm pitched arrays of PMF; the output measurements were alternated between access to all odd labeled waveguides and all even labeled waveguides. The optical chips were diced from the wafer.