# Photonic quantum simulations of coupled $\mathcal{PT}$-symmetric
Hamiltonians

Nicola Maraviglia  Quantum Engineering Technology Labs, H. H. Wills Physics
Laboratory and Department of Electrical and Electronic Engineering, University
of Bristol, Bristol BS8 1FD, UK  Patrick Yard  Quantum Engineering Technology
Labs, H. H. Wills Physics Laboratory and Department of Electrical and
Electronic Engineering, University of Bristol, Bristol BS8 1FD, UK  Ross
Wakefield  Quantum Engineering Technology Labs, H. H. Wills Physics Laboratory
and Department of Electrical and Electronic Engineering, University of
Bristol, Bristol BS8 1FD, UK  Jacques Carolan  Center for Hybrid Quantum
Networks (Hy-Q), Niels Bohr Institute, University of Copenhagen, Blegdamsvej
17, DK-2100 Copenhagen, Denmark  Chris Sparrow  Quantum Engineering Technology
Labs, H. H. Wills Physics Laboratory and Department of Electrical and
Electronic Engineering, University of Bristol, Bristol BS8 1FD, UK  Department
of Physics, Imperial College London, SW7 2AZ, UK  Levon Chakhmakhchyan
Quantum Engineering Technology Labs, H. H. Wills Physics Laboratory and
Department of Electrical and Electronic Engineering, University of Bristol,
Bristol BS8 1FD, UK  Present address: Krisp, Hrachya Kochar 5/1, Yerevan 0033,
Armenia  Chris Harrold  Quantum Engineering Technology Labs, H. H. Wills
Physics Laboratory and Department of Electrical and Electronic Engineering,
University of Bristol, Bristol BS8 1FD, UK  Toshikazu Hashimoto  NTT Device
Technology Laboratories, NTT Corporation, Atsugi, Japan  Nobuyuki Matsuda  NTT
Basic Research Laboratories, NTT Corporation, Atsugi, Japan  Currently at:
Department of Communications Engineering, Graduate School of Engineering,
Tohoku University, Sendai, Japan  Andrew K. Harter  Department of Physics,
Indiana University Purdue University Indianapolis (IUPUI), Indianapolis,
Indiana 46202, USA  Yogesh N. Joglekar<EMAIL_ADDRESS>Department of
Physics, Indiana University Purdue University Indianapolis (IUPUI),
Indianapolis, Indiana 46202, USA  Anthony Laing<EMAIL_ADDRESS>Quantum Engineering Technology Labs, H. H. Wills Physics Laboratory and
Department of Electrical and Electronic Engineering, University of Bristol,
Bristol BS8 1FD, UK

###### Abstract

