Quantum simulations with multiphoton Fock states

Quantum simulations are becoming an essential tool for studying complex phenomena, e.g. quantum topology, quantum information transfer, and relativistic wave equations, beyond the limitations of analytical computations and experimental observations. To date, the primary resources used in proof-of-principle experiments are collections of qubits, coherent states or multiple single-particle Fock states. Here we show the first quantum simulation performed using genuine higher-order Fock states, with two or more indistinguishable particles occupying the same bosonic mode. This was implemented by interfering pairs of Fock states with up to five photons on an interferometer, and measuring the output states with photon-number-resolving detectors. Already this resource-efficient demonstration reveals new topological matter, simulates non-linear systems and elucidates a perfect quantum transfer mechanism which can be used to transport Majorana fermions.


I. INTRODUCTION
Quantum simulations boost the development of topological materials [1], quantum transport [2] and quantum algorithms [3] for the benefit of low-power electronics [4], spintronics [5] and quantum computing [6]. They employ intricate quantum interference of light or matter particles. This is a challenging task: the difficulty arises from the fundamental constraint that all interfering quanta must be indistinguishable [7]. Violating this demand precludes the observation of such coherent phenomena in larger scales, in terms of particle number and duration.
So far, protocols have mainly relied on the use of three distinct quantum states: numerous qubits implemented by superconducting circuits [8] and electronic states of trapped ions [9]; coherent states of photons [10] and atoms (Bose-Einstein condensates) [11]; and multiple single-particle Fock (number) states distributed among many modes in photonic waveguides [12] and optical lattices [13]. Thus, simulations have never seriously profited from interference of multi-particle Fock states, even though the importance of this regime has been recognised [14], and the first attempt to mimic it with manybody systems was made [15].
Here we experimentally and theoretically demonstrate that multiphoton Fock state interference can be useful for quantum simulations that address applications of high impact. Remarkably, this approach grants access to a non-linearity induced by photon number detection [16] and also avoids error accumulation that weakens methods using quantum walks, built on numerous steps [17]. Our idea, shown in Fig. 1a-c, is based on overlapping two multiphoton Fock states, |l a and |S −l b (l photons in mode a and S − l in mode b), on a beam splitter with tunable reflectivity r which programs the simulation duration. We then collect photon statistics at its outputs.
The primary example of a system we can simulate is a chain of S + 1 two-level spins that initially contains just one spin excited, and that is subjected to an XY interaction. The excitation probabilities at its sites after the interaction duration are determined by the output photon statistics. These mappings are based on a solid mathematical grounding known as the Schwinger representation which links quantum harmonic oscillators with representations of spin Lie algebra su (2).
Our platform also allows us to simulate certain classes of fermionic systems, e.g. a non-linear Su-Schrieffer-Heeger (SSH) model [18], obtained from the XY spin chain by a Jordan-Wigner transformation. Furthermore, we can map to Bogoliubov-de Gennes Hamiltonians, simulating many body systems beyond the single excitation subspace e.g. a p-wave superconducting chain (Kitaev model) [19], and the transverse-field Ising model. Due to recent advances in photon-numberresolved detection, we were able to employ transition-edge sensors (TESs) [20] to count photons exiting the beam splitter.
Amazingly, TES measurements correspond to single-site-resolved detection in the chain. The use of TESs is crucial, as Fock state quantum interference is evidenced by photon bunching. For example, two identical photons impinging on a balanced beam splitter leave in a superposition of two-photon Fock states, with both always being detected in the same output port. This is known as the Hong-Ou-Mandel (HOM) effect [21] whose generalised form can be observed for higher-order Fock states if they are prepared in similar polarisation, spectral and spatio-temporal modes [22], as shown in Fig. 1d. c) The likelihood of detecting k and S − k photons in its output ports p(k), simulates the excitation probability of the kth spin in the chain as a result of the interaction. d) The statistics p(k) originates from fundamental indistinguishability of several scenarios that occur to interfering Fock states which classically are exclusive but quantum-mechanically are coexisting, and amount to the same partitioning of incoming photons into two exit ports. Events for which quantum probability amplitudes add up non-destructively are registered more often than others.

