Quantum simulations with multiphoton Fock states

TL;DR

This paper demonstrates the first quantum simulation using higher-order Fock states with multiple indistinguishable photons, enabling exploration of complex quantum phenomena beyond previous resources like qubits and coherent states.

Contribution

It introduces a novel approach to quantum simulation employing genuine multiphoton Fock states, expanding the capabilities of quantum modeling of topological matter and non-linear systems.

Findings

Revealed new topological phases

Simulated non-linear quantum systems

Demonstrated quantum transfer of Majorana fermions

Abstract

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.

Equations

$$H_{\text{BS}}=\frac{1}{2}(a^{†}b+a{b}^{†}),$$

$$p(k)=\bigl{|}\braket{k,S-k}{U_{\text{BS}}^{(r)}}{l,S-l}\bigr{|}^{2}=\bigl{|}\phi_{k}^{(r)}(l-Sr,S)\bigr{|}^{2},$$

$$S_x=\frac{a^{†}b+a{b}^{†}}{2},S_y=\frac{i\left(a{b}^{†}-a^{†}b\right)}{2},S_z=\frac{a^{†}a-b^{†}b}{2},S_0=\frac{a^{†}a+b^{†}b}{2},$$

$$[S_x,S_y]=i{S}_z,[S_y,S_z]=i{S}_x,[S_z,S_x]=i{S}_y.$$

$$H_{\text{BS}}=S_x.$$

$$S_z{|l⟩}_a{|S-l⟩}_b=(l-\frac{S}{2}){|l⟩}_a{|S-l⟩}_b,$$

$${|l⟩}_a{|S-l⟩}_b\equiv |\frac{S}{2};l-\frac{S}{2}⟩.$$

$$\boxed{{|l⟩}_a{|S-l⟩}_b\equiv |\frac{S}{2};l-\frac{S}{2}⟩\equiv |{\downarrow }_1,\dots ,{\uparrow }_{l+1},\dots ,{\downarrow }_{S+1}⟩={\sigma }_{l+1}^{+}|{\downarrow }_1,\dots ,{\downarrow }_{S+1}⟩,}$$

$$|m⟩={\sigma }_m^{+}|{\downarrow }_1,\dots ,{\downarrow }_{S+1}⟩.$$

$$H_{\text{XY}}=\sum _{n=1}^S\frac{J_n}{2}({\sigma }_n^x{\sigma }_{n+1}^x+{\sigma }_n^y{\sigma }_{n+1}^y)=\sum _{n=1}^S{J}_n\left({\sigma }_n^{+}{\sigma }_{n+1}^{-}+{\sigma }_n^{-}{\sigma }_{n+1}^{+}\right).$$

$${\left[{\mathbf{H}}_{\text{XY}}\right]}^{\text{Spin}}=$$

$$H_{\text{BS}}=\frac{1}{2}(a^{†}b+a{b}^{†}),$$

$$J_n=\frac{1}{2}\sqrt{n(S+1-n)},$$

$$\boxed{{\left[{\mathbf{H}}_{\text{BS}}\right]}_{nm}^{\text{Fock}}={\left[{\mathbf{H}}_{\text{XY}}\right]}_{nm}^{\text{Spin}}.}$$

$$\sum _{n=1}^{(S+1)/2}(-1{)}^n{e}^{-2n/\xi }{\sigma }_{2n}^{+}|{\downarrow }_1,\cdots {\downarrow }_{S+1}⟩,$$

$${\mathcal{O}}^{†}a\mathcal{O}=\frac{a-b}{\sqrt{2}},{\mathcal{O}}^{†}b\mathcal{O}=\frac{a+b}{\sqrt{2}},$$

$$H_{\text{BS}}^{\text{diag}}={\mathcal{O}}^{†}{H}_{\text{BS}}\mathcal{O}=\frac{1}{2}\left(a^{†}a-b^{†}b\right).$$

$$H_{\text{BS}}^{\text{diag}}|l,S-l⟩$$

$${\mathcal{O}}^{†}{H}_{\text{BS}}\mathcal{O}|l,S-l⟩$$

$$\displaystyle H_{\text{BS}}\bigl{(}\mathcal{O}\ket{l,S-l}\bigr{)}{}$$

$$H_{\text{BS}}\ket{\psi^{(S/2)}}=H_{\text{BS}}\bigl{(}\mathcal{O}\ket{\tfrac{S}{2},\tfrac{S}{2}}\bigr{)}=0.$$

$$\boxed{|{\psi }^{(S/2)}⟩=\mathcal{O}|\frac{S}{2},\frac{S}{2}⟩=\sum _{k=0}^S{e}^{-i\frac{\pi }{2}\frac{S}{2}}{\varphi }_k^{\left(1/2\right)}(0,S)|k,S-k⟩}$$

$${\varphi }_k^{(r)}(l-Sr,S)=(-1{)}^k\sqrt{\left(\genfrac{}{}{0ex}{}{S}{l}\right)\left(\genfrac{}{}{0ex}{}{S}{k}\right)}\sqrt{(1-r{)}^{S-l-k}{r}^{l+k}}{}_2{F}_1\left[-k,-l;-S;\frac{1}{r}\right].$$

$$X_k=$$

$$F_{k,l}^{\alpha }=$$

$${\mathcal{A}}_S(k,l)=⟨k,S-k|U_{\text{BS}}^{(r)}|l,S-l⟩=e^{i\frac{\pi }{2}(l-k-S\frac{\theta }{\pi })}{\varphi }_k^{\left(r\right)}(l-Sr,S).$$

$${\mathbf{U}}_{\text{BS}}=$$

