Photonic quantum simulations of coupled $PT$-symmetric Hamiltonians

Parity-time ($PT$) symmetric Hamiltonians are generally non-Hermitian and give rise to exotic behaviour in quantum systems at exceptional points, where eigenvectors coalesce. The recent realisation of $PT$-symmetric Hamiltonians in quantum systems has ignited efforts to simulate and investigate many-particle quantum systems across exceptional points. Here we use a programmable integrated photonic chip to simulate a model comprised of twin pairs of $PT$-symmetric Hamiltonians, with each the time reverse of its twin. We simulate quantum dynamics across exceptional points including two- and three-particle interference, and a particle-trembling behaviour that arises due to interference between subsystems undergoing time-reversed evolutions. These results show how programmable quantum simulators can be used to investigate foundational questions in quantum mechanics.

Parity-time (PT ) symmetric Hamiltonians are generally non-Hermitian and give rise to exotic behaviour in quantum systems at exceptional points, where eigenvectors coalesce. The recent realisation of PT -symmetric Hamiltonians in quantum systems has ignited efforts to simulate and investigate many-particle quantum systems across exceptional points. Here we use a programmable integrated photonic chip to simulate a model comprised of twin pairs of PT -symmetric Hamiltonians, with each the time reverse of its twin. We simulate quantum dynamics across exceptional points including two-and three-particle interference, and a particle-trembling behaviour that arises due to interference between subsystems undergoing time-reversed evolutions. These results show how programmable quantum simulators can be used to investigate foundational questions in quantum mechanics.
Dirac Hermiticity of a Hamiltonian has been a tenet of quantum theory since its inception. This constraint guarantees real eigenvalues, orthogonal eigenstates, and a unitary time evolution; it reflects the dynamics of an isolated system. Over the past two decades, non-Hermitian Hamiltonians that are symmetric under combined parity (P) and time-reversal (T ) transformations were extensively investigated, first mathematically [1,2] and then experimentally in classical wave systems [3,4]. In the classical domain, PT -symmetric systems have shown remarkable properties such as unidirectional invisibility [5], single-mode [6,7] and topological [8,9] lasing, topological energy transfer [10], mode switching [11], and enhanced sensitivity [12,13]. While generally non-Hermitian, PT -symmetric Hamiltonians exhibit exceptional points (EPs) at which two (or more) eigenvalues become degenerate and the corresponding eigenstates coalesce. For all these advances, the role of thermal and quantum noise on EPs in these semi-classical systems is not yet understood [14,15].
Realizing PT -symmetric systems in the quantum domain has been a major challenge [16], and its circumvention has relied on two methods. The first approach, known as passive PT symmetry, is based on the equivalence (after normalisation) of a PT -symmetric Hamiltonian and a dissipative Hamiltonian with mode-selective losses [17,18]. Such a dissipative Hamiltonian arises naturally from a Lindblad approach for open quantum systems when certain quantum jumps are ignored and allows the implementation of coherent, non-unitary dynamics in a minimal quantum system [19]. The second method uses Hamiltonian dilation that places the PTsymmetric Hamiltonian in a larger Hilbert space by introducing an ancillary two-level system [20][21][22], and the coherent, non-unitary dynamics of interest are recovered from the unitary dynamics in the larger Hilbert space. In the first approach, the evolution time of the system is limited by the exponentially decaying signal or postselection success probability, while the second method requires a time-dependent Hermitian Hamiltonian for the larger Hilbert space.
Here we propose a framework for the quantum simulation of PT -symmetric Hamiltonians that is particularly suited to technology platforms that allow the direct implementation of unitary transformations [23][24][25]. The non-unitary evolution operator, and a second operator that is the time reverse of the first, are embedded together into a global unitary transformation that is then implemented by the device. The overall evolution in this model allows single-or many-particle excitations to tunnel between the twin systems, with a probability that scales with the non-Hermicity of the simulated Hamiltonians. This construction permits the experimental investigation of states superposed across opposite directions of time in the context of non-Hermitian Hamiltonians.
We experimentally simulate multiparticle dynamics [26] in two-and three-mode PT -symmetric Hamiltonians using a programmable photonic chip [27],and ensembles of one, two, and three photon input states. We reproduce the dynamics in the PT -symmetric regime (defined by real eigenvalues) and across the EP into the PT -symmetry broken regime (defined by complex conju-gate eigenvalues), including the effects of mutual coherence and interference between the forward-in-time and backward-in-time subsystems. Combined with the current advances in photonic technologies, these demonstrations point towards the inclusion of non-Hermitian simulations into the realm of quantum technologies.