II. RESULTS
The Fock state quantum simulations build on a beamsplitter interaction U where a † and b † denote photonic creation operators that act on the interferometer input modes. The reflectivity r, defined as the probability of reflection of a single photon, encodes the interaction time θ(r) = 2 arcsin √ r. For entries |l a and |S − l b , the computational output from the beam-splitter and detectors is [23] k (x, S) is the Kravchuk function [24]. We selected three distinct examples of simulations, shown in Fig. 2, for experimental demonstration. The first a, uses input data initialised to | S 2 a | S 2 b and the setting of r = 0.5. For the second and third b & c, we set |0 a |S b and repeated the computation several times whilst gradually increasing r. While for the second program one can use any value of S, the third one runs exclusively for an odd number of photons.
Edge states in non-linear systems. Interpretation of the outcomes of our quantum programs be-comes straightforward if we consider matrix representations of H BS and of the Hamiltonian describing a general chiral XY 1 2 -spin chain H XY =

S+1 n=1
Jn 2 (σ x n σ x n+1 + σ y n σ y n+1 ), where σ x n and σ y n are the Pauli operators acting on the nth spin. In the single excitation subspace spanned by the states |n = σ + n |↓ 1 , . . . , ↓ S+1 , where σ + n = (1/2) (σ x n + iσ y n ) is the raising operator, the latter has matrix elements [H XY ] Spin mn = m|H XY |n = J n−1 δ n,m+1 + J m−1 δ m,n+1 , where δ i,j denotes the Kronecker delta. Spin nm , when we set the spin couplings to J n = 1 2 n(S + 1 − n). As these amplitudes are non-periodic, this chain lacks translational invariance. This precludes the usual Fourier space methods used for characterising topological insulators. Remarkably, photon statistics measured behind the beam splitter is capable of simulating this noncrystalline system. The existence of topologically nontrivial states is indicated here by the fact that the Hamiltonian belongs to the chiral orthogonal (BDI) class of Altland-Zirnbauer symmetry classes, characterised by a Z topological invariant. Our first program performs a real-space study of this system and computes probabilities that describe its zero-energy eigenmode, FIG. 2: Encoding the outcomes of quantum simulations in photon statistics. We repeatedly overlap two Fock states on a beam splitter of reflectivity r to collect the photon statistics in its exit ports, p(k). They directly provide the results of computation carried out by quantum interference. a) The first program uses | S 2 a | S 2 b and r = 0.5, revealing weakly localised edge states with a non-decaying envelope (black line), which closely resemble topological states in a non-linear SSH model. b) The second program is run for |0 a |S b and several values of r, demonstrating that perfect quantum wave packet transfer in a linear register results from mirror reflection of the input state w.r.t. the register centre. c) For an odd S, this program additionally simulates the perfect transfer of Majorana fermions an and bn in a p-wave superconductor over S+1 2 atomic sites. The bars located at even k (light green, blue, and red) correspond to the mode an with n = k/2 + 1, while those located at odd k (dark green, blue, and red) to the mode bn with n = (k + 1)/2. like the typical edge states which are exponentially peaked at the ends of a quantum wire, these two edge states are weakly localised and plateau to a constant value in the bulk, given by  1−(2k/S−1) 2 , as outlined in Fig. 2a. The intensity-dependent amplitudes J n render H XY a generalisation of the seminal Su-Schrieffer-Heeger (SSH) model [18] to the non-linear regime [25]. See Supplementary Material for details.
Perfect state transfer. The XY spin chain with these couplings has been extensively studied in the literature due to its remarkable property of facilitating the perfect transfer of an arbitrary quantum state [26]. Our quantum simulation provides new insight into this system from which the perfect transfer becomes self-evident. The equivalence of H BS and H XY matrix representations implies the correspondence between interactions generated by these Hamiltonians, U (r) BS and U XY (t) = e −itHXY , respectively. Mathematically, the beam-splitter interac-tion in the Fock state basis amounts to an α-fractional Quantum Kravchuk-Fourier transform (α-QKT) of the input state with fractionality α = 4 π arcsin √ r [23]. As 2-QKT is the spatial inversion operator [24], so is U XY (t) at t = π. Therefore, the transfer is an effect of mirror reflection of a quantum state w.r.t. the chain centre. Proving this fact was tricky within the framework of spin chains, whereas it is an evident conclusion from our photonic simulations. We note that α = 2 implies interference on a perfectly reflecting beam splitter (r = 1) which swaps input states at its outputs. To demonstrate this behaviour, in our second program, we simulated the state transfer of a strongly localised edge state, typical of e.g. the SSH model. The initial Fock state |0 a |S b is gradually transformed to |S a |0 b for increasing r, as shown in Fig. 2b.
Generalised Majorana modes. Multiphoton Fock state interference also facilitates the simulation of manybody systems that are not restricted to a single excitation subspace. For example, a p-wave superconducting chain (Kitaev model) [19] is described by the mean field where c † n and c n are creation and annihilation operators for electrons on the nth atomic site, while µ n , t n and ∆ n are site dependent chemical potentials, hopping amplitudes and superconducting pairing potentials respectively. This Hamiltonian may be expressed in the form T is a Nambu spinor and H BdG is the Bogoliubov-de Gennes Hamiltonian matrix, in the basis of Majorana operators a n = c n + c † n and b n = i(c † n − c n ). The beam splitter Hamiltonian in the Fock state basis H BS is identical to H BdG for the parameters µ n = J 2n−1 , t n = ∆ n = J2n 2 , where 2N = S + 1. This correspondence allows one to simulate the Heisenberg evolution of the Majorana operators over the interaction time θ(r), as well as the evolution of the real fermion operators c n and c † n , by using linear superpositions of Fock states as input. In particular, the evolution of the operators a n and −ib n is encoded by the evolution of the photonic modes |2(n − 1) a |S − 2(n − 1) b and |2n − 1 a |S − (2n − 1) b respectively. To evidence this, a further simulation with input |0 a |S b was performed, where S is an odd number, modelling the perfect transfer of Majorana modes between the two ends of a p-wave chain of N = S+1 2 atomic sites that is depicted in Fig. 2c. This is half the number of sites as in the XY spin chain, reflecting the fact that each physical fermion comprises a pair of Majoranas. The simulated dynamics also apply to one-dimensional arrays of photonic cavities [27] where the effective superconducting pairing and Majorana modes arise from Kerr-type non-linearities within a Bose-Hubbard model.
Non-uniform transverse-field Ising chain. One can also simulate a transverse-field Ising model, H K = 1 2 N n=1 (µ n σ z n +2t n σ x n σ x n+1 ), since this is related to the pwave superconducting chain by a Jordan-Wigner trans-formation. Due to the non-uniform field µ n and spin couplings t n , the system inherits the perfect mirror reflection from the beam splitter dynamics and allows for perfect state transfer after an interaction time θ = π. We thus highlight a new quantum spin network that allows perfect transfer, similar to the previously discussed XY model, but which has not been considered by previous authors. For an example, to simulate the transfer of an excited spin between ends of a chain, one should interfere the state 1  We characterised the set-up to confirm the high degree of indistinguishability of these Fock states, the key issue for multiphoton HOM effect. We measured the standard HOM interference dip between both sources for a small mean photon number of the order of 10 −4 , and achieved the visibility V HOM = 85.9%. Next, we took a measurement of the second order correlation function for each SPDC source separately and observed g (2) ≥ 1.86 ≈ 1 + V HOM , which corroborates the previous result. An effective Schmidt mode number of K = 1 g (2) −1 = 1.16 proves our SPDC sources nearly single-mode.
The measured simulations are presented in Fig. 4. The data shown in Fig. 4a and Fig. 4b consists of approximately 1.6 × 10 3 registered events, for each value of r, in which the total number of photons was S = 4. The data in Fig. 4c comprises 2.3 × 10 2 measurements, for each value of r, in which S = 5. We compared them with a numerical model based on Eq. (2) supplemented with the analysis of experimental imperfections, and found that they are in a good agreement. Errors were estimated as a square root inverse of the number of measurements. See Methods for details.
In Fig. 4a we show the photon statistics recorded by TES 2−3 for the coupler splitting ratio r = 0.5, conditioned on the heralded photon numbers l = 2 and S − l = 2 in modes c and d. They directly model a zeroenergy eigenmode of a non-linear SSH model described The VC allows us to set the ratio between 0 and 1 with an error of ±1.5 × 10 −2 . We used transition-edge sensors (TESs) with efficiency exceeding 90% for photon-number-resolved measurements in all modes [20]. The optimal temporal overlap at the VC is achieved by adjusting an optical path delay τ . The data is analysed with a data acquisition unit (DAQ).