$$H_{\text{K}}=\sum _{n=1}^N\left\{-{\mu }_n(c_n^{†}{c}_n-1/2)-t_n(c_{n+1}^{†}{c}_n+c_n^{†}{c}_{n+1})+{\Delta }_n(c_{n+1}^{†}{c}_n^{†}+c_n{c}_{n+1})\right\}.$$

$$H_{\text{K}}=\frac{1}{2}{\chi }^{†}{\mathbf{H}}_{\text{BdG}}\chi ,$$

Full text

I Introduction

Quantum simulations boost the development of topological materials Song2019 , quantum transport Harris2017 and quantum algorithms Jordan2012 for the benefit of low-power electronics Khang2018 , spintronics Smejkal2018 and quantum computing Nayak2008 . 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 Slussarenko2019 . 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 Houck2012 and electronic states of trapped ions Blatt2012 ; coherent states of photons Schreiber2012 and atoms (Bose–Einstein condensates) Fischer2004 ; and multiple single-particle Fock (number) states distributed among many modes in photonic waveguides Peruzzo2010 and optical lattices Bloch2012 . Thus, simulations have never seriously profited from interference of multi-particle Fock states, even though the importance of this regime has been recognised Flamini2019 , and the first attempt to mimic it with many-body systems was made Islam2015 .

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 Scheel2003 and also avoids error accumulation that weakens methods using quantum walks, built on numerous steps Preskill1998 . Our idea, shown in Fig. 1a-c, is based on overlapping two multiphoton Fock states, $\ket{l}_{a}$ and $\ket{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 $\text{su}(2)$.

Our platform also allows us to simulate certain classes of fermionic systems, e.g. a non-linear Su–Schrieffer–Heeger (SSH) model SSH1979 , 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) Kitaev2001 , and the transverse-field Ising model.

Due to recent advances in photon-number-resolved detection, we were able to employtransition-edge sensors (TESs) Gerrits2011 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 Hong1987 whose generalised form can be observed for higher-order Fock states if they are prepared in similar polarisation, spectral and spatio-temporal modes Campos1989 , as shown in Fig. 1d.

II Results

The Fock state quantum simulations build on a beam-splitter interaction $U_{\text{BS}}^{(r)}\!=\!e^{-i\theta(r)H_{\text{BS}}}$, guided by the Hamiltonian

[TABLE]

where $a^{\dagger}$ and $b^{\dagger}$ 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 $\theta(r)=2\arcsin\sqrt{r}$. For entries $\ket{l}_{a}$ and $\ket{S-l}_{b}$, the computational output from the beam-splitter and detectors is Stobinska2019

[TABLE]

where $\phi_{k}^{(r)}(x,S)$ is the Kravchuk function Atakishiyev1997 .