Model and simulation procedure
Our simulation approach begins by considering a family of PT -symmetric systems consisting of a linear chain of N modes with a coupling amplitude J > 0 between nearest-neighbour modes. In addition, an imaginary potential of opposite sign ±iγ located in the two spatially symmetric end modes determines a gain and a loss site. We fix the energy scale of our system by choosing J = +1 so that the Hamiltonian reads where |k is the state associated to the k-th mode. Using {|k } N k=1 as a basis, the parity operator we choose is represented by the antidiagonal matrix P mn = δ m,N +1−n and the time-reversal operator T corresponds to the complex conjugation operation. H N (γ) presents an EP for a critical value γ = γ c that depends on N [28,29]. The spectrum of the Hamiltonian H N (γ) is purely real for γ < γ c and the corresponding eigenvectors are simultaneous eigenvectors of the antilinear PT operator (unbroken symmetry phase). At the EP, γ = γ c , two or more states coalesce, while at γ > γ c , two eigenstates with opposite imaginary eigenvalue emerge from the previously coalescing states (broken symmetry phase). The full spectrum of the Hamiltonian can be obtained by solving a set of transcendental equations reported in [28] but γ c is given by a simple closed form: γ c /J = 1 whenever N is even and γ c /J = (N + 1)/(N − 1) when N is odd. For a static, non-Hermitian Hamiltonian H N (γ), the non-unitary time evolution operator G N (γ, t) is given by (assuming = 1) We note that these choices set the time-scale for our system, /J, to unity. To simulate its dynamical evolution with a device that implements unitary transformations, we adopt the Halmos unitary dilation [30,31]. Namely, we first rescale G N (γ, t) by its largest singular value ||G N (γ, t)|| 2 obtain-ingG so that the operator norm of the scaled operatorG N (γ, t) is unity. Then, defining the defect operator as we construct the following unitary transformation U 2N (γ, t) defined on a Hilbert space twice the size. The block matrix representation of U 2N (γ, t) is We map the modes of the dilated Hilbert space onto the spatial modes of a reprogrammable linear optical chip. The time evolution of excitations in the PT -symmetric system is simulated by injecting single photons into the optical circuit and reprogramming it to implement the unitary transfer matrix U 2N (γ, t) for a sequence of time steps {U 2N (γ, t i )} i , see Fig. 1. To simulate the evolution determined by H N (γ), we encode the initial state in the first N modes of the interferometer and consider detection events on the outputs of those same modes. Our construction is applicable for PT -symmetric systems even when the spectrum of the Hamiltonian becomes complex, allowing simulations across exceptional points. We also emphasize that since we realize the unitary, and not the (dilated) Hamiltonian, time is a parameter and so we can simulate arbitrary time scales. Specifically, we can simulate the long time dynamics expected at the EP, due to the vanishing eigenvalue gap.
Interestingly, since the Hamiltonian is transposesymmetric, i.e. H † N (γ) = H * N (γ), the global unitary U 2N (γ, t) describes two coupled PT -symmetric systems that are the time reverse of each other. Indeed, G † N (γ, t) =G * N (γ, t) = TG N (γ, t)T −1 , meaning that a system initialised in the last N modes of the interferometer transforms with a backwards-in-time (reverse) evolution. Whenever γ = 0, the two systems are coupled by the defect matrix D N and excitations can tunnel between these two systems. In addition to simulating individual PT systems, the statistics output from the photonic chip track the unitary evolution of the closed composite system.
We also note that the above dilation process is not specific to PT -Hamiltonians. Any non-Unitary evolution operator, G N (t), can be embedded in a larger unitary transformation, following the above method. This process has already been used to simulate Lindbladian dynamics [32] and could be extended to simulate anti-PT systems [33]. (a) Graphical representation of the simulated PT -symmetric non-Hermitian Hamiltonian HN (γ) that describes an open quantum system. Excitations, represented by yellow circles, can populate the N modes of a linear chain. Nearest neighbour modes are coupled by a potential J that is set to unity during the simulations reported in this work. The imaginary potential, ±iγ, has the effect of amplifying and deamplifying the probability amplitude of an excitation in the first and last mode, respectively. H † N (γ) is equivalent to HN (γ) with the first and last modes swapped. The time evolution of a system with Hamiltonian HN (γ) is described by a non-unitary operator GN (γ, t) (b) Simulation procedure based on unitary dilation. Given a non-Hermitian Hamiltonian H(γ) and an evolution time t, it is possible to create a unitary U (γ, t) that consists of G(γ, t) = G(γ, t)/||G(γ, t)|| and an auxiliary matrix D(γ, t). The simulation proceeds by implementing U (γ, t) and initialising a desired input state in modes of the device corresponding to the forward system input. The effect of the unitary device will transform the input and distribute it across the forward and reverse output modes. When the output state is detected in the forward output, the input state successfully transformed viaG(γ, t). (c) Photonic implementation of the quantum simulation. Single photon states are injected into the input ports of the first N = 3 waveguides of a configurable integrated photonic circuit. Forward and reverse output are detected at the first and last output ports of the interferometer. A triangular mesh of integrated beam-splitters and metallic phaseshifters allows us to reconfigure the integrated interferometer to implement the desired U (γ, t).
by thirty metallic thermo-optic phase shifters, and is capable of implementing any unitary linear optical transformation over its six waveguides (see Ref. [27] for additional device details). By means of the algorithm demonstrated by Reck et al. for the decomposition of any discrete optical unitary operator [34], for each U 2N (γ, t i ) we retrieve the list of phases which, when applied in our interferometer, implement the desired unitary transformation. Calibration of the on-chip thermo-optic phase shifters allows us to map a desired phase shift to a voltage bias to be set. Input photons are generated with a bulk spontaneous parametric down-conversion source, and coupled in and out of the chip with packaged optical fibres, before detection with silicon avalanche photodiodes. This setup is capable of simulating coupled PTsymmetric Hamiltonians with N = 2 and N = 3 modes in each subsystem and multi-particle evolution.
Tuning the strength of the complex potential γ corre-sponds to changing the unitary matrix U 2N (γ, t i ), and therefore the set of phases dialled on our chip. In this way, all values of γ can be implemented in the same way, allowing us to investigate the system dynamics in four regimes: uncoupled Hermitian evolution, coupled PTsymmetric evolution, evolution at the EP, and evolution in the broken symmetry phase. As we increase γ we also increase the coupling between the systems, iD(γ, t), which correspondingly decreases the probability of a photon remaining in the desired subsystem.
To distinguish between the state of the simulated system and the photonic Fock state of our simulator, we will use the subscript p for the latter. That is, for N = 2, Fock states {|1 p 0 p , |0 p 1 p }, respectively, encode the states {|1 , |2 }, while for N = 3, Fock states {|1 p 0 p 0 p , |0 p 1 p 0 p , |0 p 0 p 1 p } encode the states {|1 , |2 , |3 }. The Fock states |0 p 0 p 0 p or |0 p 0 p , corresponding to an unoccupied system, will be represented Simulation of the forward system across different regimes of the PT -Hamiltonian. Solid lines (dots) correspond to the theory (experimental data). (a,b) A single excitation is initialised in |1 F and observed in the {|vac , |1 , |2 } F basis as shown with red, black and blue data, respectively. (a) In the unbroken symmetry phase (γ = 0.25), the plot shows the tunnelling of excitations between the forward and reverse systems resulting in a periodic oscillation in the vacuum probability; the inset shows the Hermitian regime (γ = 0) with a constant vacuum probability of 0, and sinusoidal oscillations between the two populated levels. (b) In the broken symmetry phase (γ = 1.1), the evolution is not periodic and the system tends rapidly to a steady state; the inset shows the evolution of the system at the exceptional point. (c,d) Evolution of the initial state The entropy of the system is periodic for γ < 1 and goes asymptotically to 0 for γ > 1. Time is reported in units normalised as a function of γ by dividing by τ = π/ |1 − γ 2 |. (d) Signalling violation as a function of γ and for evolution time t = π/ |1 − γ 2 | with a maximum at the exceptional point γ = 1. The outcomes of measurements on one qubit are affected by local operations on a PT qubit, entangled with the first.
with |vac . Furthermore, subscript F or R will differentiate between the forward and reverse systems, or equivalently, between the first or last N modes of our interferometer, respectively.