by [H BS ]
Fock nm , with emerging two weakly localised edge states. Fig. 4b depicts the statistics gathered for l = 0 and S − l = 4 for several splitting ratios: r = 0.04 (green squares), 0.3 (orange triangles), 0.5 (blue circles) and 0.96 (red diamonds). It visualises perfect state transfer of the first spin excitation in the chain of 5 particles by means of continuous-time mirror reflection w.r.t. the chain centre. Fig. 4c shows an experimental simulation of the perfect transfer of a Majorana fermion in a pwave superconducting chain of 3 sites that is based on the statistics gathered for l = 0 and S − l = 5 for all the listed values of r.

IV. CONCLUSIONS AND OUTLOOK
Multi-particle Fock state interference is a new and compelling method in the field of quantum simulations, promising for studying non-crystalline topological materials, beyond the recently challenged bulk-edge correspondence theorem [28]. It allowed us to simulate systems as diverse as an XY spin chain and a non-linear SSH model, as well as the perfect transfer of Majorana fermions over a quantum wire, in a system that is not tied to a single-excitation subspace. Importantly, photon-number-resolved detection introduces an effective non-linearity [16] which can be harnessed in simulated models. The presented examples apply to a variety of systems such as superconducting nanowires [30], disordered graphene quasi-1D nanoribbons [31] and disordered cold atoms [32]. These may find applications in nextgeneration electronics [33] and spintronics [34] operating with almost no energy dissipation and speeds exceeding 100 GHz. Our simulations amount to computations of the Kravchuk transform that classically is expensive but in the quantum domain can be attained with a single gate [23].
Multiphoton Fock states have been utilised in quantum simulations in a very limited capacity until now. In photonics, the main focus was on successful manipulation of large numbers of single or pairs of photons in bulk optics [35], as well as in integrated platforms [36]. For example, only one-and two-photon output states of a quantum walk in coupled waveguides were measured, which are a small fraction of the total output [12]. The advantage of these photonic systems, however, lies in easily engineered waveguide layouts which can be used to e.g. model different couplings in the chains. Nevertheless, keeping a high degree of indistinguishability of photons coming from different sources remains a challenge [7].
On the contrary, our simulations are the first to be done exclusively in Fock space, with Fock states of high photon number encoding all the information from input to output. Although currently the experimental gener- ation of five-photon Fock states is already beyond the state of the art, it is soon expected to reach the level of tens of photons [37]. Moreover, our method avoids some of the error accumulation and scaling problems of the waveguide-based set-ups. It can also be extended to higher dimensions by including additional degrees of freedom such as photon frequency and polarisation. The scope of simulations could be further broaden by using input states superpositions S l=0 x l |l, S − l and altering the spin-chain couplings. Although preparation of such general superpositions poses a challenge in photonics, input states in the form of generalised Holland-Burnett states were experimentally obtained by interfering Fock states on a beam splitter [38]. Some other examples could be reached by heralding and conditional state preparation using more intricate interferometers. Merging our approach with coupled-waveguide set-ups is yet an unexplored and intriguing territory.
It would also be very interesting to implement our technique with quantum simulation platforms that are uni-versal. For example, Fock states are also available in motional states of trapped ions up to 10 excitations [39] and in the form of plaquette Fock states of atoms in optical lattices up to 4 excitations [15]. The range of accessible parameters controlling these systems could provide access to other complementary simulation models. Moreover, deterministic creation of an arbitrary superposition of Fock states has been demonstrated for trapped ions and superconducting resonators [40]. This would further expand the assortment of input states that could be used for simulation and may give birth to new fascinating results.