Parity-time ($\mathcal{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
$\mathcal{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 $\mathcal{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 ($\mathcal{P}$) and time-reversal
($\mathcal{T}$) transformations were extensively investigated, first
mathematically [1, 2] and then experimentally in classical wave systems [3,
4]. In the classical domain, $\mathcal{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, $\mathcal{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 $\mathcal{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 $\mathcal{PT}$ symmetry, is based on the
equivalence (after normalisation) of a $\mathcal{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
$\mathcal{PT}$-symmetric 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 post-selection 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
$\mathcal{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
$\mathcal{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 $\mathcal{PT}$-symmetric regime (defined by real
eigenvalues) and across the EP into the $\mathcal{PT}$-symmetry broken regime
(defined by complex conjugate 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.

##  .1 Model and simulation procedure

Our simulation approach begins by considering a family of
$\mathcal{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 $\pm i\gamma$ 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

$\displaystyle
H_{N}(\gamma)=-\sum_{k=1}^{N-1}(\left|k\right\rangle\\!\\!\left\langle
k+1\right|+\left|k+1\right\rangle\\!\\!\left\langle k\right|$ $\displaystyle)$
$\displaystyle+i\gamma(\left|1\right\rangle\\!\\!\left\langle 1\right|$
$\displaystyle-\left|N\right\rangle\\!\\!\left\langle N\right|)$ (1)

where $\left|k\right\rangle$ is the state associated to the $k$-th mode. Using
$\\{\left|k\right\rangle\\}_{k=1}^{N}$ as a basis, the parity operator we
choose is represented by the antidiagonal matrix
$\mathcal{P}_{mn}=\delta_{m,N+1-n}$ and the time-reversal operator
$\mathcal{T}$ corresponds to the complex conjugation operation.
$H_{N}(\gamma)$ presents an EP for a critical value $\gamma=\gamma_{\text{c}}$
that depends on $N$ [28, 29]. The spectrum of the Hamiltonian $H_{N}(\gamma)$
is purely real for $\gamma<\gamma_{\text{c}}$ and the corresponding
eigenvectors are simultaneous eigenvectors of the antilinear $\mathcal{PT}$
operator (unbroken symmetry phase). At the EP, $\gamma=\gamma_{\text{c}}$, two
or more states coalesce, while at $\gamma>\gamma_{\text{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
$\gamma_{\text{c}}$ is given by a simple closed form: $\gamma_{\text{c}}/J=1$
whenever $N$ is even and $\gamma_{\text{c}}/J=\sqrt{(N+1)/(N-1)}$ when $N$ is
odd. For a static, non-Hermitian Hamiltonian $H_{N}(\gamma)$, the non-unitary
time evolution operator $G_{N}(\gamma,t)$ is given by (assuming $\hbar=1$)

$G_{N}(\gamma,t)=\exp\left[-iH_{N}(\gamma)\,t\right].$ (2)

We note that these choices set the time-scale for our system, $\hbar/J$, to
unity.

Figure 1:  (a) Graphical representation of the simulated
$\mathcal{PT}$-symmetric non-Hermitian Hamiltonian $H_{N}(\gamma)$ 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, $\pm i\gamma$, has the effect
of amplifying and deamplifying the probability amplitude of an excitation in
the first and last mode, respectively. $H_{N}^{\dagger}(\gamma)$ is equivalent
to $H_{N}(\gamma)$ with the first and last modes swapped. The time evolution
of a system with Hamiltonian $H_{N}(\gamma)$ is described by a non-unitary
operator $G_{N}(\gamma,t)$ (b) Simulation procedure based on unitary dilation.
Given a non-Hermitian Hamiltonian $H(\gamma)$ and an evolution time $t$, it is
possible to create a unitary $U(\gamma,t)$ that consists of
$\tilde{G}(\gamma,t)=G(\gamma,t)/||G(\gamma,t)||$ and an auxiliary matrix
$D(\gamma,t)$. The simulation proceeds by implementing $U(\gamma,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 via $\tilde{G}(\gamma,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(\gamma,t)$.

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}(\gamma,t)$ by its largest singular value
$||G_{N}(\gamma,t)||_{2}$ obtaining

$\tilde{G}_{N}(\gamma,t):=G_{N}(\gamma,t)/||G_{N}(\gamma,t)||_{2},$ (3)

so that the operator norm of the scaled operator $\tilde{G}_{N}(\gamma,t)$ is
unity. Then, defining the defect operator as

$D_{N}(\tilde{G}_{N}):=\left[{\mathbb{1}}-\tilde{G}_{N}(\gamma,t)\tilde{G}_{N}^{\dagger}(\gamma,t)\right]^{1/2},$
(4)

we construct the following unitary transformation $U_{2N}(\gamma,t)$ defined
on a Hilbert space twice the size. The block matrix representation of
$U_{2N}(\gamma,t)$ is

$U_{2N}(\gamma,t)=\begin{bmatrix}\tilde{G}_{N}(\gamma,t)&iD_{N}(\tilde{G}_{N})\\\
iD_{N}(\tilde{G}_{N}^{\dagger})&\tilde{G}^{\dagger}_{N}(\gamma,t)\end{bmatrix}.$
(5)

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
$\mathcal{PT}$-symmetric system is simulated by injecting single photons into
the optical circuit and reprogramming it to implement the unitary transfer
matrix $U_{2N}(\gamma,t)$ for a sequence of time steps
$\\{U_{2N}(\gamma,t_{i})\\}_{i}$, see Fig. 1. To simulate the evolution
determined by $H_{N}(\gamma)$, 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 $\mathcal{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 transpose-symmetric, i.e.
$H_{N}^{\dagger}(\gamma)=H_{N}^{*}(\gamma)$, the global unitary
$U_{2N}(\gamma,t)$ describes two coupled $\mathcal{PT}$-symmetric systems that
are the time reverse of each other. Indeed,
$\tilde{G}^{\dagger}_{N}(\gamma,t)=\tilde{G}^{*}_{N}(\gamma,t)=\mathcal{T}\tilde{G}_{N}(\gamma,t)\mathcal{T}^{-1}$,
meaning that a system initialised in the last $N$ modes of the interferometer
transforms with a backwards-in-time (reverse) evolution. Whenever $\gamma\neq
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
$\mathcal{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
$\mathcal{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-$\mathcal{PT}$ systems [33].

Figure 2:  Simulation of the forward system across different regimes of the
$\mathcal{PT}$-Hamiltonian. Solid lines (dots) correspond to the theory
(experimental data). (a,b) A single excitation is initialised in
$\left|1\right\rangle_{\text{F}}$ and observed in the
$\left\\{\left|\text{vac}\right\rangle,\left|1\right\rangle,\left|2\right\rangle\right\\}_{\text{F}}$
basis as shown with red, black and blue data, respectively. (a) In the
unbroken symmetry phase ($\gamma=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
($\gamma=0$) with a constant vacuum probability of $0$, and sinusoidal
oscillations between the two populated levels. (b) In the broken symmetry
phase ($\gamma=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
$\rho_{0}=0.5\left|1\right\rangle_{\text{F}}\\!\left\langle
1\right|_{\text{F}}+0.5\left|2\right\rangle_{\text{F}}\\!\left\langle
2\right|_{\text{F}}$. (c) The entropy of the system is periodic for $\gamma<1$
and goes asymptotically to 0 for $\gamma>1$. Time is reported in units
normalised as a function of $\gamma$ by dividing by
$\tau=\pi/\sqrt{|1-\gamma^{2}|}$. (d) Signalling violation as a function of
$\gamma$ and for evolution time $t=\pi/\sqrt{|1-\gamma^{2}|}$ with a maximum
at the exceptional point $\gamma=1$. The outcomes of measurements on one qubit
are affected by local operations on a $\mathcal{PT}$ qubit, entangled with the
first.

##  .2 Experimental procedure

Our photonic chip comprises a triangular mesh of 15 interconnected Mach-
Zehnder interferometers, controlled 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}(\gamma,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
photo-diodes. This setup is capable of simulating coupled
$\mathcal{PT}$-symmetric Hamiltonians with $N=2$ and $N=3$ modes in each
subsystem and multi-particle evolution.

Tuning the strength of the complex potential $\gamma$ corresponds to changing
the unitary matrix $U_{2N}(\gamma,t_{i})$, and therefore the set of phases
dialled on our chip. In this way, all values of $\gamma$ can be implemented in
the same way, allowing us to investigate the system dynamics in four regimes:
uncoupled Hermitian evolution, coupled $\mathcal{PT}$-symmetric evolution,
evolution at the EP, and evolution in the broken symmetry phase. As we
increase $\gamma$ we also increase the coupling between the systems,
$iD(\gamma,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
$\\{\left|1_{p}0_{p}\right\rangle,\left|0_{p}1_{p}\right\rangle\\}$,
respectively, encode the states
$\\{\left|1\right\rangle,\left|2\right\rangle\\}$, while for $N=3$, Fock
states
$\\{\left|1_{p}0_{p}0_{p}\right\rangle,\left|0_{p}1_{p}0_{p}\right\rangle,\left|0_{p}0_{p}1_{p}\right\rangle\\}$
encode the states
$\\{\left|1\right\rangle,\left|2\right\rangle,\left|3\right\rangle\\}$. The
Fock states $\left|0_{p}0_{p}0_{p}\right\rangle$ or
$\left|0_{p}0_{p}\right\rangle$, corresponding to an unoccupied system, will
be represented with $\left|\text{vac}\right\rangle$. 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.

##  .3 Simulation of two-mode $\mathcal{PT}$ symmetric systems

Figure 3:  Simulation of Zitterbewegung-like interference effects between
coupled, time-reversed, two-mode $\mathcal{PT}$-systems. Expectation value of
$S_{\text{F(R)}}=\left|1\right\rangle_{\text{F(R)}}\left\langle
2\right|_{\text{F(R)}}+\left|2\right\rangle_{\text{F(R)}}\left\langle
1\right|_{\text{F(R)}}$ for the system initialised as
${\rho_{\text{in}}=1/2(\left|1\right\rangle_{\text{F}}\left\langle
1\right|_{\text{F}}+\left|1\right\rangle_{\text{R}}\left\langle
1\right|_{\text{R}})+a(\left|1\right\rangle_{\text{F}}\left\langle
1\right|_{\text{R}}+\left|1\right\rangle_{\text{R}}\left\langle
1\right|_{\text{F}})}$ are shown. The parameter $a$ characterises the degree
of coherence between the time-reversed systems. (a) In the unbroken symmetry
phase ($\gamma=0.25$), the amplitude of the oscillations of
$\langle{S_{\text{F(R)}}}\rangle$ are proportional to the initial level of
coherence of the input state. (b) In the broken symmetry phase ($\gamma=1.1$),
the initial level of coherence of the input state dictates the saturation
level of $\langle{S_{\text{F(R)}}}\rangle$.

We first simulate 2-mode $\mathcal{PT}$-symmetric systems. With $N=2$, the
eigenvalues of the Hamiltonian are $\epsilon=\pm\sqrt{1-\gamma^{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}(\gamma,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
$\left|1\right\rangle_{\text{F}}\otimes\left|\text{vac}\right\rangle_{\text{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 to $\tilde{G}$, 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, $\left|\text{vac}\right\rangle_{\text{F}}$, and as such is
reported in Figs 2(a,b). A formal derivation of the forward state description
as 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
$\gamma=0.25$, the time-evolution of the system is periodic and the finite
probability of detecting $\left|\text{vac}\right\rangle_{\text{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
$\gamma=0$. Since the system is closed and Hermitian, only sinusoidal
oscillations of the probabilities of detecting
$\left|1\right\rangle_{\text{F}}$ and $\left|2\right\rangle_{\text{F}}$ are
observed. Figure 2(b) shows, instead, the evolution in the broken symmetry
phase with $\gamma=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
$\gamma=1$. Here the period of the system oscillation
$\tau=\pi/\sqrt{1-\gamma^{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,
$\tau^{\prime}=\pi/\sqrt{\gamma^{2}-1}$, decreases for increasing $\gamma$.

Figures 2(c,d) show the evolution of a completely mixed input state,
$\rho_{0}=0.5\left|1\right\rangle_{\text{F}}\\!\left\langle
1\right|_{\text{F}}+0.5\left|2\right\rangle_{\text{F}}\\!\left\langle
2\right|_{\text{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}(\gamma,t)$, in such a way that one photon is injected into the 2
modes of the forward system of $U_{4}(\gamma,t)$ while the second photon of
the entangled pair, is transmitted by two ancillary waveguides not affected by
the $U_{4}(\gamma,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
$\mathcal{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[\rho(t)]$,
of the initially mixed $\mathcal{PT}$ qubit $\rho_{0}$, for different values
of $\gamma$. 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
$\tau=\pi/\sqrt{|1-\gamma^{2}|}$, the entropy has a minimum when $t=0.5\tau$.
Instead, when $\gamma>\gamma_{\text{c}}$, the entropy decays to 0. The
periodic behavior of the entropy in the $\mathcal{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:

$s.v.=\Big{[}\text{P}(1|S)-\text{P}(2|S)\Big{]}-\Big{[}\text{P}(1|I)-\text{P}(2|I)\Big{]},$
(6)

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 $\mathcal{PT}$-symmetric Hamiltonian.
Meanwhile, the qubit $B$ is encoded by the photon transmitted over the two
ancillary waveguides. Figure 2(d) shows that the probabilities of the two
alternative outputs for qubit $B$ are conditioned by the a local unitary
transformation applied to the $\mathcal{PT}$ qubit $A$ prior to its non-
unitary evolution. The no-signaling condition is violated for all $\gamma\neq
0$, and is maximal for $\gamma=\gamma_{\text{c}}$. To maximise the signalling
violation, the evolution time for the $\mathcal{PT}$ system is set to
$t=\pi/\sqrt{|1-\gamma^{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 $\sigma_{x}$ operator for
the forward (reverse) system

$S_{\text{F(R)}}=\begin{pmatrix}0&0&0\\\ 0&0&1\\\ 0&1&0\end{pmatrix},$ (7)

as represented in the
$\\{\left|\text{vac}\right\rangle_{\text{F(R)}},\left|1\right\rangle_{\text{F(R)}},\left|2\right\rangle_{\text{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
$(\left|1\right\rangle_{\text{F}}+\left|1\right\rangle_{\text{R}})/\sqrt{2}$,
both $\langle{S_{\text{F}}}\rangle$ and $\langle{S_{\text{R}}}\rangle$ exhibit
oscillatory behaviour.

The amplitude of these oscillations is reduced if the coherence between the
states $\left|1\right\rangle_{\text{F}}$ and $\left|1\right\rangle_{\text{R}}$
decreases and it vanishes if the input is a statistical mixture of
$\left|1\right\rangle_{\text{F}}$ and $\left|1\right\rangle_{\text{R}}$, Fig.
3(a). We therefore associate this oscillatory behaviour with the coherent
interference between forward and reverse modes of the $\mathcal{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),
$\langle{S_{\text{F}\left(\text{R}\right)}}\rangle$ 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.

Figure 4:  (a-c) Two-photon correlation measurements. Data representing the
evolution of the input states
$\left|\Omega_{1}\right\rangle=\left|1_{p}1_{p}0_{p}\right\rangle_{\text{F}}\otimes\left|0_{p}0_{p}0_{p}\right\rangle_{\text{R}}$
and
$\left|\Omega_{2}\right\rangle=\left|0_{p}1_{p}1_{p}\right\rangle_{\text{F}}\otimes\left|0_{p}0_{p}0_{p}\right\rangle_{\text{R}}$.
The plot shows the probabilities of detecting the output patterns
$\left|1_{p}1_{p}0_{p}\right\rangle_{\text{F}}\otimes\left|0_{p}0_{p}0_{p}\right\rangle_{\text{R}}$,
$\left|1_{p}0_{p}1_{p}\right\rangle_{\text{F}}\otimes\left|0_{p}0_{p}0_{p}\right\rangle_{\text{R}}$,
and
$\left|0_{p}1_{p}1_{p}\right\rangle_{\text{F}}\otimes\left|0_{p}0_{p}0_{p}\right\rangle_{\text{R}}$,
with dots corresponding to experimental data. Data collected for the state
$\left|\Omega_{2}\right\rangle$ (in red) are plotted in the reverse temporal
direction. Data for input state $\left|\Omega_{2}\right\rangle$ follow the
same evolution as for input state $\left|\Omega_{1}\right\rangle$, provided
that the output patterns $\left|1_{p}1_{p}0_{p}\right\rangle_{\text{F}}$ and
$\left|0_{p}1_{p}1_{p}\right\rangle_{\text{F}}$ are exchanged. (d-f) Three-
particle evolution. When $\gamma>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 $\gamma=0$. Error bars are
obtained via bootstrapping methods [45].

##  .4 Three-mode and many-particle $\mathcal{PT}$-symmetric systems

We have up to now discussed the non-unitary evolution of a single excitation
over 2 modes, thereby simulating single-particle $\mathcal{PT}$-symmetric
systems. We now fix $N=3$. In addition to gain and loss modes, the three-mode
$\mathcal{PT}$-symmetric systems possess a neutral mode. With $N=3$, the
eigenvalues of $H_{N}(\gamma)$ become
$\epsilon\in\\{0,\sqrt{2-\gamma^{2}},-\sqrt{2-\gamma^{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}(\gamma,t)$ can be used to study the
interference of many identical photons undergoing non-unitary time evolution
under a $\mathcal{PT}$-symmetric Hamiltonian.

Using pairs of indistinguishable photons, we observe the evolution of the
following two-particle input states,

$\begin{split}\left|\Omega_{1}\right\rangle=\left|1_{p}1_{p}0_{p}\right\rangle_{\text{F}}\otimes\left|0_{p}0_{p}0_{p}\right\rangle_{\text{R}},\\\
\left|\Omega_{2}\right\rangle=\mathcal{PT}\left|\Omega_{1}\right\rangle=\left|0_{p}1_{p}1_{p}\right\rangle_{\text{F}}\otimes\left|0_{p}0_{p}0_{p}\right\rangle_{\text{R}}.\end{split}$
(8)

$\left|\Omega_{1}\right\rangle$ has one photon in the gain and neutral modes
each, while $\left|\Omega_{2}\right\rangle$ 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
$\left|1_{p}1_{p}0_{p}\right\rangle_{\text{F}}$,
$\left|1_{p}0_{p}1_{p}\right\rangle_{\text{F}}$ and
$\left|0_{p}1_{p}1_{p}\right\rangle_{\text{F}}$ respectively as a function of
time. We choose four values of $\gamma$, in the Hermitian,
$\mathcal{PT}$-symmetric ($\gamma=\sqrt{2}/2$,$\gamma=3\sqrt{2}/4$), and
$\mathcal{PT}$-symmetry broken regimes ($\gamma=1.1\sqrt{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 $\mathcal{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 $\left|1_{p}1_{p}0_{p}\right\rangle$ and
$\left|0_{p}1_{p}1_{p}\right\rangle$. The probability of detecting one photon
in gain and neutral modes each at time $t$ starting from state
$\left|\Omega_{1}\right\rangle$ is equal to the probability, at time $-t$, of
detecting one photon in neutral and loss modes each when starting from state
$\left|\Omega_{2}\right\rangle$. For these reasons, data collected for the
input state $\left|\Omega_{2}\right\rangle$ have been plotted in the reverse
temporal direction to simplify the comparison with the input state
$\left|\Omega_{1}\right\rangle$. By increasing the value of $\gamma$, the
system undergoes the transition from the symmetric to the broken symmetry
phase characterized by an increase of the period of the evolution, for
$\gamma<\sqrt{2}$, and by a non-periodic evolution, for $\gamma>\sqrt{2}$.

Lastly, we demonstrate the effect of non-hermiticity on photon-bunching by
starting with a maximally anti-bunched three-photon input state,

$\left|\chi_{1}\right\rangle=\left|1_{p}1_{p}1_{p}\right\rangle_{\text{F}}\otimes\left|0_{p}0_{p}0_{p}\right\rangle_{\text{R}}.$
(9)

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(d-f) show the evolution of the probability of different photon bunching
events in the forward system modes. Due to the underlying parity-time
symmetry, the probability distribution satisfies
$\mathrm{Prob}(t,\left|\psi\right\rangle)=\mathrm{Prob}(T(\gamma)-t,\mathcal{P}\left|\psi\right\rangle)$
where $T(\gamma)=2\pi/\sqrt{|2-\gamma^{2}|}$ is the fundamental period of a
$\mathcal{PT}$-symmetric trimer. Additionally, in the Hermitian limit, due to
the underlying parity symmetry, the probability of two photons in mode 1 and
one photon in mode 2 (i.e. pattern
$\left|2_{p},1_{p},0_{p}\right\rangle_{\text{F}}$) is equal to that of its
parity-symmetric pattern $\left|0_{p},1_{p},2_{p}\right\rangle_{\text{F}}$.
When $\gamma>0$, the exchange symmetry between modes 1 and 3 is broken; the
weight at smaller times $t<T(\gamma)/2$ shifts to the gain mode and while the
weight at times $t>T(\gamma)/2$ shifts to the loss mode. For $\gamma=0.2$,
parity symmetric output patterns such as
$\left|2_{p}0_{p}1_{p}\right\rangle_{\text{F}}$ and
$\left|1_{p}0_{p}2_{p}\right\rangle_{\text{F}}$ show probability curves that
are exchanged under time-inversion $t\rightarrow T(\gamma)-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.

##  .5 Discussion

The theory of $\mathcal{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 $\mathcal{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 many-particle 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 $\mathcal{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.

##  .6 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.

##  .7 Appendix A: Unitary dilation procedure

To dilate a non-unitary operator $G_{N}\left(t\right)$ of size $N\times N$ we
must first ensure that its maximum singular value is less than 1,
$||G(t)||_{2}\leqslant 1$. We enforce this condition by scaling the operator
$G_{N}\left(t\right)$ by a multiplicative factor $\alpha$. Such normalisation
does not alter quantity of the form

$\frac{\text{Tr}\left[\Pi\alpha G(\gamma,t)\rho
G^{\dagger}(\gamma,t)\alpha^{*}\right]}{\text{Tr}\left[\alpha G(\gamma,t)\rho
G^{\dagger}(\gamma,t)\alpha^{*}\right]}=\frac{|\alpha|^{2}}{|\alpha|^{2}}\frac{\text{Tr}\left[\Pi
G(\gamma,t)\rho G^{\dagger}(\gamma,t)\right]}{\text{Tr}\left[G(\gamma,t)\rho
G^{\dagger}(\gamma,t)\right]}$

where $\Pi$ is a linear operator. Therefore $\alpha$ 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}\left(t\right)$ is normalised as

$\tilde{G}_{N}\left(t\right)=\frac{G_{N}\left(t\right)}{\max_{t}\left[||G_{N}\left(t\right)||_{2}\right]}\quad\text{or}\quad\tilde{G}_{N}\left(t\right)=\frac{G_{N}\left(t\right)}{\lVert
G_{N}\left(t\right)\rVert_{2}}.$ (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. $\tilde{G}_{N}\left(t\right)$ can be decomposed by
singular value decomposition as $U\cos{\overrightarrow{\theta}}V^{\dagger}$,
where $U,V$ are $N$ dimensional unitary matrices and
$\cos{\overrightarrow{\theta}}=\text{diag}\left(\cos{\theta_{1}},\ldots,\cos{\theta_{N}}\right)$.
This allows for $\tilde{G}_{N}\left(t\right)$ to be embedded in the upper left
corner of a $2N\times 2N$ unitary matrix. Indeed, possible dilations of
$\tilde{G}_{N}\left(t\right)$ can be written as

$U_{2N}\left(\gamma,t\right)=\begin{pmatrix}U&&0\\\
0&&A\end{pmatrix}\begin{pmatrix}\cos{\overrightarrow{\theta}}&&i\sin{\overrightarrow{\theta}}\\\
i\sin{\overrightarrow{\theta}}&&\cos{\overrightarrow{\theta}}\end{pmatrix}\begin{pmatrix}V^{\dagger}&&0\\\
0&&B\end{pmatrix}$ (11)

for any $N$ dimensional unitary matrix $A$ and $B$ and for
$\sin{\overrightarrow{\theta}}=\text{diag}\left(\sqrt{1-\cos{\theta_{1}}^{2}},\ldots,\sqrt{1-\cos{\theta_{N}}^{2}}\right)$.
Setting $A=V$ and $B=U^{\dagger}$ gives

$U_{2N}\left(\gamma,t\right)=\begin{pmatrix}U\cos{\overrightarrow{\theta}}V^{\dagger}&&iU\sqrt{I-\cos^{2}{\overrightarrow{\theta}}}U^{\dagger}\\\
iV\sqrt{I-\cos^{2}{\overrightarrow{\theta}}}V^{\dagger}&&V\cos{\overrightarrow{\theta}}U^{\dagger}\end{pmatrix}$
(12)

Furthermore, if $C$ is a semi-positive definite diagonal matrix, by using that
for a unitary matrix $W$ and a matrix
$E=WCW^{\dagger},\,\sqrt{E}=W\sqrt{C}W^{\dagger}$, we obtain

$U_{2N}\left(\gamma,t\right)=\begin{pmatrix}\tilde{G}_{N}\left(t\right)&&i\sqrt{I-\tilde{G}_{N}\left(t\right)\tilde{G}_{N}^{\dagger}\left(t\right)}\\\
i\sqrt{I-\tilde{G}_{N}^{\dagger}\left(t\right)\tilde{G}_{N}\left(t\right)}&&\tilde{G}_{N}^{\dagger}\left(t\right)\end{pmatrix}.$
(13)

Since in our model $G_{N}\left(t\right)=G_{N}^{\mathrm{T}}\left(t\right)$, we
have $V=U^{*}$ and
$\left(\sqrt{I-\tilde{G}_{N}^{\dagger}\left(t\right)\tilde{G}_{N}\left(t\right)}\right)^{*}=\sqrt{I-\tilde{G}_{N}\left(t\right)\tilde{G}_{N}^{\dagger}\left(t\right)}$.
Then, defining the defect operator as
$D_{N}\left(t\right)=\sqrt{I-\tilde{G}_{N}\left(t\right)\tilde{G}_{N}^{\dagger}\left(t\right)}$,
we can rewrite $U_{2N}$ as

$U_{2N}\left(\gamma,t\right)=\begin{pmatrix}\tilde{G}_{N}\left(t\right)&&iD_{N}\left(t\right)\\\
iD_{N}^{*}\left(t\right)&&\tilde{G}_{N}^{\dagger}\left(t\right)\end{pmatrix}$
(14)

In the case when $H_{N}$ is Hermitian, $G_{N}\left(t\right)$ is unitary and
the singular values are $1$, so $D_{N}=0$, giving two independently evolving
systems.

Figure 5: non-zero elements of the effective Hamiltonian $H_{\text{eff}}$ of
the dilated unitary for the $2-$ mode $\mathcal{PT}$-symmetric Hamiltonian
$H_{2}$

##  .8 Appendix B: Effective Hamiltonian

As the evolution of the pair of coupled systems is unitary, it can be
generated by a Hermitian Hamiltonian $H_{\text{eff}}$ over the $2N$ modes
satisfying the equation

$H_{\text{eff}}U_{2N}=i\frac{\text{d}}{\text{d}t}U_{2N}.$ (15)

From this, knowing $U_{2N}$ it is possible to obtain
$H_{\text{eff}}=i\left(\frac{d}{dt}U_{2N}\right)U_{2N}^{\dagger}$ and
numerically solve for $H_{\text{eff}}$, which turns out to be time dependent.
In Fig. 5 we show the results of the numerical calculation of the
$H_{\text{eff}}$ elements for $N=2$ and $\gamma\in\\{0.25,0.9,1.1\\}$. In the
symmetric regime, the elements oscillate with a recurrence time
$\tau=2\pi/\sqrt{1-\gamma^{2}}$ so we set an equivalent ‘relaxation time’
$\tau=2\pi/\sqrt{|1-\gamma^{2}|}$ in the broken regime. The unplotted elements
in figure are all (numerically) $0$ including all diagonal elements $H_{k,k}$
for $1\leq k\geq 4$. The spikes appearing in $H_{1,3}$ and $H_{2,4}$ represent
asymptotes.

##  .9 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 $2N$-mode system:

$\begin{split}\rho_{\text{F}}(t)&=\text{Tr}_{\text{R}}\,\rho(t):=\text{Tr}_{\text{R}}\big{[}\mathcal{U}_{2N}(\gamma,t)\rho(0)\mathcal{U}_{2N}^{{\dagger}}(\gamma,t)\big{]},\\\
\rho_{\text{R}}(t)&=\text{Tr}_{\text{F}}\,\rho(t):=\text{Tr}_{\text{F}}\big{[}\mathcal{U}_{2N}(\gamma,t)\rho(0)\mathcal{U}_{2N}^{{\dagger}}(\gamma,t)\big{]}.\end{split}$
(16)

Here, $\rho(0)$ is the input photonic state and $\mathcal{U}_{2N}$ is the
unitary operator acting on the photonic Fock space according to the transfer
matrix $U_{2N}(\gamma,t)$, which transforms the input-mode creation operators
$a_{i}^{\dagger}$ as

$\mathcal{U}_{2N}(\gamma,t)a_{i}^{{\dagger}}\mathcal{U}^{{\dagger}}_{2N}(\gamma,t)=\sum_{j}[U_{2N}(\gamma,t)]_{ji}a^{{\dagger}}_{j}.$
(17)

To simulate a single $N$-dimensional $\mathcal{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

$\left\\{\left|k\right>=\left|0_{p}\right>^{\otimes\left(k-1\right)}\left|1_{p}\right>\left|0_{p}\right>^{\otimes\left(2N-k\right)}\right\\}_{k=1}^{2N}$
(18)

where $\left|0_{p}\right>$ is the single mode vacuum state. Doing this, the
density matrix at time $t$ can then be denoted as

$\rho\left(t\right)=U_{2N}\left(\gamma,t\right)\rho\left(0\right)U_{2N}^{\dagger}\left(\gamma,t\right).$
(19)

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 $\rho_{F}$ we obtain:

$\displaystyle\rho_{F}=\mathrm{tr}_{R}\rho$
$\displaystyle=\mathrm{tr}_{R}\sum_{k,l=1}^{2N}p_{kl}\left|k\rangle\langle
l\right|$
$\displaystyle=\mathrm{tr}_{R}\sum_{k,l=1}^{2N}p_{kl}\left|0_{p}\right>^{\otimes\left(k-1\right)}\left|1_{p}\right>\left|0_{p}\right>^{\otimes\left(2N-k\right)}\left<0_{p}\right|^{\otimes\left(l-1\right)}\left<1_{p}\right|\left<0_{p}\right|^{\otimes\left(2N-l\right)}$
$\displaystyle=\mathrm{tr}_{R}\biggl{(}\sum_{k,l=1}^{N}p_{kl}\left|0_{p}\right>^{\otimes\left(k-1\right)}\left|1_{p}\right>\left|0_{p}\right>^{\otimes\left(2N-k\right)}\left<0_{p}\right|^{\otimes\left(l-1\right)}\left<1_{p}\right|\left<0_{p}\right|^{\otimes\left(2N-l\right)}$
$\displaystyle\qquad\qquad+\sum_{k,l=N+1}^{2N}p_{kl}\left|0_{p}\right>^{\otimes\left(k-1\right)}\left|1_{p}\right>\left|0_{p}\right>^{\otimes\left(2N-k\right)}\left<0_{p}\right|^{\otimes\left(l-1\right)}\left<1_{p}\right|\left<0_{p}\right|^{\otimes\left(2N-l\right)}$
$\displaystyle\qquad\qquad+\sum_{k=1}^{N}\sum_{l=N+1}^{2N}p_{kl}\left|0_{p}\right>^{\otimes\left(k-1\right)}\left|1_{p}\right>\left|0_{p}\right>^{\otimes\left(2N-k\right)}\left<0_{p}\right|^{\otimes\left(l-1\right)}\left<1_{p}\right|\left<0_{p}\right|^{\otimes\left(2N-l\right)}$
$\displaystyle\qquad\qquad+\sum_{k=N+1}^{2N}\sum_{l=1}^{N}p_{kl}\left|0_{p}\right>^{\otimes\left(k-1\right)}\left|1_{p}\right>\left|0_{p}\right>^{\otimes\left(2N-k\right)}\left<0_{p}\right|^{\otimes\left(l-1\right)}\left<1_{p}\right|\left<0_{p}\right|^{\otimes\left(2N-l\right)}\biggr{)}$
$\displaystyle=\sum_{k,l=1}^{N}p_{kl}\left|0_{p}\right>^{\otimes\left(k-1\right)}\left|1_{p}\right>\left|0_{p}\right>^{\otimes\left(N-k\right)}\left<0_{p}\right|^{\otimes\left(l-1\right)}\left<1_{p}\right|\left<0_{p}\right|^{\otimes\left(N-l\right)}+\sum_{k,l=N+1}^{2N}p_{kl}\left|0_{p}\right>^{\otimes
N}\left<0_{p}\right|^{\otimes N}\delta_{k,l}$
$\displaystyle=\sum_{k,l=1}^{N}p_{kl}\left|k\rangle\langle
l\right|+\sum_{k=N+1}^{2N}p_{kk}\left|\text{vac}\rangle\langle\text{vac}\right|$

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
$\rho\left(t\right)$ 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
$\mathcal{F}^{\left(F/R\right)}=\mathcal{H}^{\left(F/R\right)}_{0}\oplus\mathcal{H}^{\left(F/R\right)}_{1}$.
According to this, we can choose the basis states of the forward or reverse
system as
$\left\\{\left|k\right>^{\left(F/R\right)}=\left|0_{p}\right>^{\otimes\left(k-1\right)}\left|1_{p}\right>\left|0_{p}\right>^{\otimes\left(N-k\right)}\right\\}_{k=1}^{N}$
and an additional basis state in each system
$\left|\text{vac}\right\rangle=\left|0_{p}\right>^{\otimes N}$. Furthermore,
no coherence is present between the 1 photon and 0 photon subspaces.

Whenever there is interest in considering only the original $\mathcal{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

$\tilde{\rho}_{F}=\frac{\sum_{k,l=1}^{N}p_{kl}\left|k\rangle\langle
l\right|}{\sum_{k=1}^{N}p_{kk}}$ (20)

This approach has been used in the simulation of the no-signalling violation.

##  .10 Appendix D: Experimental methods

###  .10.1 _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 mode-locked
laser. The pairs of photons generated in the down-conversion 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.

###  .10.2 _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\times 2$ interferometers in order to
sequentially implement state preparation, evolution or projection in
physically different parts of the chip.

###  .10.3 _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.

Figure 6: Plots showing the evolution of a single excitation in the forward
subsystem. Solid lines (points) show the theoretical curve (experimental
data). We show the probabilities for the excitation being detected in the
forward subsystem ($\left|1\right\rangle_{\mathrm{F}}$,
$\left|2\right\rangle_{\mathrm{F}}$) and the reverse subsystem
($\left|1\right\rangle_{\mathrm{R}}$, $\left|2\right\rangle_{\mathrm{R}}$).
Both panels follow the legend in (b). (a) Plots for a
$\mathcal{PT}$-Hamiltonian in the unbroken regime ($\gamma=0.25$). We see
periodic oscillations of the excitation in all four modes. (b) Plots for the
broken symmetry regime. All probabilities evolve exponentially to a steady
state ($\gamma=1.1$).

###  .10.4 _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 $\nu_{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

$\Gamma((n_{1,1},n_{1,2}),(n_{2,1},n_{2,2}),(n_{3,1},n_{3,2}))\\\
=\prod_{i=1}^{i=3}\nu_{i,1}^{n_{i,1}}\nu_{i,2}^{n_{i,2}}(n_{i,1}+n_{i,2})!$
(21)

The number of events associated to an output state
$\left|N_{1p}N_{2p}N_{3p}\right\rangle$ 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 $\left|2_{p}1_{p}0_{p}\right\rangle$. 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.

###  .10.5 _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_{\text{Proj}}\
U_{4}\left(\gamma,t\right)\ U_{\text{Ent}}$, where :

$U_{\text{Ent}}=\begin{pmatrix}0&0&\frac{\sqrt{2}}{2}&\frac{\sqrt{2}}{2}&0&0\\\
0&0&0&0&\frac{\sqrt{2}}{2}&\frac{\sqrt{2}}{2}\\\ 1&0&0&0&0&0\\\ 0&1&0&0&0&0\\\
0&0&\frac{\sqrt{2}}{2}&-\frac{\sqrt{2}}{2}&0&0\\\
0&0&0&0&\frac{\sqrt{2}}{2}&-\frac{\sqrt{2}}{2}\\\ \end{pmatrix},$ (22)

$U_{4}\left(\gamma,t\right)$ transforms the bottom four modes of the
interferometer, and $U_{\text{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.

Figure 7:  Results for the simulation of the 3-level $\mathcal{PT}$-symmetric
Hamiltonian. Plots show evolution of the probability, in the unbroken (a)
($\gamma=0.5\gamma_{\text{c}}$) and broken (b) ($\gamma=1.1\gamma_{\text{c}}$)
symmetry regimes, for three input states:
$\left|\psi_{1}(t=0)\right\rangle=(\left|1\right\rangle_{\text{F}}-i\left|2\right\rangle_{\text{F}})/\sqrt{2}$,
$\left|\psi_{3}(t=0)\right\rangle=(\left|1\right\rangle_{\text{R}}+i\left|2\right\rangle_{\text{R}})/\sqrt{2}$
(top row), and
$\left|\psi_{2}(t=0)\right\rangle=(i\left|2\right\rangle_{\text{F}}+\left|3\right\rangle_{\text{F}})/\sqrt{2}$
(bottom row). Solid lines in the top rows represent theory curves for the time
evolution of $\left|\psi_{1}(t)\right\rangle$ while solid lines in the bottom
rows represent theory curves for the time evolution of
$\left|\psi_{2}(t)\right\rangle$. The data for $\left|\psi_{3}\right\rangle$
show the same temporal behaviour as for $\left|\psi_{1}\right\rangle$.
Conversely, the data for $\left|\psi_{2}\right\rangle$ follow the opposite
temporal evolution as $\left|\psi_{1}\right\rangle$, while the probabilities
of observing $\left|1\right\rangle$ and $\left|3\right\rangle$ are exchanged.

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}\left(\gamma,t\right)$ 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}(\gamma,t)$ while the input state is
obtained injecting either
$\left|0_{p}1_{p}1_{p}\right\rangle_{\text{F}}\otimes\left|0_{p}0_{p}0_{p}\right\rangle_{\text{R}}$
or
$\left|1_{p}1_{p}0_{p}\right\rangle_{\text{F}}\otimes\left|0_{p}0_{p}0_{p}\right\rangle_{\text{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.

##  .11 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 $\mathcal{PT}$-system. We included a combined vacuum
probability which contained 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 ($\gamma=0.25$) and broken ($\gamma=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.

##  .12 Appendix F: Three-mode system: single particle dynamics

Since, for $N=3$, the third-order exceptional point occurs at
$\gamma_{\text{c}}=\sqrt{2}$, to analyse unbroken and broken symmetry regimes
we perform simulations for $\gamma=\sqrt{2}/2$ and $\gamma=1.1\,\sqrt{2}$,
respectively. Transfer matrices, concatenating state preparation and unitary
evolution described by $U_{6}(\gamma,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 {$\left|\psi_{1}(t=0)\right\rangle_{\text{F}}$,
$\left|\psi_{2}(t=0)\right\rangle_{\text{F}}$,
$\left|\psi_{3}(t=0)\right\rangle_{\text{R}}$}. These are related by symmetry
transformations as follows:

$\left|\psi_{1}(0)\right\rangle_{\text{F}}=(\left|1\right\rangle_{\text{F}}-i\left|2\right\rangle_{\text{F}})/\sqrt{2}$
(23)

is an example of a state coherently superposed across gain and neutral sites
of the forward system;

$\displaystyle\left|\psi_{2}(0)\right\rangle_{\text{F}}$
$\displaystyle=(i\left|2\right\rangle_{\text{F}}+\left|3\right\rangle_{\text{F}})/\sqrt{2}=\mathcal{PT}\left|\psi_{1}(0)\right\rangle_{\text{F}}$
(24)

is the result of a $\mathcal{PT}$ transformation of
$\left|\psi_{1}(0)\right\rangle_{\text{F}}$;

$\displaystyle\left|\psi_{3}(0)\right\rangle_{\text{R}}$
$\displaystyle=(\left|1\right\rangle_{\text{R}}+i\left|2\right\rangle_{\text{R}})/\sqrt{2}=\mathcal{T}\left|\psi_{1}(0)\right\rangle_{\text{R}}$
(25)

is prepared in the reverse system and, via time inversion, is analogous to the
state $\mathcal{T}\left|\psi_{1}(0)\right\rangle_{\text{F}}$.

Due to the symmetries of our model, the mode occupation of these states during
their evolution is closely related. Specifically, we expect:

$\displaystyle\left\langle k|\psi_{2}(t)\right\rangle$
$\displaystyle=\left\langle
k|\mathcal{PT}G_{3}(-t)\mathcal{PT}\psi_{2}(0)\right\rangle$
$\displaystyle=\left\langle 4-k|G_{3}(-t)\psi_{1}(0)\right\rangle^{*},$ (26)

and, since $H_{\text{3}}(\gamma)=H_{\text{3}}^{\text{T}}(\gamma)$, we also
have

$\displaystyle\left\langle k|G_{3}(t)^{\dagger}\psi_{3}(0)\right\rangle$
$\displaystyle=\left\langle
k|\mathcal{T}G_{3}(t)\mathcal{T}\psi_{3}(0)\right\rangle$
$\displaystyle=\left\langle k|\mathcal{T}G_{3}(t)\psi_{1}(0)\right\rangle$
$\displaystyle=\left\langle k|G_{3}(t)\psi_{1}(0)\right\rangle^{*}$ (27)

Figures 7(a,b) show that the (forward) evolution of state
$\left|\psi_{1}\right\rangle$ with Hamiltonian $H_{3}$ is the same as the
(reverse) evolution of the state $\left|\psi_{3}\right\rangle$ with
Hamiltonian $H_{3}^{*}$. These results at time $-t$ match the evolution of
state $\left|\psi_{2}\right\rangle$ with Hamiltonian $H_{3}$ at time $t$,
demonstrating the mirror symmetry about $t=0$.

## References

  * [1] C. M. Bender, “Making sense of non-hermitian hamiltonians,” Reports on Progress in Physics, vol. 70, pp. 947–1018, may 2007.
  * [2] A. Mostafazadeh, “Pseudo-hermitian representation of quantum mechanics,” International Journal of Geometric Methods in Modern Physics, vol. 07, pp. 1191–1306, Nov. 2010.
  * [3] R. El-Ganainy, K. G. Makris, M. Khajavikhan, Z. H. Musslimani, S. Rotter, and D. N. Christodoulides, “Non-hermitian physics and PT symmetry,” Nature Physics, vol. 14, pp. 11–19, Jan. 2018.
  * [4] M.-A. Miri and A. Alù, “Exceptional points in optics and photonics,” Science, vol. 363, p. eaar7709, Jan. 2019.
  * [5] L. Feng, Y.-L. Xu, W. S. Fegadolli, M.-H. Lu, J. E. B. Oliveira, V. R. Almeida, Y.-F. Chen, and A. Scherer, “Experimental demonstration of a unidirectional reflectionless parity-time metamaterial at optical frequencies,” Nature Materials, vol. 12, pp. 108–113, Nov. 2012.
  * [6] L. Feng, Z. J. Wong, R.-M. Ma, Y. Wang, and X. Zhang, “Single-mode laser by parity-time symmetry breaking,” Science, vol. 346, pp. 972–975, Oct. 2014\.
  * [7] H. Hodaei, M.-A. Miri, M. Heinrich, D. N. Christodoulides, and M. Khajavikhan, “Parity-time-symmetric microring lasers,” Science, vol. 346, pp. 975–978, Oct. 2014.
  * [8] G. Harari, M. A. Bandres, Y. Lumer, M. C. Rechtsman, Y. D. Chong, M. Khajavikhan, D. N. Christodoulides, and M. Segev, “Topological insulator laser: Theory,” Science, vol. 359, p. eaar4003, Feb. 2018.
  * [9] M. A. Bandres, S. Wittek, G. Harari, M. Parto, J. Ren, M. Segev, D. N. Christodoulides, and M. Khajavikhan, “Topological insulator laser: Experiments,” Science, vol. 359, p. eaar4005, Feb. 2018.
  * [10] H. Xu, D. Mason, L. Jiang, and J. G. E. Harris, “Topological energy transfer in an optomechanical system with exceptional points,” Nature, vol. 537, pp. 80–83, July 2016.
  * [11] J. Doppler, A. A. Mailybaev, J. Böhm, U. Kuhl, A. Girschik, F. Libisch, T. J. Milburn, P. Rabl, N. Moiseyev, and S. Rotter, “Dynamically encircling an exceptional point for asymmetric mode switching,” Nature, vol. 537, pp. 76–79, July 2016.
  * [12] Y.-H. Lai, Y.-K. Lu, M.-G. Suh, Z. Yuan, and K. Vahala, “Observation of the exceptional-point-enhanced sagnac effect,” Nature, vol. 576, pp. 65–69, Dec. 2019.
  * [13] M. P. Hokmabadi, A. Schumer, D. N. Christodoulides, and M. Khajavikhan, “Non-hermitian ring laser gyroscopes with enhanced sagnac sensitivity,” Nature, vol. 576, pp. 70–74, Dec. 2019.
  * [14] H. Wang, Y.-H. Lai, Z. Yuan, M.-G. Suh, and K. Vahala, “Petermann-factor sensitivity limit near an exceptional point in a brillouin ring laser gyroscope,” Nature Communications, vol. 11, Mar. 2020.
  * [15] J. Wiersig, “Prospects and fundamental limits in exceptional point-based sensing,” Nature Communications, vol. 11, May 2020.
  * [16] S. Scheel and A. Szameit, “PT-symmetric photonic quantum systems with gain and loss do not exist,” EPL (Europhysics Letters), vol. 122, p. 34001, 06 2018.
  * [17] L. Xiao, X. Zhan, Z. H. Bian, K. K. Wang, X. Zhang, X. P. Wang, J. Li, K. Mochizuki, D. Kim, N. Kawakami, W. Yi, H. Obuse, B. C. Sanders, and P. Xue, “Observation of topological edge states in parity-time-symmetric quantum walks,” Nature Physics, vol. 13, no. 11, pp. 1117–1123, 2017\.
  * [18] J. Li, A. K. Harter, J. Liu, L. de Melo, Y. N. Joglekar, and L. Luo, “Observation of parity-time symmetry breaking transitions in a dissipative Floquet system of ultracold atoms,” Nature Communications, no. 10, p. 855, 2019.
  * [19] M. Naghiloo, M. Abbasi, Y. N. Joglekar, and K. W. Murch, “Quantum state tomography across the exceptional point in a single dissipative qubit,” Nature Physics, vol. 15, pp. 1232 – 1236, 2019.
  * [20] U. Günther and B. F. Samsonov, “Naimark-Dilated PT-Symmetric Brachistochrone,” Physical Review Letters, vol. 101, no. 23, p. 230404, 2008.
  * [21] L. Xiao, K. Wang, X. Zhan, Z. Bian, K. Kawabata, M. Ueda, W. Yi, and P. Xue, “Observation of critical phenomena in parity-time-symmetric quantum dynamics,” pp. 1–8, 2018.
  * [22] Y. Wu, W. Liu, J. Geng, X. Song, X. Ye, C.-K. Duan, X. Rong, and J. Du, “Observation of parity-time symmetry breaking in a single-spin system,” Science, vol. 364, no. 6443, pp. 878–880, 2019.
  * [23] H.-K. Lau and D. F. V. James, “Proposal for a scalable universal bosonic simulator using individually trapped ions,” Phys. Rev. A, vol. 85, p. 062329, Jun 2012.
  * [24] C. Shen, Z. Zhang, and L.-M. Duan, “Scalable implementation of boson sampling with trapped ions,” Phys. Rev. Lett., vol. 112, p. 050504, Feb 2014.
  * [25] X. Gu, A. F. Kockum, A. Miranowicz, Y. xi Liu, and F. Nori, “Microwave photonics with superconducting quantum circuits,” Physics Reports, vol. 718-719, pp. 1 – 102, 2017.
  * [26] F. Klauck, L. Teuber, M. Ornigotti, M. Heinrich, S. Scheel, and A. Szameit, “Observation of PT-symmetric quantum interference,” Nature Photonics, vol. 13, pp. 883–887, Sept. 2019.
  * [27] J. Carolan, C. Harrold, C. Sparrow, E. Martín-López, N. J. Russell, J. W. Silverstone, P. J. Shadbolt, N. Matsuda, M. Oguma, M. Itoh, G. D. Marshall, M. G. Thompson, J. C. F. Matthews, T. Hashimoto, J. L. O’Brien, and A. Laing, “Universal linear optics,” Science, vol. 349, no. 6249, pp. 711–716, 2015.
  * [28] L. Jin and Z. Song, “Solutions of PT -symmetric tight-binding chain and its equivalent Hermitian counterpart,” Physical Review A - Atomic, Molecular, and Optical Physics, vol. 80, no. 5, pp. 1–7, 2009.
  * [29] Y. N. Joglekar, D. Scott, M. Babbey, and A. Saxena, “Robust and fragile $\mathcal{PT}$-symmetric phases in a tight-binding chain,” Phys. Rev. A, vol. 82, p. 030103, Sep 2010.
  * [30] P. R. Halmos, “Normal dilations and extensions of operators,” Summa Bras. Math., vol. 2, p. 125, 1950.
  * [31] K. Wang, Y. Shi, L. Xiao, J. Wang, Y. N. Joglekar, and P. Xue, “Experimental realization of continuous-time quantum walks on directed graphs and their application in pagerank,” Optica, vol. 7, pp. 1524–1530, Nov 2020.
  * [32] C. Sparrow, E. Martín-López, N. Maraviglia, A. Neville, C. Harrold, J. Carolan, Y. N. Joglekar, T. Hashimoto, N. Matsuda, J. L. O’Brien, D. P. Tew, and A. Laing, “Simulating the vibrational quantum dynamics of molecules using photonics,” Nature, vol. 557, no. 7707, pp. 660–667, 2018.
  * [33] P. Peng, W. Cao, C. Shen, W. Qu, J. Wen, L. Jiang, and Y. Xiao, “Anti-parity–time symmetry with flying atoms,” Nature Physics, vol. 12, no. 12, pp. 1139–1145, 2016.
  * [34] M. Reck, A. Zeilinger, H. J. Bernstein, and P. Bertani, “Experimental realization of any discrete unitary operator,” Phys. Rev. Lett., vol. 73, pp. 58–61, Jul 1994.
  * [35] Z. Y. Ou and L. Mandel, “Violation of bell’s inequality and classical probability in a two-photon correlation experiment,” Phys. Rev. Lett., vol. 61, pp. 50–53, Jul 1988.
  * [36] K. Kawabata, Y. Ashida, and M. Ueda, “Information retrieval and criticality in parity-time-symmetric systems,” Phys. Rev. Lett., vol. 119, p. 190401, Nov 2017.
  * [37] Y. Lee, M. Hsieh, S. T. Flammia, and R. Lee, “Local $\mathcal{P}\mathcal{T}$ symmetry violates the no-signaling principle,” Physical Review Letters, vol. 112, no. 13, p. 130404, 2014.
  * [38] J. Tang, Y. Wang, S. Yu, D. He, J. Xu, B. Liu, G. Chen, Y. Sun, K. Sun, Y. Han, C. Li, and G. Guo, “Experimental investigation of the no-signalling principle in parity-time symmetric theory using an open quantum system,” Nature Photonics, vol. 10, p. 642, 2016.
  * [39] Y. Aharonov and L. Susskind, “Charge superselection rule,” Physical Review, vol. 155, p. 1428, March 1967.
  * [40] E. Schrödinger, “Über die kräftefreie bewegung in der relativistischen quantenmechanik,” Preuss. Akad. Wiss. Phys.-Math. Kl., vol. 24, p. 418, 1930.
  * [41] O. Klein, “Die reflexion von elektronen an einem potentialsprung nach der relativistischen dynamik von dirac,” Z. Phys., vol. 53, p. 157, 1929.
  * [42] L. J. LeBlanc, M. C. Beeler, K. Jiménez-García, A. R. Perry, S. Sugawa, R. A. Williams, and I. B. Spielman, “Direct observation of zitterbewegung in a bose-einstein condensate,” New J. Phys., vol. 15, no. 7, p. 073011, 2013\.
  * [43] R. Gerritsma, G. Kirchmair, F. Zähringer, E. Solano, R. Blatt, and C. Roos, “Quantum simulation of the dirac equation,” Nature, vol. 463, no. 7277, p. 68, 2010.
  * [44] For continuous-variable systems, the coordinate operator is defined as $\hat{x}=(\hat{a}^{\dagger}+\hat{a})/\sqrt{2}$ where $\hat{a}^{\dagger}$ and $\hat{a}$ are the particle creation and annihilation operators. In turn, for discrete-variable systems, one can write down Pauli $\hat{\sigma}^{x}$ matrix as $\hat{\sigma}^{x}=(\hat{S}^{+}+\hat{S}^{-})/\sqrt{2}$, where $\hat{S}^{\pm}$ are Pauli ladder operators. Since the creation and annihilation operators are similar to ladder operators in their action in the state space, in this context, we consider $\hat{\sigma}^{x}$ as the discrete-variable analog of $\hat{x}$.
  * [45] G. Simpson and H. Mayer-Hasselwander, “Bootstrap sampling - Applications in gamma-ray astronomy,” Astronomy and Astrophysics, vol. 162, pp. 340–348, July 1986.
  * [46] S. Aaronson and A. Arkhipov, “The computational complexity of linear optics,” Theory of Computing, vol. 9, no. 4, pp. 143–252, 2013.
  * [47] J. J. Renema, V. Shchesnovich, and R. Garcia-Patron, “Boson sampling with linear loss is classically simulable,” in Quantum Information and Measurement (QIM) V: Quantum Technologies, p. F3A.3, Optical Society of America, 2019.
  * [48] C. M. Bender, D. C. Brody, and H. F. Jones, “Complex extension of quantum mechanics,” Phys. Rev. Lett., vol. 89, p. 270401, Dec 2002.
  * [49] Y. Ashida, S. Furukawa, and M. Ueda, “Parity-time-symmetric quantum critical phenomena,” Nature Communications, vol. 8, 2017.
  * [50] B. Dóra and C. P. Moca, “Quantum quench in $\mathcal{P}\mathcal{T}$-symmetric luttinger liquid,” Phys. Rev. Lett., vol. 124, p. 136802, Mar 2020.