Simulation of two-mode PT symmetric systems
We first simulate 2-mode PT -symmetric systems. With N = 2, the eigenvalues of the Hamiltonian are = ± 1 − γ 2 . Since the global model in this instance comprises 4 modes in total, there is sufficient redundancy in the six-mode chip to implement a separate tomography stage. The first section of the chip implements the state evolution described by {U 4 (γ, t i )} i . The remaining part of the chip is used to perform projective measurements allowing the reconstruction of the density matrices for evolved single photon states. Figure 2 shows the results from the simulations of the unbroken and broken symmetry phases. Experimental points are presented next to theory lines. In Figs 2(a,b) we report the probability for an excitation initialised in the state |1 F ⊗ |vac R to be detected in any of the forward system modes but we enforce that the cumulative outcome statistic from the forward and reverse subsystems is normalised to one. To do so, we retained both the events where the propagating excitation was detected in the forward subsystem modes heralding its evolution according toG, and those where the excitation has left the forward modes to reappear in the reverse subsystem modes. If the excitation moves to the reverse system, the forward system is left in its vacuum state with no excitation, |vac F , and as such is reported in Figs 2(a,b). A formal derivation of the forward state description as Simulation of Zitterbewegung-like interference effects between coupled, time-reversed, two-mode PT -systems. result of a partial trace over the modes of the reverse system is available in the AppendixC. In the unbroken symmetry phase, shown in Fig. 2(a) with γ = 0.25, the time-evolution of the system is periodic and the finite probability of detecting |vac F indicates that the excitation can tunnel from the forward system to its time reversed twin. In contrast, the inset in the same plot shows the evolution when γ = 0. Since the system is closed and Hermitian, only sinusoidal oscillations of the probabilities of detecting |1 F and |2 F are observed. Figure 2(b) shows, instead, the evolution in the broken symmetry phase with γ = 1.1. The system rapidly tends to a steady state with no oscillations and with a non-zero probability of finding the excitation in either the forward or reverse systems. Such a qualitative distinction between the two regimes is caused by the characteristic phase transition at the EP. In the inset of Fig. 2(b), we report the evolution observed at the EP when γ = 1. Here the period of the system oscillation τ = π/ 1 − γ 2 tends to infinity. The probabilities evolve towards a never-reached turning point and the systems effectively tends to a steady state with coefficients that follows a polynomial power law. As opposite, in the broken symmetry phase, the system approaches its steady state with an exponential power law whose characteristic time, τ = π/ γ 2 − 1, decreases for increasing γ. Figures 2(c,d) show the evolution of a completely mixed input state, ρ 0 = 0.5 |1 F 1| F + 0.5 |2 F 2| F , which was realised by first entangling two photons using a photonic gate [35]. The total unitary implemented by the chip is obtained by multiplying together the transfer matrix that describes the entangling gate, and U 4 (γ, t), in such a way that one photon is injected into the 2 modes of the forward system of U 4 (γ, t) while the second photon of the entangled pair, is transmitted by two ancillary waveguides not affected by the U 4 (γ, t). To herald the successful generation of the entangled photon pair, the tomographic reconstruction of the state is performed by retaining only the coincidental events of a photon detected in the forward PT -symmetric system and a photon detected in the two ancillary waveguides.
In Fig. 2(c) we report the evolution of the Von Neumann entropy S[ρ(t)], of the initially mixed PT qubit ρ 0 , for different values of γ. As opposite to a closed Hermitian system whose entropy is a constant of motion, in the unbroken symmetry phase, the entropy is a non-monotonic function of time. By reporting the time parameter in units of τ = π/ |1 − γ 2 |, the entropy has a minimum when t = 0.5τ . Instead, when γ > γ c , the entropy decays to 0. The periodic behavior of the entropy in the PT -symmetric phase indicates the flow of information from the forward system to its time-reversed twin and back [36].
Figure 2(d) shows the results of a signalling violation test, as described in [37,38]. The test refers to the possiblity of affecting the state of a system B by performing only local operation on a system A. For a pair of qubits, the signalling violation, s.v., expressed in terms of the measurement probability P on the qubit B, reads: where 1(2) corresponds to the first (second) level of the qubit B, and I(S) refers to a local identity (swap) transformation applied to the qubit A. Our simulation proceeds by entangling the qubits A and B, then we apply either the identity or the swap transformation, and finally we simulate the evolution of qubit A with the PT -symmetric Hamiltonian. Meanwhile, the qubit B is encoded by the photon transmitted over the two ancillary waveguides. Figure 2(d) shows that the probabili-ties of the two alternative outputs for qubit B are conditioned by the a local unitary transformation applied to the PT qubit A prior to its non-unitary evolution. The no-signaling condition is violated for all γ = 0, and is maximal for γ = γ c . To maximise the signalling violation, the evolution time for the PT system is set to t = π/ |1 − γ 2 |.
Since in our framework excitations initialised in one system naturally tunnel to its time-reversed twin, we investigate the interference effects that arise from an initial excitation coherently superposed across the forward and the reverse systems. Interference between systems evolving in opposite directions in time is an effect pertinent to the foundations of quantum mechanics [39], and gives rise to peculiar quantum effects predicted in the dynamics of a relativistic electron [40,41], recently investigated with trapped ions and Bose-Einstein condensates [42,43].
We choose to measure an observable that manifests a qualitative change depending both on the coherence between the two systems and their deviation from hermiticity. We define a generalisation of the σ x operator for the forward (reverse) system as represented in the {|vac F(R) , |1 F(R) , |2 F(R) } basis. Such an observable is a constant of motion in the Hermitian regime and can be considered as a discrete-variable equivalent of a particle coordinate operator [44]. However, if we prepare an input state that is a coherent superposition of the forward and reverse subsystems, such as (|1 F + |1 R )/ √ 2, both S F and S R exhibit oscillatory behaviour.
The amplitude of these oscillations is reduced if the coherence between the states |1 F and |1 R decreases and it vanishes if the input is a statistical mixture of |1 F and |1 R , Fig. 3(a). We therefore associate this oscillatory behaviour with the coherent interference between forward and reverse modes of the PT -symmetric doublet, a situation that is analogous to the Zitterbewegung effect [40]. This predicts an oscillatory or "trembling" behaviour of the position expectation value for a relativistic electron, stemming from interference between positive and negative energy solutions of the Dirac equation. For this experiment, the mixed state spread across the forward and reverse system is simulated by statistically averaging the results from two orthogonal pure single particle input states. We finally note that in the broken symmetry phase, Fig. 3(b), S F(R) no longer exhibits oscillatory behaviour since the system loses its periodicity, but the final configuration of the system changes greatly depending on the initial coherence between the two subsystems.
Three-mode and many-particle PT -symmetric systems We have up to now discussed the non-unitary evolution of a single excitation over 2 modes, thereby simulating single-particle PT -symmetric systems. We now fix N = 3. In addition to gain and loss modes, the threemode PT -symmetric systems possess a neutral mode. With N = 3, the eigenvalues of H N (γ) become ∈ {0, 2 − γ 2 , − 2 − γ 2 }. We will refer to the AppendixF for an analysis of the evolution of a single excitation over 3 modes. Instead, here we show how the same interferometer set to realise the evolution described by U 2N (γ, t) can be used to study the interference of many identical photons undergoing non-unitary time evolution under a PT -symmetric Hamiltonian.
Using pairs of indistinguishable photons, we observe the evolution of the following two-particle input states, |Ω 1 has one photon in the gain and neutral modes each, while |Ω 2 has one photon in the loss and neutral modes each. Figures 4(a-c) show the probability for the system to collapse onto gain-neutral, gain-loss, and neutral-loss modes, i.e. states |1 p 1 p 0 p F , |1 p 0 p 1 p F and |0 p 1 p 1 p F respectively as a function of time. We choose four values of γ, in the Hermitian, PT -symmetric (γ = √ 2/2,γ = 3 √ 2/4), and PT -symmetry broken regimes (γ = 1.1 √ 2). A similar procedure can be adopted to simulate the non-Hermitian two-photon (Hong-Ou-Mandel) interference [26]. The probabilities reported in Fig 4(a-c) are normalised over all the possible antibunched patterns over the six modes of the inteferometer. We note that, due to the PT symmetry between these input states, we expect their evolution to proceed in opposite directions in time and to present an exchange of probabilities between the observation of the patterns |1 p 1 p 0 p and |0 p 1 p 1 p . The probability of detecting one photon in gain and neutral modes each at time t starting from state |Ω 1 is equal to the probability, at time −t, of detecting one photon in neutral and loss modes each when starting from state |Ω 2 . For these reasons, data collected for the input state |Ω 2 have been plotted in the reverse temporal direction to simplify the comparison with the input state |Ω 1 . By increasing the value of γ, the system undergoes the transition from the symmetric to the broken symmetry phase characterized by an increase of the period of the evolution, for γ < √ 2, and by a non-periodic evolution, for γ > √ 2. Lastly, we demonstrate the effect of non-hermiticity on photon-bunching by starting with a maximally antibunched three-photon input state, The plot shows the probabilities of detecting the output patterns |1p1p0p F ⊗ |0p0p0p R , |1p0p1p F ⊗ |0p0p0p R , and |0p1p1p F ⊗ |0p0p0p R , with dots corresponding to experimental data. Data collected for the state |Ω2 (in red) are plotted in the reverse temporal direction. Data for input state |Ω2 follow the same evolution as for input state |Ω1 , provided that the output patterns |1p1p0p F and |0p1p1p F are exchanged. (d-f) Three-particle evolution. When γ > 0, the output patterns related by a spatial inversion symmetry evolve in opposite directions in time. The curves for these related output patterns become identical when γ = 0. Error bars are obtained via bootstrapping methods [45].
To prepare the desired input state, we simultaneously collect two photon pairs, then use one of these photons to herald the generation of a three-photon anti-bunched state, which we inject into the chip. The data in Figs 4(df) show the evolution of the probability of different pho-ton bunching events in the forward system modes. Due to the underlying parity-time symmetry, the probability distribution satisfies Prob(t, |ψ ) = Prob(T (γ) − t, P |ψ ) where T (γ) = 2π/ |2 − γ 2 | is the fundamental period of a PT -symmetric trimer. Additionally, in the Hermi-tian limit, due to the underlying parity symmetry, the probability of two photons in mode 1 and one photon in mode 2 (i.e. pattern |2 p , 1 p , 0 p F ) is equal to that of its parity-symmetric pattern |0 p , 1 p , 2 p F . When γ > 0, the exchange symmetry between modes 1 and 3 is broken; the weight at smaller times t < T (γ)/2 shifts to the gain mode and while the weight at times t > T (γ)/2 shifts to the loss mode. For γ = 0.2, parity symmetric output patterns such as |2 p 0 p 1 p F and |1 p 0 p 2 p F show probability curves that are exchanged under time-inversion t → T (γ) − t.
It is computationally complex to simulate the evolution of quantum Hamiltonians [32] that produce particle statistics which are governed by intractable matrix functions. In the case of boson sampling [46], in which photons propagate in a circuit that is ideally described by a unitary matrix, classical intractability arises because the detection statistics are given by permanents of complex matrices. Yet photon loss reduces the computational complexity in experimental implementations of boson sampling [47]. Similarly, a photonic simulation of a traditional model for an open system, such as the Lindblad approach, could become classically tractable via the loss that is necessarily built into the model. However, in our simulation model for twin PT symmetric Hamiltonians, while each subsystem is open and lossy, the overall model is unitary and an ideal implementation is lossless; generally, statistics provide information about both the forward and backward subsystems, and the relation between them.