A. Characterisation of the set-up
Each integrated SPDC source produced a two-mode weakly squeezed vacuum state |Ψ = ∞ n=0 λ n |n, n s,i , where s and i denote two output modes, named the signal and idler, λ n = tanh n g cosh g , |λ n | 2 is a probability of creation of a pair of n photons and g is the parametric gain. The average photon number in each mode of |Ψ is n = sinh 2 g. The observed average photon number of n ≈ 0.2 amounts to g = 0.44, which was sufficient to ensure the emission of multiphoton pairs. In this regime one can approximate cosh g ≈ 1 and thus, λ n ≈ sinh n g = n n . The transition-edge sensors (TESs) were operated at 70 mK, which allowed photon-number resolved measurements in all modes [20].
The transmission losses in the set-up were estimated by means of Klyshko efficiency measurements. To this end, we set the reflectivity of variable coupler at r = 0.5, and pumped each of the two SPDC sources separately at successively lower power values. The registered four-mode photon statistics were then binned into 'photon(s)/nophoton' datasets to mimic the use of standard binary detectors, e.g. avalanche photo-diodes, and we concluded the total efficiencies of the heralding modes c and d to be η c = 50.3% and η d = 48.5%, respectively. The variablecoupler modes a and b exhibited a total efficiency of η a = 21.6% and η b = 20.6%, respectively. These values result from the fact that each mode carried a 3 dB loss from the coupler itself and another 1 dB due to coupler insertion and fibre-to-fibre coupling losses. We estimated the transmission losses to be approximately of 50% ≈ 3 dB. Here 1 dB stands for the initial fibre in-coupling loss due to spatial mode mismatch, while 0.25 dB stems from detectors inefficiencies, and the remaining loss is from three FC/PC fibre-to-fibre couplers per mode as well as bending losses in the transmission fibres between the set-up and the detectors.