We selected three distinct examples of simulations, shown in Fig. 2, for experimental demonstration. The first a, uses input data initialised to $\ket{\tfrac{S}{2}}_{a}\ket{\tfrac{S}{2}}_{b}$ and the setting of $r=0.5$. For the second and third b & c, we set $\ket{0}_{a}\ket{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 becomes straightforward if we consider matrix representations of $H_{\text{BS}}$ and of the Hamiltonian describing a general chiral XY $\tfrac{1}{2}$-spin chain $H_{\text{XY}}=\sum_{n=1}^{S+1}\tfrac{J_{n}}{2}\,(\sigma^{x}_{n}\sigma^{x}_{n+1}+\sigma^{y}_{n}\sigma^{y}_{n+1})$, where $\sigma_{n}^{x}$ and $\sigma_{n}^{y}$ are the Pauli operators acting on the $n$th spin. In the single excitation subspace spanned by the states $\ket{n}=\sigma_{n}^{+}\ket{\downarrow_{1},\dots,\downarrow_{S+1}}$, where $\sigma_{n}^{+}=(1/2)\left(\sigma_{n}^{x}+i\sigma_{n}^{y}\right)$ is the raising operator, the latter has matrix elements $\left[\mathbf{H}_{\text{XY}}\right]^{\text{Spin}}_{mn}=\braket{m}{H_{\text{XY}}}{n}=J_{n-1}\delta_{n,m+1}+J_{m-1}\delta_{m,n+1}$, where $\delta_{i,j}$ denotes the Kronecker delta. The elements of $H_{\text{BS}}$ in the Fock state basis are given by $\left[\mathbf{H}_{\text{BS}}\right]^{\text{Fock}}_{nm}=\braket{n,S-n}{H_{\text{BS}}}{m,S-m}$. The two representations are identical, $\left[\mathbf{H}_{\text{BS}}\right]^{\text{Fock}}_{nm}=\left[\mathbf{H}_{\text{XY}}\right]^{\text{Spin}}_{nm}$, when we set the spin couplings to $J_{n}=\frac{1}{2}\sqrt{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 non-crystalline system. The existence of topologically non-trivial 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 $\mathbb{Z}$ topological invariant. Our first program performs a real-space study of this system and computes probabilities that describe its zero-energy eigenmode, $\lvert\Psi_{0}\rangle=\sum_{k=0}^{S}e^{-i\frac{\pi}{2}(S/2-k)}\,\phi_{k}^{(1/2)}(0,S)\,\sigma_{k+1}^{+}\ket{\downarrow_{1},\dots,\downarrow_{S+1}}$. Unlike 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 $\frac{4}{\pi S\sqrt{1-(2k/S-1)^{2}}}$, as outlined in Fig. 2a. The intensity-dependent amplitudes $J_{n}$ render $H_{\text{XY}}$ a generalisation of the seminal Su–Schrieffer–Heeger (SSH) model SSH1979 to the non-linear regime Gorlach2017 . 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 Christandl2004 . Our quantum simulation provides new insight into this system from which the perfect transfer becomes self-evident. The equivalence of $H_{\text{BS}}$ and $H_{\text{XY}}$ matrix representations implies the correspondence between interactions generated by these Hamiltonians, $U_{\text{BS}}^{(r)}$ and $U_{\text{XY}}(t)=e^{-itH_{\text{XY}}}$, respectively. Mathematically, the beam-splitter interaction in the Fock state basis amounts to an $\alpha$-fractional Quantum Kravchuk–Fourier transform ($\alpha$-QKT) of the input state with fractionality $\alpha=\frac{4}{\pi}\arcsin\sqrt{r}$ Stobinska2019 . As $2$-QKT is the spatial inversion operator Atakishiyev1997 , so is $U_{\text{XY}}(t)$ at $t=\pi$. 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 $\alpha=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 $\ket{0}_{a}\ket{S}_{b}$ is gradually transformed to $\ket{S}_{a}\ket{0}_{b}$ for increasing $r$, as shown in Fig. 2b.

Generalised Majorana modes. Multiphoton Fock state interference also facilitates the simulation of many-body systems that are not restricted to a single excitation subspace. For example, a p-wave superconducting chain (Kitaev model) Kitaev2001 is described by the mean field Hamiltonian $H_{\textrm{K}}=\sum_{n=1}^{N}\{-\mu_{n}(c_{n}^{\dagger}c_{n}-1/2)-t_{n}(c_{n+1}^{\dagger}c_{n}+c_{n}^{\dagger}c_{n+1})$ $+\Delta_{n}(c_{n+1}^{\dagger}c_{n}^{\dagger}+c_{n}c_{n+1})\}$, where $c_{n}^{\dagger}$ and $c_{n}$ are creation and annihilation operators for electrons on the $n$th atomic site, while $\mu_{n}$, $t_{n}$ and $\Delta_{n}$ are site dependent chemical potentials, hopping amplitudes and superconducting pairing potentials respectively. This Hamiltonian may be expressed in the form $H_{\textrm{K}}=\tfrac{1}{2}\chi^{\dagger}\mathbf{H}_{\textrm{BdG}}\chi$ where $\chi=\frac{1}{\sqrt{2}}(a_{1},-ib_{1},a_{2},-ib_{2},\ldots a_{N},-ib_{N})^{\textrm{T}}$ is a Nambu spinor and $\mathbf{H}_{\textrm{BdG}}$ is the Bogoliubov–de Gennes Hamiltonian matrix, in the basis of Majorana operators $a_{n}=c_{n}+c_{n}^{\dagger}$ and $b_{n}=i(c_{n}^{\dagger}-c_{n})$. The beam splitter Hamiltonian in the Fock state basis $\mathbf{H}_{\text{BS}}$ is identical to $\mathbf{H}_{\textrm{BdG}}$ for the parameters $\mu_{n}=J_{2n-1},t_{n}=\Delta_{n}=\tfrac{J_{2n}}{2}$, where $2N=S+1$. This correspondence allows one to simulate the Heisenberg evolution of the Majorana operators over the interaction time $\theta(r)$, as well as the evolution of the real fermion operators $c_{n}$ and $c_{n}^{\dagger}$, 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 $\ket{2(n-1)}_{a}\ket{S-2(n-1)}_{b}$ and $\ket{2n-1}_{a}\ket{S-(2n-1)}_{b}$ respectively. To evidence this, a further simulation with input $\ket{0}_{a}\ket{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=\tfrac{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 Bardyn2012 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_{\textrm{K}}=\frac{1}{2}\sum_{n=1}^{N}(\mu_{n}\sigma_{n}^{z}+2t_{n}\sigma_{n}^{x}\sigma_{n+1}^{x})$, since this is related to the p-wave superconducting chain by a Jordan–Wigner transformation. Due to the non-uniform field $\mu_{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 $\theta=\pi$. 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 $\frac{1}{\sqrt{2}}(\ket{0}_{a}\ket{S}_{b}+\ket{1}_{a}\ket{S-1}_{b})$ on a balanced beam splitter. See Supplementary Material for details.

III Experimental study

Fig. 3 shows the experimental integrated-photonics schema used for Fock state quantum simulations. Two pulsed spontaneous parametric down-conversion sources (SPDC) each generated independent two-mode photon-number-entangled states $\ket{\Psi}=\sum_{n=0}^{\infty}\sqrt{\braket{n}^{n}}\,\ket{n,n}$ with an average photon number $2\braket{n}=0.4$. For the pump repetition rate of $75\kern 2.5pt\mathrm{kHz}$ this led to approximately $0.46$ five-photon ($12$ four-photon) Fock states created per minute in each arm of the SPDC, of which about 0.2 (6) reached the detectors due to ca. $50\%$ losses in the set-up. One mode from each $\ket{\Psi}$ (the idlers, $c$ and $d$) was sent to a TES. Due to photon-number entanglement in $\ket{\Psi}$ states, the outcomes of TESs, $l$ and $S-l$, heralded the creation of Fock states $\ket{l}_{a}$ and $\ket{S-l}_{b}$ in the signal modes $a$ and $b$.

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_{\textrm{HOM}}=85.9\%$. Next, we took a measurement of the second order correlation function for each SPDC source separately and observed $g^{(2)}\geq 1.86\approx 1+V_{\textrm{HOM}}$, which corroborates the previous result. An effective Schmidt mode number of $K=\frac{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\times 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\times 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 TES2-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 zero-energy eigenmode of a non-linear SSH model described by $\left[\mathbf{H}_{\text{BS}}\right]^{\text{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 p-wave 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 Downing2019 . 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 Scheel2003 which can be harnessed in simulated models. The presented examples apply to a variety of systems such as superconducting nanowires Diez2012 , disordered graphene quasi-1D nanoribbons Han2007 and disordered cold atoms Pinheiro2015 . These may find applications in next-generation electronics McCaughan2014 and spintronics Nautiyal2004 operating with almost no energy dissipation and speeds exceeding 100GHz. 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 Stobinska2019 .

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 Zhong2018 , as well as in integrated platforms Wang2019 . 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 Peruzzo2010 . 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 Slussarenko2019 .

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 generation 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 Harder2016 . 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 $\sum_{l=0}^{S}x_{l}\ket{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 Thekkadath2020 . 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 universal. For example, Fock states are also available in motional states of trapped ions up to 10 excitations Ding2017 and in the form of plaquette Fock states of atoms in optical lattices up to 4 excitations Islam2015 . 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 Hofheinz2009 . This would further expand the assortment of input states that could be used for simulation and may give birth to new fascinating results.

V Methods

V.1 Characterisation of the set-up

Each integrated SPDC source produced a two-mode weakly squeezed vacuum state $\ket{\Psi}=\sum_{n=0}^{\infty}\lambda_{n}\,\ket{n,n}_{s,i}$, where $s$ and $i$ denote two output modes, named the signal and idler, $\lambda_{n}=\tfrac{\tanh^{n}g}{\cosh g}$, $\lvert\lambda_{n}\rvert^{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 $\ket{\Psi}$ is $\braket{n}=\sinh^{2}g$. The observed average photon number of $\braket{n}\approx 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\approx 1$ and thus, $\lambda_{n}\approx\sinh^{n}g=\sqrt{\braket{n}^{n}}$.

The transition-edge sensors (TESs) were operated at $70\kern 2.5pt\mathrm{mK}$, which allowed photon-number resolved measurements in all modes Gerrits2011 .

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)/no-photon’ 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 $\eta_{c}=50.3\%$ and $\eta_{d}=48.5\%$, respectively. The variable-coupler modes $a$ and $b$ exhibited a total efficiency of $\eta_{a}=21.6\%$ and $\eta_{b}=20.6\%$, respectively. These values result from the fact that each mode carried a $3\kern 2.5pt\mathrm{dB}$ loss from the coupler itself and another $1\kern 2.5pt\mathrm{dB}$ due to coupler insertion and fibre-to-fibre coupling losses. We estimated the transmission losses to be approximately of $50\%\approx 3\kern 2.5pt\mathrm{dB}$. Here $1\kern 2.5pt\mathrm{dB}$ stands for the initial fibre in-coupling loss due to spatial mode mismatch, while $0.25\kern 2.5pt\mathrm{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.

The HOM visibility is computed using the formula $v^{(2)}=\frac{n_{\text{max}}-n_{\text{min}}}{n_{\text{max}}+n_{\text{min}}}$, where $n_{\text{max}}$ and $n_{\text{min}}$ are the maximal and minimal number of events registered by the TES detectors for the given photon number $S$. In the experiment for input $\ket{2,2}$ and $r=0.5$, we obtained $v^{(2)}=50.6\%\pm 1.2\%$, whereas for input $\ket{0,4}$ ($\ket{0,5}$) they were $99.1\%\pm 2.5\%$ ($97.8\%\pm 6.2\%$) for $r=0.04$, $87.6\%\pm 2.2\%$ ($96.7\%\pm 7.2\%$) for $r=0.3$, $65.7\%\pm 1.7\%$ ($71.4\%\pm 4.6\%$) for $r=0.5$ and $99.9\%\pm 0.8\%$ ($98.6\%\pm 7.2\%$) for $r=0.96$.

V.2 Error estimation

In the experiment, each measurement results in a 4-tuple consisting of the number of photons registered by TES1-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 $\Delta p=1/\sqrt{N}$.

V.3 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^{\dagger}$ with a superposition of the same mode $b^{\dagger}$ and an orthogonal one $b^{\dagger}_{\perp}$, i.e. $b^{\dagger}\to\cos y\,b^{\dagger}+\sin y\,b^{\dagger}_{\perp}$, where the parameter $y\in(0,\tfrac{\pi}{2})$ introduced weights and ‘tuned’ the distinguishability. This transformation led to the interference of $\ket{l}_{a}$ with a two-mode Fock state superposition$\sum_{n=0}^{S-l}\binom{S-l}{n}^{-1/2}\cos^{n}y\,(\sin y)^{S-l-n}\ket{n}_{b}\ket{S-l-n}_{b\perp}$ instead of the single-mode Fock state$\ket{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\sqrt{(K-1)/K}$, where $K$ denoted the effective Schmidt mode number measured during the set-up characterisation. For $K=1.16$, we used $y=0.38$.

Realistic model of photodetection requires taking into account a probability of detecting $n_{d}$ photons when a Fock state $\ket{n_{\text{in}}}$ reaches a TES. It is given by $p_{\text{TES}}(n_{\text{in}},n_{d})=\binom{n_{\text{in}}}{n_{d}}(1-\eta)^{n_{\text{in}}-n_{d}}\,\eta^{n_{d}}$ where $n_{d}\leq n_{\text{in}}$ and $\eta$ is the detector efficiency. In our computations we first used a starting value of $\eta=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.

Acknowledgements.

Funding:

TS, TM, AB and MS were supported by the Foundation for Polish Science ‘First Team’ project No. POIR.04.04.00-00-220E/16-00 (originally: FIRST TEAM/2016-2/17) and the National Science Centre ‘Sonata Bis’ project No. 2019/34/E/ST2/00273. AE and IW were supported by the Engineering and Physical Sciences Research Council project No. EP/K034480/1.

Authors contributions:

TS, TM, AB and MS developed the theory while AE, WRC, WSK, and IAW were responsible for realisation of the experiment. 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 $\text{su}(2)$ can be represented in terms of annihilation and creation operators of a harmonic oscillator. This is known as the Schwinger representation Chruscinski2004 . It allows one to associate two independent quantum-harmonic oscillator modes with spin operators as follows

[TABLE]

where $S_{0}$ is the Casimir operator $S_{0}(S_{0}+1)=S_{x}^{2}+S_{y}^{2}+S_{z}^{2}$ and the standard $\text{su}(2)$ commutation relations hold

[TABLE]

Therefore, we can immediately identify that

[TABLE]

In this picture, a two-mode Fock state $\ket{l}_{a}\ket{S-l}_{b}$ corresponds to spin-$\tfrac{S}{2}$ particle that is prepared in an eigenstate of $S_{z}$ with eigenvalue $l-\tfrac{S}{2}$, known as a Dicke state $\ket{\tfrac{S}{2};l-\tfrac{S}{2}}$

[TABLE]

[TABLE]

Furthermore, one can one-to-one map the Dicke states $\ket{\tfrac{S}{2};m}$ to the basis states that span the single excitation subspace of a spin-$\tfrac{1}{2}$ chain. To this end, we employ the following relabelling $m=-\tfrac{S}{2}+n-1$, where $1\leq n\leq S+1$ denotes $n$th spin-$\tfrac{1}{2}$ in the chain with $S+1$ sites. Then, $\ket{\tfrac{S}{2};l-\tfrac{S}{2}}$ corresponds to the chain where $(l\!+\!1)$th spin-$\tfrac{1}{2}$ is excited

[TABLE]

where $\sigma_{m}^{+}=\frac{1}{2}\left(\sigma_{m}^{x}+i\sigma_{m}^{y}\right)$ is the raising operator acting on $m$th spin, with $\sigma_{m}^{x}$ and $\sigma_{m}^{y}$ denoting the Pauli operators. For a concise notation, we denote such spin-chain states as follows

[TABLE]

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

[TABLE]

In the single excitation subspace that is spanned by the states shown in Eq. (9), this Hamiltonian has the matrix representation

[TABLE]

where the matrix elements equal $\left[\mathbf{H}_{\text{XY}}\right]^{\text{Spin}}_{mn}=\bra{m}H_{\text{XY}}\ket{n}=J_{n-1}\,\delta_{n,m+1}+J_{m-1}\,\delta_{m,n+1}$ and $\delta_{i,j}$ is the Kronecker delta; $\delta_{i,j}=1$ if $i=j$ and $\delta_{i,j}=0$ otherwise.

The beam splitter Hamiltonian Kim2002

[TABLE]

where $a^{\dagger}$ ($a$) and $b^{\dagger}$ ($b$) denote photonic creation (annihilation) operators which act on the interferometer input modes, features the following matrix representation elements in the Fock state basis $\left[\mathbf{H}_{\text{BS}}\right]^{\text{Fock}}_{nm}=\bra{n,S-n}H_{\text{BS}}\ket{m,S-m}=\frac{1}{2}\sqrt{n(S+1-n)}\,\delta_{n,m+1}+\frac{1}{2}\sqrt{m(S+1-m)}\,\delta_{m,n+1}$. Please note that in the notation used the photonic state $\ket{m,S-m}$ corresponds to spin chain state $\ket{m+1}$. Therefore, if we set the spin-couplings to

[TABLE]

then

[TABLE]

Appendix C Quantum program 1: simulation of weakly localised edge states

Input data initialised to: $\ket{l=\tfrac{S}{2}}_{a}\ket{S-l=\tfrac{S}{2}}_{b}$.

Program setting: $r=0.5$.

C.1 Short introduction to edge states: the Su–Schrieffer–Heeger (SSH) model

The SSH model is the seminal example of a 1D system where edge (topologically non-trivial) states can be observed. The system is a dimerised chain of atoms that is described by a Hamiltonian of the form shown in Eq. (10) with periodic couplings of $J_{n}=2J(1+\delta(-1)^{n})$, where $-1\leq\delta\leq 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 ($\delta>0$) or the stronger coupling ($\delta<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 $\delta=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

[TABLE]

where \xi=2/\ln\bigl{(}\frac{1+\delta}{1-\delta}\bigr{)} is the localisation length Nevado . Usually, the boundary states are studied for $\delta$ 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 $\delta=0.005$ (blue). Interestingly, the conventional SSH model can present weakly localised edge states that are nonetheless topological.

C.2 Edge states in a non-linear SSH model

An interesting generalisation of the SSH model to the non-linear domain was achieved by setting intensity-dependent site couplings Gorlach2017 . This model was theoretically implemented with an array of cavities with the tunneling constants equal $J_{2n}=J_{2}+\alpha\left(\lvert E_{2n}\rvert^{2}+\lvert E_{2n+1}\rvert^{2}\right)$ and $J_{2n+1}=J_{1}=\textrm{const}$, as shown in Fig. 7. Here $E_{n}$ denotes the field amplitude in $n$th resonator. Fig. 8 depicts a self-induced topological edge state that arises in this system. It is plotted for $J_{1}>J_{2}$, a regime where the linear SSH model shows no boundary states. Interestingly, its envelope reveals no exponential decay.

C.3 Our simulation: edge states in a generalised non-linear SSH model

Interference of Fock states on a beam splitter can simulate an interacting spin chain with non-periodic next-neighbour spin couplings shown in Eq. (13) and Hamiltonian given by Eq. (10). This system is depicted in Fig. 9. Interestingly, the couplings (13) also are intensity dependent, as $n$ is the intensity of Fock state $\ket{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 Kitaev2009 , a system is categorised according to the properties of its time-reversal operator $T=U_{T}\otimes\mathcal{K}$, charge-conjugation operator $C=U_{C}\otimes\mathcal{K}$, and chiral-symmetry operator $\Gamma=T\otimes C$, where $U_{T(C)}$ is a unitary operator and $\mathcal{K}$ is complex conjugation. If there exists a $T$ ($C$ or $\Gamma$) 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 beam-splitter Hamiltonian matrix representation $\left[\mathbf{H}_{\text{BS}}\right]^{\text{Fock}}$ is real-valued and thus, its time-reversal symmetry operator is simply $T=1\otimes\mathcal{K}$. Due to the absence of couplings beyond nearest-neighbour the chiral symmetry operator is given by $\Gamma_{ij}=(-\delta_{ij})^{2}$, and finally $C=T\otimes\Gamma$. Our photonic system possesses all three symmetries with $(T^{2},C^{2},\Gamma^{2})=(1,1,1)$. Thus it belongs to the BDI (chiral orthogonal) class of topological insulators which in one-dimension is characterised by a $\mathbb{Z}$ topological invariant. Therefore, the simulated spin chain does so too.

Weakly localised zero-energy edge states

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_{\text{BS}}$ that can simulate it. To this end, we will employ the following transformation

[TABLE]

where $\mathcal{O}=\exp\{\tfrac{\pi}{4}(ab^{\dagger}-a^{\dagger}b)\}$, which diagonalises $H_{\text{BS}}$ in the basis of Fock states

[TABLE]

Thus, the eigenstates of $H_{\text{BS}}^{\text{diag}}$ are two-mode Fock states

[TABLE]

From the above we learn about the eigenstates of the original $H_{\text{BS}}$

[TABLE]

Thus, the states $\ket{\psi^{(l)}}=\mathcal{O}\ket{l,S-l}$ for $l\in\{0,\dots,S\}$ are the eigenstates of $H_{\text{BS}}$ with corresponding eigenvalues of $l-\frac{S}{2}$. In particular, the eigenstate defined by $l=\tfrac{S}{2}$ corresponds to the zero eigenvalue and thus, simulates the zero-energy mode

[TABLE]

The explicit form of this photonic state is as follows

[TABLE]

where $\phi_{k}^{\left(1/2\right)}(0,S)$ are symmetric orthonormal Kravchuk functions. They may be expressed by means of the Gauss hypergeometric function

[TABLE]

The state $\lvert\psi^{(S/2)}\rangle$ differs from the zero-energy eigenmode $\lvert\Psi_{0}\rangle$ that we simulated and discussed in the main text by a phase factor. Nevertheless its photon statistics, which reads p(k)=\bigl{\lvert}\phi_{k}^{(1/2)}(0,S)\bigr{\rvert}^{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 $\frac{4}{\pi S\sqrt{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: $\ket{l=0}_{a}\ket{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$.

D.1 The Kravchuk transform

The $\alpha$-fractional Kravchuk–Fourier transform of an input sequence $\{x_{k}\}$ for $k=0,1,\ldots,S$ is defined as follows

[TABLE]

where $\phi_{k}^{(p)}(l-Sp,S)$ is a Kravchuk function and p=\sin^{2}\bigl{(}\frac{\pi\alpha}{4}\bigr{)}.

Mathematically, the beam-splitter interaction in the Fock state basis amounts to an $\alpha$-fractional quantum Kravchuk transform ($\alpha$-QKT) of the input state with fractionality $\alpha=\tfrac{2\theta}{\pi}=\frac{4}{\pi}\arcsin\sqrt{r}$, where $r$ is the beam splitter reflectivity Stobinska2019 . In the supplementary material Stobinska2019 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

[TABLE]

2-QKT is the spatial inversion operator. This becomes clear if we consider the matrix representation of a general beam-splitter interaction

[TABLE]

If we set $r=\sin^{2}\tfrac{\theta}{2}$

[TABLE]

As $\alpha=2$ corresponds to $r=1$, for $\varphi=\tfrac{\pi}{2}$ we obtain an inversion operation that swaps the input modes

[TABLE]

D.2 Our simulation: perfect state transfer as a result of mirror reflection

The beam splitter interaction $U_{\text{BS}}^{(r)}=e^{-i\theta(r)H_{\text{BS}}}$ can simulate the dynamics $U_{\text{XY}}(t)=e^{-itH_{\text{XY}}}$ of the spin chain (10) in its single-excitation subspace. $U_{\text{BS}}^{(1)}$ performs a spatial inversion operation of the sequence $\{x_{l}\}$, where $l=0,\ldots,S$, for input state $\ket{\Psi}=\sum_{l=0}^{S}x_{l}\ket{l}\ket{S-l}$ with respect to the point $l=\frac{S}{2}$ Hakioglu . 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 $\theta=\pi$, $U_{\text{XY}}(t)$ performs the same operation at $t=\pi$, regardless the input state.

For any beam-splitter reflectivity, $U_{\text{BS}}^{(r)}$ corresponds to an $\alpha$-QKT which is additive Stobinska2019 , $U_{\text{BS}}^{(r(\theta_{1}+\theta_{2}))}=U_{\text{BS}}^{(r(\theta_{1}))}U_{\text{BS}}^{(r(\theta_{2}))}$, where $r(\theta)=\sin^{2}\frac{\theta}{2}$. Thus, $U_{\text{BS}}^{(1)}$ can be decomposed into $N$ infinitesimal evolutions $U_{\text{BS}}^{(1)}=\bigotimes_{i=1}^{N}U_{\text{BS}}^{(r_{i})}$, where $\sum_{i=1}^{N}\theta(r_{i})=\theta(1)=\pi$. This property is demonstrated by our quantum simulations performed for subsequent values of $r$, shown in Figs. 2b and 4b in the main text.

Appendix E Additional simulations: Majorana modes and the Ising model

E.1 Generalised Kitaev model

The quantum simulations based on Fock-state interference may be reinterpreted in the language of Bogoliubov–de Gennes Hamiltonians to simulate systems that are not restricted to a single excitation subspace.

A one-dimensional p-wave superconducting chain of $N$ atomic sites is described by the following second quantized Hamiltonian

[TABLE]

This is the Kitaev model Kitaev2001 but with site dependent chemical potential $\mu_{n}$, hopping amplitudes $t_{n}$ and energy gap $\Delta_{n}$. This Hamiltonian may be diagonalised using the Bogoliubov–de Gennes trick as follows. First, each term is written in a symmetric form using the fermion anticommutation relations, for example $c_{n}^{\dagger}c_{n}=\frac{1}{2}(c_{n}^{\dagger}c_{n}-c_{n}c_{n}^{\dagger})+\frac{1}{2}$. Secondly, we introduce the Nambu spinor of all annihilation and creation operators $\chi=(c_{1},c_{2},\ldots c_{1}^{\dagger},c_{2}^{\dagger},\ldots)^{\textrm{T}}$. Finally, we write the Hamiltonian as

[TABLE]

where the Bogoliubov–de Gennes (BdG) Hamiltonian $\mathbf{H}_{\textrm{BdG}}$ is a $2N\times 2N$ matrix that may be interpreted as an effective single particle Hamiltonian for the systems quasiparticles. Notice that the existence of terms like $c_{n+1}^{\dagger}c_{n}^{\dagger}$ and $c_{n}c_{n+1}$ that don’t conserve total particle number force us to include ‘hole operators’ $c_{n}^{\dagger}$ in the ‘vector’ $\chi$. 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}^{\dagger}$ and $b_{n}=i(c_{n}^{\dagger}-c_{n})$ as

[TABLE]

Here $\chi^{\prime}=\frac{1}{\sqrt{2}}(a_{1},-ib_{1},a_{2},-ib_{2},\ldots a_{N},-ib_{N})^{\textrm{T}}$ and

[TABLE]

(written for $N=3$ for simplicity) are expressed in the basis of Majorana operators by the following unitary transformation $\chi^{\prime}=\mathbf{M}\chi$, $\mathbf{H}_{\textrm{BdG}}^{\prime}=\mathbf{M}\mathbf{H}_{\textrm{BdG}}\mathbf{M}^{\dagger}$ where

[TABLE]

If we assign the values $\mu_{n}=J_{2n-1},t_{n}=\Delta_{n}=\tfrac{J_{2n}}{2}$, where $J_{n}$ is defined by Eq. (13) and $2N=S+1$, $\mathbf{H}_{\textrm{BdG}}^{\prime}$ is identical to our non-linear SSH and beam splitter Hamiltonian Eq. (11). 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\in\{0,\dots,S\}$. The doubling of the energy spectrum is an artefact of the particle-hole symmetry imposed when constructing the Nambu spinor $\Psi$. 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 Cobanera2015 .

E.2 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}^{\dagger}\ket{\textrm{GS}}$, will evolve after a time $\theta$ into the state $U_{\textrm{K}}^{\dagger}c_{n}^{\dagger}U_{\textrm{K}}\ket{\textrm{GS}}$ where $U_{\textrm{K}}=\exp(-i\theta H_{\textrm{K}})$ is the unitary operator for the second quantized Hamiltonian $H_{\textrm{K}}$ and $\ket{\textrm{GS}}$ represents the ground state. Thus by taking combinations of the operators $c_{n}$ and $c_{n}^{\dagger}$ (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_{\textrm{BdG}}(\theta)=\exp(-i\theta H_{\textrm{BdG}})$. We simply write the operator using the Nambu basis $(c_{1},c_{2},\ldots c_{1}^{\dagger},c_{2}^{\dagger},\ldots)$ and then operate on it with $U_{\textrm{BdG}}$. For an example let’s take $N=3$ and find the evolution of $c_{1}^{\dagger}$. The operator $c_{1}^{\dagger}$ is written in the Nambu basis as $(0,0,0,1,0,0)$, its evolution after a time $\theta=\pi$ is then determined by

[TABLE]

indicating that $c_{1}^{\dagger}$ evolves into $U_{\textrm{K}}^{\dagger}c_{1}^{\dagger}U_{\textrm{K}}=-ic_{3}^{\dagger}$ as can be confirmed by direct calculation.

This can also be done using the Majorana basis $\frac{1}{\sqrt{2}}(a_{1},-ib_{1},\ldots a_{N},-ib_{N})$ and $U_{\textrm{BdG}}^{\prime}(\theta)=\exp(-i\theta H_{\textrm{BdG}}^{\prime})$. The same calculation as above can then be performed, the operator $c_{1}^{\dagger}=\frac{1}{2}(a_{1}-ib_{1})$ is now written as $(1/\sqrt{2},1/\sqrt{2},0,0,0,0)$ and we evaluate

[TABLE]

indicating that $c_{1}^{\dagger}=\frac{1}{2}(a_{1}-ib_{1})$ evolves into $\frac{-i}{2}(a_{3}-ib_{3})=-ic_{3}^{\dagger}$ as before.

For our chosen parameters $U_{\textrm{BdG}}^{\prime}(\theta)$ is identical to the beam splitter unitary operator where $\theta=2\arcsin(\sqrt{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

[TABLE]

To find how the operators evolves after a time $\theta$ we interfere the corresponding two mode Fock state on a beam splitter of reflectivity $r=\sin^{2}(\theta/2)$. To find how the real fermion operators transform, e.g. $c_{1}^{\dagger}=\frac{1}{2}(a_{1}-ib_{1})$, we can take superpositions $U_{\textrm{BS}}\frac{1}{\sqrt{2}}(\ket{0}_{a}\ket{S}_{b}+\ket{1}_{a}\ket{S-1}_{b})$. In a similar way, the evolution of multiply excited states e.g. $U_{\textrm{K}}^{\dagger}c_{1}^{\dagger}c_{2}^{\dagger}U_{\textrm{K}}$ may be determined by decomposing into $(U_{\textrm{K}}^{\dagger}c_{1}^{\dagger}U_{\textrm{K}})(U_{\textrm{K}}^{\dagger}c_{2}^{\dagger}U_{\textrm{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 $\frac{1}{\sqrt{2}}(\ket{S-1}_{a}\ket{1}_{b}\pm\ket{S}_{a}\ket{0}_{b})$ and thus don’t know if we obtain $c_{3}=\frac{1}{2}(a_{3}+ib_{3})$ or $c_{3}^{\dagger}=\frac{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 $\theta=\pi$ an initial state localised at one end of the chain $c_{1}^{\dagger}\ket{\textrm{GS}}$ is perfectly transferred to the other end $U_{\textrm{K}}^{\dagger}c_{1}^{\dagger}U_{\textrm{K}}\ket{\textrm{GS}}=-ic_{N}^{\dagger}\ket{\textrm{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 $\theta=\pi$

[TABLE]

reminiscent of a Majorana braiding operation Beenakker2013 . 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, $\mu_{n}=0$, $t_{n}=\Delta_{n}=\textrm{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 Bardyn2012 .

E.3 Transverse-field Ising Model

The Hamiltonian Eq. (30) may be mapped to a spin chain using the Jordan–Wigner transformation

[TABLE]

where $\sigma_{m}^{z}$ acts on the $m$th spin of the chain and similarly for the other Pauli operators. The resulting Hamiltonian is

[TABLE]

For the couplings we consider $\mu_{n}=J_{2n-1},t_{n}=\Delta_{n}=\tfrac{J_{2n}}{2}$, this becomes

[TABLE]

which is an Ising chain in a non-uniform transverse field. Thus in principle we can simulate this system using our Fock-state interference platform. One inverts the transformations Eq. (39) and then uses the correspondence Eq. (37) to associate the two-mode Fock states with combinations of spin operators. For a concrete example we again take $N=3$ $(S=2N-1=5)$, the Jordan–Wigner transformation is then

[TABLE]

while, using Eq. (37) and $c_{n}=\frac{1}{2}(a_{n}+ib_{n})$, the mapping to Fock states is

[TABLE]

To simulate the evolution of a spin initially at the first site, we consider the initial state $\sigma_{1}^{+}\ket{\downarrow_{1},\downarrow_{2},\downarrow_{3}}$. The final state after a time $\theta(r)$ is then found from interfering the two-mode state $\frac{1}{\sqrt{2}}(\ket{0}_{a}\ket{5}_{b}-\ket{1}_{a}\ket{4}_{b})$ on a beam splitter of reflectivity $r$. Taking again $\theta=\pi$ $(r=1)$ as an example, we know that the perfect reflectivity implies that output state is $\frac{1}{\sqrt{2}}(\ket{4}_{a}\ket{1}_{b}-\ket{5}_{a}\ket{0}_{b})$. Using the table Eq. (43), we thus conclude that the final state of the spin system is $\sigma_{1}^{z}\sigma_{2}^{z}\sigma_{3}^{+}\ket{\downarrow_{1},\downarrow_{2},\downarrow_{3}}$. The leading operators $\sigma_{n}^{z}=1-2\sigma_{n}^{-}\sigma_{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 Yao2013 , 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 $\theta=\pi$.

E.4 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 Bogoliubov–de 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.

Appendix F Extension to multiple dimensions

As an example of how to extend our simulations to multiple dimensions, we take the XY spin model (or equivalently, generalised SSH model) in a two dimensional rectangular lattice. The Hamiltonian for a lattice of size $(S_{x}+1)\times(S_{y}+1)$ is

[TABLE]

Here the operators $\sigma_{m,n}^{+}$ etc. act on the spin at site $(m,n)$, the terms $J_{n}^{(S_{x})}$ denote couplings along the x-axis of the lattice, while $J_{n}^{(S_{y})}$ denote couplings along the y-axis. When $S_{y}=0$, this reduces to the 1-dimensional case Eq. (10). The system is depicted in Fig. (13). Just as in the 1-dimensional case we look at the Hamiltonian matrix elements in the single excitation subspace. Denoting the state $\ket{i,j}=\sigma_{i,j}^{+}\ket{\downarrow_{1,1}\ldots\downarrow_{S_{x}+1,S_{y}+1}}$, these are

[TABLE]

To perform the simulation, we take four photonic modes which we label $a,a^{\prime},b,b^{\prime}$. The modes $a$ and $b$ are interfered on a beam splitter, and similarly with $a^{\prime}$ and $b^{\prime}$. This could be achieved with two different beam splitters, or the same beam splitter in different frequency, polarisation or spatio-temporal modes. The total Hamiltonian is $H_{BS}=\frac{1}{2}(a^{\dagger}b+ab^{\dagger}+a^{\prime\dagger}b^{\prime}+a^{\prime}b^{\prime\dagger})$ and in the Fock state basis $\ket{i}_{a}\ket{S_{x}-i}_{b}\ket{j}_{a^{\prime}}\ket{S_{y}-j}_{b^{\prime}}$ one can check that the matrix elements are identical to Eq. (45) provided that $J_{n}^{(S)}=\tfrac{1}{2}\sqrt{n(S+1-n)}$. To simulate a system of size $(S_{x}+1)\times(S_{y}+1)$ one simply interferes the corresponding Fock states, with $S_{x}$ total photons in modes $a$ and $b$, and $S_{y}$ total photons in modes $a^{\prime}$ and $b^{\prime}$.