Discussion
The theory of PT symmetric Hamiltonians emerged two decades ago as a complex extension of quantum mechanics [48]. Over the past decade, it has galvanized research primarily in the classical (wave) domain [3,4]. With recent realizations of effective PT -symmetry and exceptional points in minimal quantum systems, the subject of non-Hermitian quantum systems is at the forefront again. Although most studies to date have been limited to single-particle systems, non-Hermitian quantum many-particle systems are an emerging and challenging frontier. This is true because most approximate methods -variational principle, tensor networks, perturbation theory -developed to reduce the exponential complexity of a many-body system, are designed for and work only with Hermitian Hamiltonians.
The task of simulating non-Hermitian, quantum manyparticle systems, including non-interacting bosons propagating with a non-Hermitian Hamiltonian, brings new challenges [49,50]. We have shown that a programmable, unitary simulator is suited to addressing these challenges. Combined with the increasing sophistication of integrated photonic devices, our simulation techniques will allow the investigation of PT -symmetric quantum systems beyond the computational capability of classical computers, paving the way for quantum technologies that can simulate and harness the potential advantages of non-Hermitian quantum systems.

Acknowledgments
We acknowledge support from the Engineering and Physical Sciences Research Council (EPSRC) Hub in Quantum Computing and Simulation (EP/T001062/1). Fellowship support from EPSRC is acknowledged by A.L. (EP/N003470/1). YJ was supported by NSF grant DMR-1054020 and thanks Institute for Advanced Studies, University of Bristol for their support.