B. Error estimation
In the experiment, each measurement results in a 4tuple consisting of the number of photons registered by TES 1−4 , corresponding to photon-number states in modes a-d (Fig. 3). The tuple counts are stored in a database. The probability of detecting k and S − k photons in modes a and b is computed as p(k) = N k /N , where N k is the database value retrieved for the key (k, S − k, l, S − l) and N is the total count of events characterised by the given total number of photons S. The measurement errors for each mode were estimated to ∆p = 1/ √ N .

C. Numerical model of experimental outcomes
To assess the experimental results we developed a theoretical model which extended Eq. (2) by taking into account the influence of losses, multi-modeness of beams as well as inefficient photodetection.
Decoherence resulting from the first two effects was modelled by replacing the mode b † with a superposition of the same mode b † and an orthogonal one b † ⊥ , i.e. b † → cos y b † + sin y b † ⊥ , where the parameter y ∈ (0, π 2 ) introduced weights and 'tuned' the distinguishability. This transformation led to the interference of |l a with a two-mode Fock state superposition S−l n=0 S−l n −1/2 cos n y (sin y) S−l−n |n b |S − l − n b⊥ instead of the single-mode Fock state |S − l b , as before. Thus, effectively, some of the multiphoton states interfered with the vacuum state and this implemented the usual model describing particle loss. In our computations, we took y = arcsin (K − 1)/K, where K denoted the effective Schmidt mode number measured during the set-up characterisation.
Realistic model of photodetection requires taking into account a probability of detecting n d photons when a Fock state |n in reaches a TES. It is given by p TES (n in , n d ) = nin n d (1 − η) nin−n d η n d where n d ≤ n in and η is the detector efficiency. In our computations we first used a starting value of η = 0.7 and then numerically optimised efficiencies for individual TESs to compensate for the uneven photon number distribution p(k) seen in Fig. 3a. The programs were written in Python using mpmath library. JJR, SWN, TG, and AL delivered and maintained the transition edge sensor detection system. AB, TS and TM developed the software and performed numerical computations. AB prepared the plots. All the co-authors wrote up the manuscript.
Competing interests: The authors declare that they have no competing financial interests.
Data and materials availability: All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Material. Additional data available from authors upon request.