Appendix A: Unitary dilation procedure
To dilate a non-unitary operator G N (t) of size N × N we must first ensure that its maximum singular value is less than 1, ||G(t)|| 2 1. We enforce this condition by scaling the operator G N (t) by a multiplicative factor α. Such normalisation does not alter quantity of the form where Π is a linear operator. Therefore α does not affect the resulting statistics of our simulations. We can either choose the scaling factor to be the maximum singular value over both time and spectrum of singular values, or choose to use a time varying scaling factor that is the maximum over the spectrum of singular values at time t. That is, G N (t) is normalised as (10) The second solution was preferred since it better adapts to the broken symmetry regime cases where the singular value of G N (t) increases exponentially with time.G N (t) can be decomposed by singular value decomposition as U cos − → θ V † , where U, V are N dimensional unitary matrices and cos − → θ = diag (cos θ 1 , . . . , cos θ N ). This allows forG N (t) to be embedded in the upper left corner of a 2N × 2N unitary matrix. Indeed, possible dilations of G N (t) can be written as (12) Furthermore, if C is a semi-positive definite diagonal matrix, by using that for a unitary matrix W and a ma- . Then, defining the defect operator as D N (t) = In the case when H N is Hermitian, G N (t) is unitary and the singular values are 1, so D N = 0, giving two independently evolving systems.