Appendix A: The Schwinger representation: mapping between photonic and spin platforms
It is a well-known fact in the representation theory of groups and algebras that the Lie algebra su(2) can be represented in terms of annihilation and creation operators of a harmonic oscillator. This is known as the Schwinger representation [41]. It allows one to associate two independent quantum-harmonic oscillator modes with spin operators as follows where S 0 is the Casimir operator S 0 (S 0 + 1) = S 2 x + S 2 y + S 2 z and the standard su(2) commutation relations hold Therefore, we can immediately identify that In this picture, a two-mode Fock state |l a |S − l b corresponds to spin-S 2 particle that is prepared in an eigenstate of S z with eigenvalue l − S 2 , known as a Dicke state Furthermore, one can one-to-one map the Dicke states | S 2 ; m to the basis states that span the single excitation subspace of a spin-1 2 chain. To this end, we employ the following relabelling m = − S 2 + n − 1, where 1 ≤ n ≤ S + 1 denotes nth spin-1 2 in the chain with S + 1 sites. Then, | S 2 ; l − S 2 corresponds to the chain where (l+1)th spin-1 2 is excited where σ + m = 1 2 (σ x m + iσ y m ) is the raising operator acting on mth spin, with σ x m and σ y m denoting the Pauli operators. For a concise notation, we denote such spin-chain states as follows Appendix B: Introduction to Fock-state photonic quantum simulations The key observation that provides the basis for quantum simulations based on Fock-state interference is the formal mathematical mapping between the Hamiltonian matrix representations of an XY-type of interacting spin chain and that of a beam splitter. As we pointed out in the main text, a general chiral XY spin chain is represented by the Hamiltonian In the single excitation subspace that is spanned by the states shown in Eq. (A7), this Hamiltonian has the matrix representation periodic couplings of J n = 2J(1 + δ(−1) n ), where −1 ≤ δ ≤ 1 is the dimerisation parameter. Therefore the chain consists of alternating couplings, one weaker and one stronger. The dimerisation can be chosen in such a way that the atoms at the ends of the chain experience either the weaker coupling (δ > 0) or the stronger coupling (δ < 0), shown in Fig. 5. The first option results in the formation of topologically non-trivial states and the second one trivial states. The topological phase transition takes place at δ = 0.
In the topologically non-trivial phase, the system is a topological insulator with two zero-energy edge states. For example, the left boundary state is of the form where ξ = 2/ ln 1+δ 1−δ is the localisation length [43]. Usually, the boundary states are studied for δ close to 1, when they are exponentially peaked at the ends of the chain, shown red in Fig. 6a. Fig. 6b also shows the zero-energy states but in a less studied regime, close to the topological phase transition, for δ = 0.005 (blue). Interestingly, the conventional SSH model can present weakly localised edge states that are nonetheless topological. An interesting generalisation of the SSH model to the non-linear domain was achieved by setting intensity-dependent site couplings [25]. This model was theoretically implemented with an array of cavities with the tunneling constants 1 Interference of Fock states on a beam splitter can simulate an interacting spin chain with non-periodic nextneighbour spin couplings shown in Eq. (B4) and Hamiltonian given by Eq. (B1). This system is depicted in Fig. 9. Interestingly, the couplings (B4) also are intensity dependent, as n is the intensity of Fock state |n and thus, we also work with a non-linear SSH-type of model. However, unlike in Section C 2, the dependency is not linear.

The BDI Altland-Zirnbauer symmetry class
In the periodic table of topological insulators defined by the Altland-Zirnbauer symmetry classes [44], a system is categorised according to the properties of its time-reversal operator T = U T ⊗ K, charge-conjugation operator C = U C ⊗ K, and chiral-symmetry operator Γ = T ⊗ C, where U T (C) is a unitary operator and K is complex conjugation. If there exists a T (C or Γ) that commutes (anti-commutes) with the system Hamiltonian, than the system is said to possess the respective symmetry and is classified according to the square of that operator.

The beamsplitter Hamiltonian matrix representation [H BS ]
Fock is real-valued and thus, its time-reversal symmetry operator is simply T = 1 ⊗ K. Due to the absence of couplings beyond nearest-neighbour the chiral symmetry operator is given by Γ ij = (−δ ij ) 2 , and finally C = T ⊗Γ. Our photonic system possesses all three symmetries with (T 2 , C 2 , Γ 2 ) = (1, 1, 1). Thus it belongs to the BDI (chiral orthogonal) class of topological insulators which in one-dimension is characterised by a Z topological invariant. Therefore, the simulated spin chain does so too. The zero-energy eigenmode for a chain of length S + 1 = 51 is shown in Fig. 2a of the main text. Let us find the eigenstate of H BS that can simulate it. To this end, we will employ the following transformation where O = exp{ π 4 (ab † − a † b)}, which diagonalises H BS in the basis of Fock states Thus, the eigenstates of H diag BS are two-mode Fock states From the above we learn about the eigenstates of the original H BS Thus, the states |ψ (l) = O |l, S − l for l ∈ {0, . . . , S} are the eigenstates of H BS with corresponding eigenvalues of l − S 2 . In particular, the eigenstate defined by l = S 2 corresponds to the zero eigenvalue and thus, simulates the zero-energy mode The explicit form of this photonic state is as follows (0, S) are symmetric orthonormal Kravchuk functions. They may be expressed by means of the Gauss hypergeometric function The state |ψ (S/2) differs from the zero-energy eigenmode |Ψ 0 that we simulated and discussed in the main text by a phase factor. Nevertheless its photon statistics, which reads p(k) = φ (1/2) k (0, S) 2 and is shown in Fig. 10, is identical to these shown in Fig. 2 and 6 in the main text. For large S, the function 4 πS √ 1−(2k/(S−1)) 2 provides an asymptotic envelope to p(k), which is indicated as a solid curve in this figure.
Appendix D: Quantum program 2: simulation of perfect quantum state transfer Input data initialised to: |l = 0 a |S − l = S b . Program setting: we run this program for five different settings: r = 0.02, 0.15, 0.5, 0.85, and 0.98.

The Kravchuk transform
The α-fractional Kravchuk-Fourier transform of an input sequence {x k } for k = 0, 1, . . . , S is defined as follows where φ (p) k (l − Sp, S) is a Kravchuk function and p = sin 2 πα 4 . Mathematically, the beam-splitter interaction in the Fock state basis amounts to an α-fractional quantum Kravchuk transform (α-QKT) of the input state with fractionality α = 2θ π = 4 π arcsin √ r , where r is the beam splitter reflectivity [23]. In the supplementary material [23] we have provided analytical computations proving that the probability amplitude of detecting k and S − k photons behind the beam splitter provided that l and S − l were injected into it evaluates the Kravchuk transform 2-QKT is the spatial inversion operator. This becomes clear if we consider the matrix representation of a general beam-splitter interaction If we set r = sin 2 θ As α = 2 corresponds to r = 1, for ϕ = π 2 we obtain an inversion operation that swaps the input modes

Our simulation: perfect state transfer as a result of mirror reflection
The beam splitter interaction U  [45]. This leads to mapping {x l } to {X l } = {x S−l } thus, to the mirror reflection of the input sequence w.r.t. the centre of the domain. Since r = 1 corresponds to θ = π, U XY (t) performs the same operation at t = π, regardless the input state.
For any beam-splitter reflectivity, U BS corresponds to an α-QKT which is additive [23], U and c n c n+1 that don't conserve total particle number force us to include 'hole operators' c † n in the 'vector' χ. This doubles the dimension of our effective single particle Hamiltonian, and leads to the quasiparticles being a combination of particle and hole operators. The BdG Hamiltonian may be directly diagonalised to obtain the energy spectrum of the system.
The second quantized Hamiltonian can equally well be expressed in terms of Majorana operators a n = c n + c † n and b n = i(c † n − c n ) as {−µ n a n b n + (∆ n + t n )b n a n+1 + (∆ n − t n )a n b n+1 } = 1 2 Here (written for N = 3 for simplicity) are expressed in the basis of Majorana operators by the following unitary transfor- If we assign the values µ n = J 2n−1 , t n = ∆ n = J2n 2 , where J n is defined by Eq. (B4) and 2N = S + 1, H BdG is identical to our non-linear SSH and beam splitter Hamiltonian Eq. (B2). This generalised Kitaev chain of N sites thus has the same quasiparticle energy spectrum as the non-linear SSH chain of 2N sites i.e. l − S/2 for l ∈ {0, . . . , S}. The doubling of the energy spectrum is an artefact of the particle-hole symmetry imposed when constructing the Nambu spinor Ψ. A single site of the SSH chain is in effect mapped to a single Majorana fermion of the Kitaev chain, as can be seen in Fig. 11. The correspondence between the conventional SSH and Kitaev models is well known, being part of the broader equivalence between topological insulators and superconductors [46]. 2 may be simulated by our photonic system. Here the system is depicted in terms of a real Fermion operators and b Majorana operators. In this latter representation, the system shows a strong similarity to the generalised SSH model Fig. 9.