Appendix B: Effective Hamiltonian
As the evolution of the pair of coupled systems is unitary, it can be generated by a Hermitian Hamiltonian H eff over the 2N modes satisfying the equation From this, knowing U 2N it is possible to obtain H eff = i d dt U 2N U † 2N and numerically solve for H eff , which turns out to be time dependent. In Fig. 5 we show the results of the numerical calculation of the H eff elements for N = 2 and γ ∈ {0.25, 0.9, 1.1}. In the symmetric regime, the elements oscillate with a recurrence time τ = 2π/ 1 − γ 2 so we set an equivalent 'relaxation time' τ = 2π/ |1 − γ 2 | in the broken regime. The unplotted elements in figure are all (numerically) 0 including all diagonal elements H k,k for 1 ≤ k ≥ 4. The spikes appearing in H 1,3 and H 2,4 represent asymptotes.
Appendix C: Forward and reverse systems with vacuum contribution The dynamics of the forward (F) and reverse (R) systems are recovered by the reduced states of the full 2Nmode system: Here, ρ(0) is the input photonic state and U 2N is the unitary operator acting on the photonic Fock space according to the transfer matrix U 2N (γ, t), which transforms the input-mode creation operators a † i as To simulate a single N -dimensional PT system, we use 2N optical modes and we assume that the optical system is constrained to the one-photon subspace. In this space a basis is where |0 p is the single mode vacuum state. Doing this, the density matrix at time t can then be denoted as However we can treat the modes 1 to N as composing the forward system and modes from N + 1 to 2N as composing the reverse system leading to the definition of two reduced density matrices obtained performing the partial trace over one of the systems. For the reduced density matrix of the forward system ρ F we obtain: p kk |vac vac| Similar results are obtained when the partial trace is performed over the forward system, showing that the partial trace operation over the reverse (forward) system sums the diagonal elements of the bottom (top) N modes of ρ (t) and places this value in the vacuum state of the forward (reverse) system. That is, the Fock spaces of the systems are decomposed into particle number subspaces as . According to this, we can choose the basis states of the forward or reverse system as and an additional basis state in each system |vac = |0 p ⊗N . Furthermore, no coherence is present between the 1 photon and 0 photon subspaces.
Whenever there is interest in considering only the original PT system without introducing the reverse system and the interaction terms between them, we need to redefine the density matrix of the forward system without including the vacuum contribution as This approach has been used in the simulation of the no-signalling violation.