Interpreting the simulations
From the quantum simulation we can infer how an initial state evolves in this system. An initial state e.g. c † n |GS , will evolve after a time θ into the state U † K c † n U K |GS where U K = exp(−iθH K ) is the unitary operator for the second quantized Hamiltonian H K and |GS represents the ground state. Thus by taking combinations of the operators c n and c † n (or equivalently the Majorana operators a n and b n ) we can find the evolution of any state. Alternatively, one may find the evolution of these operators using the smaller BdG matrix U BdG (θ) = exp(−iθH BdG ). We simply write the operator using the Nambu basis (c 1 , c 2 , . . . c † 1 , c † 2 , . . .) and then operate on it with U BdG . For an example let's take N = 3 and find the evolution of c † 1 . The operator c † 1 is written in the Nambu basis as (0, 0, 0, 1, 0, 0), its evolution after a time θ = π is then determined by 3 as can be confirmed by direct calculation. This can also be done using the Majorana basis 1 √ 2 (a 1 , −ib 1 , . . . a N , −ib N ) and U BdG (θ) = exp(−iθH BdG ). The same calculation as above can then be performed, the operator c † 1 = 1 2 (a 1 − ib 1 ) is now written as (1/ √ 2, 1/ √ 2, 0, 0, 0, 0) and we evaluate as before. For our chosen parameters U BdG (θ) is identical to the beam splitter unitary operator where θ = 2 arcsin( √ r) and r is the beam splitter reflectivity. The above calculation then tells us how to perform the photonic simulation. To simulate a p-wave chain of N sites we take a photonic system of S = 2N − 1 photons and associate the two mode Fock states with the Majorana operators as To find how the operators evolves after a time θ we interfere the corresponding two mode Fock state on a beam splitter of reflectivity r = sin 2 (θ/2). To find how the real fermion operators transform, e.g. c † 1 = 1 2 (a 1 − ib 1 ), we can take superpositions U BS 1 √ 2 (|0 a |S b + |1 a |S − 1 b ). In a similar way, the evolution of multiply excited states e.g. U † K c † 1 c † 2 U K may be determined by decomposing into (U † K c † 1 U K )(U † K c † 2 U K ) and finding the evolution of the bracketed terms individually, which determines the evolution of the final state up to the global phase. One note of caution is that if we just perform photon number measurements we cannot distinguish between the states 1 √ 2 (|S − 1 a |1 b ± |S a |0 b ) and thus don't know if we obtain c 3 = 1 2 (a 3 +ib 3 ) or c † 3 = 1 2 (a 3 −ib 3 ). This problem may be easily resolved by extending the experimental detection to obtain more tomographically complete information, using e.g. homodyne detection.

Perfect transfer of Majorana fermions
Just like the generalised SSH model, this new system allows the perfect transfer of a quantum state due to its correspondence with the beam-splitter dynamics. After an interaction time θ = π an initial state localised at one end of the chain c † 1 |GS is perfectly transferred to the other end U † K c † 1 U K |GS = −ic † N |GS . This transfer can be viewed in terms of Majorana operators, in particular the two Majoranas at the ends of the chain swap after the interaction time θ = π reminiscent of a Majorana braiding operation [47]. If the system is left to evolve further these operators will keep evolving as, unlike the original Kitaev model, they are not zero energy modes. However, the above transfer could be viewed as an intermediate process, with the parameters being engineered from the original Kitaev values, µ n = 0, t n = ∆ n = constant, to the above values only during the transfer process. The physical engineering of such a system is extremely difficult in a condensed matter setting, but is possible in systems of 1-dimensional photonic cavities where the effective Majorana modes emerge due to Kerr-type non-linearities within a Bose-Hubbard model [27].
To simulate the evolution of a spin initially at the first site, we consider the initial state σ + 1 |↓ 1 , ↓ 2 , ↓ 3 . The final state after a time θ(r) is then found from interfering the two-mode state 1 √ 2 (|0 a |5 b − |1 a |4 b ) on a beam splitter of reflectivity r. Taking again θ = π (r = 1) as an example, we know that the perfect reflectivity implies that output state is 1 √ 2 (|4 a |1 b − |5 a |0 b ). Using the table Eq. (E14), we thus conclude that the final state of the spin system is σ z 1 σ z 2 σ + 3 |↓ 1 , ↓ 2 , ↓ 3 . The leading operators σ z n = 1 − 2σ − n σ + n only provide a phase factor and do not flip any spins so that the initial spin is transferred to the other end of the chain, just as we saw in section IV. The Ising model with uniform parameters has been previously studied in the context of quantum state transfer [48], our simulations suggest that modifying the parameters as described above could improve the situation, as the system inherits the perfect reflection property of the beam splitter when θ = π.

Relationship between simulated systems
We have discussed four example systems that may be simulated by our platform. The XY spin chain and non-linear SSH model are, strictly speaking, related by a Jordan-Wigner transformation, with the latter being expressed in terms of Fermionic operators in a crystal lattice. However, due to the conservation of total spin / particle number we can restrict ourselves to the single excitation subspace whereupon these two systems have identical interpretations and the terms are used interchangeably in the text. The other two systems, the non-uniform Ising and Kitaev models, are obtained from the XY spin chain and SSH models by a correspondence between topological insulating and topological superconducting systems. In simpler terms, we map the beam splitter Hamiltonian to Bogoliubovde Gennes Hamiltonians rather than single particle ones, which gives the simulations a slightly more complicated interpretation. We have depicted the relationships between these simulated systems in Fig. 12 for clarity.