Photon generation
For all the experiments described in this paper, the single photons are generated via spontaneous parametric down-conversion in a Bismuth Borate non-linear crystal satisfying type I phase-matching conditions. The crystal is pumped with light pulses at 404 nm obtained from a second harmonic generation process induced in a Barium Borate crystal by a 808 nm Ti:Sapphire modelocked laser. The pairs of photons generated in the downconversion are spectrally filtered with a 3 nm bandwidth interferometric filter centred around 808 nm and are collected from two diametrically opposite positions of the generation cone. The pair of photons, coupled to single mode polarisation maintaining fibres, are used either as a two-photon Fock state, by injecting them both into our interferometer, or as heralded single photon by connecting one of the two fibres to a single photon detector and injecting the other photon into the interferometer. When interested in the two-photon Fock state we record all the events where two single photons are detected at the output of the chip. For the experiments in which we only need a single photon source, we record all the events when both the heralding photon and the photon at the output of the chip are detected. Three-photon experiments are performed collecting two separate pairs of photons from the emission cone and recording the coincidental detection of a heralding photon from one pair and of three photons passing through the optical modes of the interferometer.

Integrated interferometer
The programmable interferometer used to implement the multiple unitary transfer matrices of the experiments was produced by the NTTGroup in Japan and is described in [27]. The triangular arrangement of 15 integrated Mach-Zehnder interferometers with 30 thermally tunable phase-shifters allows us to prepare any unitary optical transfer matrix across 6 modes. Alternatively, the same device can be used as a cascade of 2×2 interferometers in order to sequentially implement state preparation, evolution or projection in physically different parts of the chip.

Probabilistic number resolving
Three photons simulations reported in this article require the ability to resolve multiple occupancy in the output modes to detect bunching photon events. To measure the probability of these events, we performed probabilistic number resolving detection of up to two photons in the same mode using auxiliary fibre beamsplitters (FBS) and detectors. By inserting a FBS at each of the first 3 output modes of the inteferometer and connecting both output modes of each FBS to detectors, there is a finite probability that photons bunched in the same output mode separate at the FBS and generate a signal at both detectors to which the FBS is connected.

Detection calibration
We experimentally estimate the probability of collapsing the quantum state of our simulator onto a particulat Fock state by multiplying the number of occurrence of each detected pattern by a detector calibration factor associated to that pattern, and then normalising to unity the sum over all relevent output patterns. The relative detector efficiency at the different output modes, for single photon events, is obtained by testing 500 random phases configurations of our interferometer and comparing the experimental statistic with the expected photon statistics due to the phase settings. The efficiency corrections of each mode that maximises the statistical fidelity between the experimental statistics and the theoretical distribution is chosen as numerical correction for the data collected during single photon experiments. For two-photon anti-bunched events, the occurence of each pattern is multiplied by the product of the correction factor of the individual output modes populated with photons.
In the pseudo-number-resolving configuration, the detection calibration factor is obtained by using the single photon detector efficiency ν i recorded at each output of the FBS. To take into account the probabilistic nature of the pseudo-number resolving scheme, the experimental occurrence of a detection pattern such as ((n 1,1 , n 1,2 ), (n 2,1 , n 2,2 ), (n 3,1 , n 3,2 )) is divided by Γ((n 1,1 , n 1,2 ), (n 2,1 , n 2,2 ), (n 3,1 , n 3,2 )) The number of events associated to an output state |N 1p N 2p N 3p with N kp = n k,1 + n k,2 is given as the average of the occurrences of equivalent detection patterns. For example, ((1, 1), (1, 0), (0, 0)) is equivalent to ((1, 1), (0, 1), (0, 0)) since both correspond to the output state |2 p 1 p 0 p . The final probabilities are normalised over all the output states characterised by three photons in the forward system modes but with less than three photons bunched in the same mode.

State preparation and reconstruction
The preparation of fully mixed single photon states is obtained by interfering two indistinguishable single photons. The Hong Ou Mandel visibility of their inteference is measured to be 0.97. During the entropy evolution simulation, the inteferometer is set to implement U Proj U 4 (γ, t) U Ent , where : U 4 (γ, t) transforms the bottom four modes of the interferometer, and U Proj projects the bottom two modes on the bases of the Pauli operators. Physically separate components of the optical circuit are involved in the implementation of the three cascaded unitaries. The evolution of excitations superposed across the forward and reverse systems is performed by compiling the two transformations required for preparing and evolving the states. The resulting unitary is implemented in the first part of the chip. The expectation values of S F and S R are measured using 3 further concatenated interferometers to project each of the systems, independently, on the eigenstates of the Pauli operators. When N = 3, to simulate the evolution of a single particle state prepared in a superposition across multiple modes of the forward or reverse system, we multiplied the transfer matrix U 6 (γ, t) with the state rotation necessary to transform a localised single photon into the desired superposition. The resulting unitary transformation is then dialed on the optical chip.
The experiments with two-photon Fock states are performed setting the transfer matrix of the chip to implement U 6 (γ, t) while the input state is obtained injecting either |0 p 1 p 1 p F ⊗|0 p 0 p 0 p R or |1 p 1 p 0 p F ⊗|0 p 0 p 0 p R into our integrated interferometer. The production of the desired input state is assumed everytime a coincidence detection happens at two of the detectors connected to the output modes of our device neglecting higher order pair generation occurring with a two order of magnitude lower rate.
Appendix E: Two-mode system: reverse subsystem In the main text we show the dynamics of an excitation in the forward subsystem for a two-mode PT -system. We included a combined vacuum probability which contained , and |ψ2(t = 0) = (i |2 F + |3 F )/ √ 2 (bottom row). Solid lines in the top rows represent theory curves for the time evolution of |ψ1(t) while solid lines in the bottom rows represent theory curves for the time evolution of |ψ2(t) . The data for |ψ3 show the same temporal behaviour as for |ψ1 . Conversely, the data for |ψ2 follow the opposite temporal evolution as |ψ1 , while the probabilities of observing |1 and |3 are exchanged.
the probability of the excitation tunnelling into the reverse subsystem. Using the same experimental data, in Fig. 6 we show the evolution of both forward and both reverse modes. This captures the full dynamics of a single excitation in the forward subsystem. We show the data for both the unbroken (γ = 0.25) and broken (γ = 1.1) symmetry regimes. We see a periodic oscillation in the probability of all four modes in the unbroken regime. In the broken regime all four mode probabilities evolve exponentially towards a steady state value.
Appendix F: Three-mode system: single particle dynamics Since, for N = 3, the third-order exceptional point occurs at γ c = √ 2, to analyse unbroken and broken symmetry regimes we perform simulations for γ = √ 2/2 and γ = 1.1 √ 2, respectively. Transfer matrices, concatenating state preparation and unitary evolution described by U 6 (γ, t i ), are implemented by using the entirety of the six mode chip. Figure 7 shows data reproducing the evolution of three single-particle initial states {|ψ 1 (t = 0) F , |ψ 2 (t = 0) F , |ψ 3 (t = 0) R }. These are related by symmetry transfor-mations as follows: is an example of a state coherently superposed across gain and neutral sites of the forward system; is the result of a PT transformation of |ψ 1 (0) F ; is prepared in the reverse system and, via time inversion, is analogous to the state T |ψ 1 (0) F . Due to the symmetries of our model, the mode occupation of these states during their evolution is closely related. Specifically, we expect: and, since H 3 (γ) = H T 3 (γ), we also have k|G 3 (t) † ψ 3 (0) = k|T G 3 (t)T ψ 3 (0) = k|T G 3 (t)ψ 1 (0) = k|G 3 (t)ψ 1 (0) * Figures 7(a,b) show that the (forward) evolution of state |ψ 1 with Hamiltonian H 3 is the same as the (reverse) evolution of the state |ψ 3 with Hamiltonian H * 3 . These results at time −t match the evolution of state |ψ 2 with Hamiltonian H 3 at time t, demonstrating the mirror symmetry about t = 0.