Experimental Simulation of Symmetry-Protected Higher-Order Exceptional Points with Single Photons

Exceptional points (EPs) of non-Hermitian (NH) systems have recently attracted increasing attention due to their rich phenomenology and intriguing applications. Compared to the predominantly studied second-order EPs, higher-order EPs have been assumed to play a much less prominent role because they generically require the tuning of more parameters. Here we experimentally simulate two-dimensional topological NH band structures using single-photon interferometry, and observe topologically stable third-order EPs obtained by tuning only two real parameters in the presence of symmetry. In particular, we explore how different symmetries stabilize qualitatively different third-order EPs: the parity-time symmetry leads to a generic cube-root dispersion, while a generalized chiral symmetry implies a square-root dispersion coexisting with a flat band. Additionally, we simulate fourfold degeneracies, composed of the non-defective twofold degeneracies and second-order EPs. Our work reveals the abundant and conceptually richer higher-order EPs protected by symmetries and offers a versatile platform for further research on topological NH systems.

The simplest case of EPs is the second-order EP, that is, two-fold degeneracy, which intuitively occurs in a two-dimensional (2D) system, but can be promoted to knotted exceptional lines in three dimensions [26][27][28][29]. Generally, an nth-order EP is stable in a (2n − 2)dimensional NH system [30], which makes them qualitatively more abundant than degeneracies in Hermitian systems. When exploring the NH systems with symmetries, the dimension for the occurrence of generic secondorder EPs can be further reduced from two to one, leading to the observation of stable second-order EPs even in one-dimensional systems [31][32][33][34]. Similarly, generic NH symmetries have been found to reduce the dimension for the occurrence of the third-order EPs from four to two [35][36][37]. Moreover, different symmetries may also entail qualitatively different phenomenology for EPs with the same order [35]. Although higher-order EPs have been studied by a number of works theoretically, their experimental realization and direct observation appear to be rather difficult, as ingenious designs are required to simulate the NH dynamics [38,39] and explore the band structures.
Here, by using single-photon interferometry, we overcome the difficulty of building NH systems with a large number of tunable parameters, and simulate a 2D NH system in the reciprocal space. With interferometric measurements, complex eigenenergies are measured directly, enabling us to construct the band structures of the 2D NH system. In particular, we experimentally observe and characterize two distinct types of symmetry-protected third-order EPs. We also experimentally confirm the stability of the third-order EPs with respect to perturbations. Because they are protected by either the parity-time (PT) or a generalized chiral P symmetry, these EPs disappear upon introducing symmetry-breaking perturbations. Our experimental results demonstrate that the energy near the third-order EP exhibits a generic ∼ k 1/3 dispersion enforced by PT symmetry. The coexistent exceptional ring (ER) composed by the second-order EP bounds an open Fermi surface, which is also the boundary between the PT-unbroken and PT-broken regimes. By contrast, an anomalous ∼ k 1/2 dispersion is observed away from the chiral-P-symmetry-protected third-order EPs.
We further extend our experiment to simulate 2D NH models with four bands, where PT and P symmetries have radically different implications. For the PT symmetry, the third-order EPs and ER are present regardless of the additional bands. In contrast, for the P-symmetry case, the third-order EPs are forbidden by the additional bands. We also observe fourfold degeneracies, each composed by a non-defective twofold degeneracy and a second-order EP [40,41]  tally expose the abundance of higher-order EPs whose codimensions are reduced in the presence of NH symmetries, thus offering potential applications in efficient device design.

Results
Experimental setup. We explore the symmetryprotected third-order EPs by simulating the dynamics of the corresponding three-band NH Bloch Hamiltonians H(k). The complex eigenenergies are measured via interferometric measurements [26].
As illustrated in Fig. 1, our experiment involves three stages: state preparation, nonunitary evolution, and measurement. The basis states of the three-band system are encoded into the hybrid polarization-spatial modes of single photons. In the preparation stage, the states of single photons are initialized in the eigenstate |ψ j (k)⟩ of H(k) with the corresponding eigenenergy E j . After passing through a non-polarizing beam splitter (NPBS), the transmitted photons go through the nonunitary evolution governed by H(k) and the reflected photons remain unchanged. They interfere at the second PBS for interferometric measurements. To circumvent the quantum limit on achieving gain for single photons in the nonunitary-evolution stage, we map the time-evolution operator U = e −iH(k)τ to a dissipative oneŨ = U/ √ Λ. Here we fix the evolution time τ = 1 and ℏ = 1, and take Λ = max j |ζ j |, where ζ j is the eigenvalue of e −iH e iH † [26,42,43]. The corresponding effective NH Hamiltonian for the mapped dissipative system is thus H(k) = H(k) + i ln 1/Λσ 0 , where σ 0 is a 3 × 3 identity matrix. It follows thatH(k) and H(k) have the same eigenstates, and their eigenenergies are related and sat-isfyẼ j = E j − i ln √ Λ. In the measurement stage, the complex inner products of ξ = e −iẼj = ⟨ψ j |Ũ |ψ j ⟩ are obtained through the interferometric measurements. The corresponding eigenenergies are then calculated from the experimentally measured ξ (see Materials and Methods).
As illustrated in Fig. 2B, at ϵ = 0.5, the degenerate point splits into two third-order EPs, which occur at {k x = −0.5929, k y = −0.5694} and {k x = 0.4701, k y = 0.6241}, respectively. Focusing on one of the third-order EPs, we fix k x = −0.5929 and sample 11 different values of k y . The measured real and imaginary parts of the   eigenenergies are shown in the right column of Fig. 2B. By fitting the power exponents with the formula ∼ k β , we confirm that the real eigenenergy exhibits a generic cube-root dispersion (∼ k 1/3 ) near the third-order EP.
Furthermore, an ER emerges in the momentum space, which is particularly interesting as the ER signals the PT transitions, and bounds an open Fermi surface [35]. Because the ER is a collection of second-order EPs, we reveal its existence in Fig. 2C, where we measure the eigenenergies of H P T along the lines of k x = 0 and k y = 0, respectively. The second-order EPs are observed as the bifurcations of the real eigenenergies. In between the two EPs, the imaginary parts of all the eigenenergies are nonzero, indicating the PT-symmetry broken region. Therein, two of the measured eigenergies are approximately equal, indicating the emergence of an open Fermi surface.
As a feature of symmetry protection, the existence of the third-order EPs is robust against symmetrypreserving perturbations. As shown in Fig. 2D, under a small PT-symmetry-preserving perturbation iπ(λ 4 + λ 5 )/20, both third-order EPs persist, but are shifted in parameter space. By contrast, the EPs disappear as illustrated in Fig. 2E, by introducing a general, symmetrybreaking perturbation 8 i=1 δ i λ i (δ i ∈ [−π/20, π/20]). The perturbed spectrum exhibits branch cuts that are π 0 3π 2 2π 2π H P (ε = 0.5) with symmetry-preserving perturbation H P (ε = 0.5) with symmetry-breaking perturbation terminated by paired second-order EPs [41]. These paired second-order EPs further reveal a transition from regions where the real parts of the eigenenergies are degenerate, to those with degeneracy in the imaginary parts.
P-symmetry protected third-order EPs. To demonstrate the features of NH model with P symmetry, we consider a NH Lieb lattice [35,46] which satisfies the relation H P = −P H P P −1 . Note that we consider the case where this symmetry acts locally in momentum space. For ϵ = 0, H P corresponds to the Hermitian Lieb-lattice model. In the momentum space, similar to the case of H P T , a triple degeneracy of the eigenenergy exists at k x = k y = π, from which emerges two dispersive bands with linear scaling [see Fig. 3A]. We introduce the NH term by setting ϵ = 0.5. The triple degeneracy then splits into four third-order EPs that locate at (k x = ±2.6362, k y = ±2.6362). Focusing on one of the EPs, we fix k x = 2.6362 and sample 11 different k y . In Fig. 3B, we show the measured real and imaginary components of the eigenenergies, where the third-order EP is visible near k y = 2.6362. Different from the PT-symmetry-protected case, the two nonzero real bands exhibit an anomalous square-root dispersion near the third-order EP, that is, ∼ k 1/2 . This is a typ- ical behavior associated with the second-order EPs. As shown in Fig. 3C, we introduce a P-symmetry-preserving perturbation in the form of iπλ 1 /20. The four thirdorder EPs persist, with a small shift in the parameter space. Similar to H P T , third-order EPs are destroyed by the general, symmetry-breaking perturbation Accompanying the P-symmetry-protected third-order EPs of H P , there are arc-like degeneracies along the lines of k x = ±k y , which correspond to the real and imaginary Fermi arcs with Re[E ± ]=0 and Im[E ± ]=0, respectively [35]. Experimentally, we again sample 11 different parameters along the line of k x = k y for H P with ϵ = 0.5. In Fig. 4A, we show the measured real and imaginary components of the bands, where both the real and imaginary Fermi arcs can be observed. The boundaries between them give the locations of the third-order EPs. With increasing ϵ, as illustrated in Figs. 4B and C, the third-order EPs move from the center of the Brillouin zone near (k x = π, k y = π) to the edges. At ϵ = 2.0, the imaginary Fermi arcs give way to real Fermi arcs, as pairs of third-order EPs recombine, giving rise to linear dispersions at the EPs. Further increasing ϵ beyond 2, as shown in Fig. 4D with ϵ = 2.1, the third-order EPs completely disappear, as the spectrum opens up complex gaps. The real Fermi arcs develop into closed line-degeneracies with Re[E ± ]=0, slicing through the entire Brillouin zone.
Compared to H P T in Eq. 1, the eigenspectrum of H ′ P T possesses an additional flat band at E = α. As shown in Fig. 5A, third-order EPs persist under α = 1 and ϵ = 0.5, at the same locations as those of H P T . A manifest difference is the emergence of an additional ER composed by the second-order EPs in the PT-symmetric four-band model. In Fig. 5B, by sampling the momentum of k x along the line of k y = 0, we experimentally confirm the existence of the extra real eigenenergy and the additional ER.
For P-symmetry-protected four band models, the restricted NH Hamiltonian becomes where β and γ are arbitrary complex numbers.
Here Hamiltonian H ′ P satisfies the relation with E ′ ± = ± 4 + βγ − 2ϵ 2 + (2 − 2iϵ) cos k x + (2 + 2iϵ) cos k y , where the third-order EPs are destroyed by the additional band. The two degenerate flat bands at zero energy correspond to the non-defective degeneracies with two distinct eigenstates [40,41]. Setting ϵ = 0.5 and tuning k x and k y , we confirm that the two dispersive bands E ′ ± touch each other at E ′ ± = 0, to form the second-order EPs, where the local degeneracies become four-fold (see Figs. 5C and D).
Such four-fold degeneracies are not protected by the P symmetry. Changing the P-symmetry operator to P ′′ = diag(1, −1, 1, −1) [35], we choose the non-Hermitian Hamiltonian in the form As shown in Fig. 5E, the four-fold degeneracies are gapped out, with β = γ = sin k x and ϵ = 0.5.

Conclusion
By experimentally studying the symmetry-protected higher-order EPs of a series of 2D NH models, our work provides the first experimental confirmation that symmetries qualitatively enrich the phenomenology of higherorder degeneracies in the NH systems. By revealing the abundance of higher-order EPs whose codimensions are reduced in the presence of NH symmetries, our results have rich implications for efficient device design.
EPs are generic features of the effective description of a vast range of complex systems ranging from mechanical systems to strongly interacting quantum materials. In contrast, our simple experimental scheme is highly controllable, scalable, and can be readily extended to the detailed study of NH models with arbitrary design and physical origin. It thus constitutes a versatile and efficient platform for the systematic exploration of eigenspectrum and dynamical properties in NH settings.
Materials and Methods Initial state preparation. As illustrated in Fig. 1, our experimental setup can be used to prepare arbitrary pure qutrit states. The basis states of the three-band systems are encoded by the polarizations and spatial modes of single photons, with |0⟩ ⇔ |H 1 ⟩, |1⟩ ⇔ |H 2 ⟩, and |2⟩ ⇔ |V 2 ⟩. Here the subscripts denote the different spatial modes and H (V ) denotes the horizontal (vertical) polarization of single photons. To prepare the generic qutrit state a |H 1 ⟩ + be iφ1 |H 2 ⟩ + ce iφ2 |V 2 ⟩, the photons are first initialized to the horizontal polarization by passing them through a PBS. The two spatial modes of photons are introduced after the photons pass through a beam displacer. The real coefficients {a, b, c} with a 2 +b 2 +c 2 = 1, can be adjusted by controlling HWPs with the setting angles of H 1 = (arcsin a)/2 and H 3 = (arctan c/b)/2. As for the relative phases {φ 1 , φ 2 } ranging from 0 to 2π, they can be adjusted by the sandwich-type set of HWP and QWPs, denoted as QWP-HWP-QWP. Setting the angles of QWPs at 45 • and HWP at 45 • − 180 • φ/π, one can achieve the phase operation of diag(e 2iφ , e −2iφ ). Thus, by setting H 2 = 45 • (φ 1 + φ 2 )/π + 56.25 • and H 4 = 45 • (φ 2 − φ 1 ) + 11.25 • , we can tune the relative phases accordingly.
It follows that the initial state can be exactly prepared in one of the eigenstates of the NH Hamiltonian with prior information. Even if the eigenstates are unknown, we can still generate them as initial states by maximizing the probability of P = | ⟨Ψ|Ũ |Ψ⟩ | 2 / ⟨Ψ|Ũ †Ũ |Ψ⟩. If and only if |Ψ⟩ is one of the eigenstates of the NH Hamiltonian with its dissipative nonunitary unit-time evolution ofŨ , P is maximized to 1. Thus, tuning the angles of the waveplates in state preparation to maximize the measured P , the prepared state must tend to one of the target eigenstates (see section S3 and the Supplementary Materials for more details)).
Experimental implementation ofŨ . To implement the nonunitary operationŨ , we first decompose it intoŨ = V D 3 W using singular value decomposition [47]. The operations of V and W are unitary and We further decompose the two unitary matrices W and V by the established methods in [48,49], i.e., W = W 23 W 13 W 12 and V = V 12 V 13 V 23 , where W ij and V ij are the unitary operations acting on the 2D subspaces of the qutrit system, with the complementary subspace unchanged. As shown in Fig. 1, the unitary operations of W ij and V ij are realized by recombining the photons into the certain spatial mode depending on their polarizations and applying a 2×2 transformation via waveplates. The diagonal matrix D 3 can be realized by introducing the mode-selective losses of photons, where the photon losses can be controlled by the HWPs with the setting angles of H 5 = (arccos ν)/2 and H 6 = (arccos µ)/2.
Interferometric measurements. As shown in Fig.  1, the single photons are generated by spontaneous parametric down-conversion, and prepared in one of the eigenstates |ψ j ⟩ of H(k) with the corresponding eigenenergies of E j . After photons pass through the NPBS, their state is (|t⟩ |ψ j ⟩ + |r⟩ |ψ j ⟩) / √ 2, where t and r denote the transmitted and reflected modes of the single photons, respectively. Applying the nonunitary operation governed byH(k) on the photons in t, the state evolves to e −iẼj |t⟩ |ψ j ⟩ + |r⟩ |ψ j ⟩ / √ 2. After this stage, the interferometric measurements are performed to obtain the overlap between the states of the photons in the transmitted and reflected modes, and thus to extract the complex phase shift ofẼ j .
In the measurement stage, an HWP at 45 • is first applied on the polarization of the photons in the transmitted mode. After the photons pass through a PBS, their state evolves into Here t ′ 1(2) (r ′ 2 ) denotes the transmitted (reflected) mode of the photons which are reflected by the NPBS after passing through the second PBS. Projective measurements are performed on the polarization of the photons in different spatial modes with the bases of {|±⟩ = (|H⟩ ± |V ⟩) corresponds to the spatial mode. Then, we have where N tot denotes the number of the input photons to achieve the state preparation. Note that our setup can be readily extended to simulate the dynamic evolution and construct the energy spectrum of arbitrary NH Hamiltonians, by taking advantage of the extendable degrees of freedom of photons and specially designed interferometric network (see section S5 and the Supplementary Materials for more details).
and2019.0068), and the Göran Gustafsson Foundation for Research in Natural Sciences and Medicine. Author contributions: K.W. performed the experiments with contributions from L.X. P.X. designed the experiments, analyzed the results, and wrote part of the paper. E.J.B. and W.Y. developed the theoretical aspects and revised the paper. H.L. participated in the discussion and revised the paper. Competing interests: The authors declare that they have no competing interests. Data and materials availability: All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials.

Supplementary Materials for "Experimental Simulation of Symmetry-Protected
Higher-Order Exceptional Points with Single Photons" S1 PT-symmetry-protected third-order exceptional points.
In our experiment, we first consider the three-band linearized, higher-spin Dirac-like non-Hermitian (NH) Hamiltonians with parity-time (PT) symmetry in two-dimensional (2D) reciprocal space [35], . (S1) Here λ i with i = 1, 2, · · · , 8 denote the Gell-Mann matrices [44] , i.e., For the spectrum of H P T , all the three eigenenergies can be uniquely defined as The third-order EPs occur when p = q = 0 is taken. The exceptional ring (ER) is composed by second-order EPs and exists on the curve of p 3 + q 2 = 0, which specifies the PT transitions, with all the eigenvalues being real only in the exact PT region with p 3 + q 2 > 0. For p 3 + q 2 < 0, PT-symmetry is broken with E 2 = E * 3 . In Figs. 2A and B of the main text, we characterize the energy dispersion along k y for H P T with ϵ = 0 and ϵ = 0.5. As compensation, we fix k y = 0 with ϵ = 0 and k y = −0.5694 with ϵ = 0.5 for H P T with various k x experimentally. As shown Fig. S1A, the measured energy spectrum for the Hermitian model with ϵ = 0 presents a triple degeneracy  with the same linear dispersion along k x . The PT-symmetric NH model with ϵ = 0.5 has the energy spectrum which exhibits a cube-root dispersion near the third-order EP along k x as well [see Fig. S1B]. Additionally, we also study the energy dispersion near the other third-order EP for H P T with ϵ = 0.5. The same features can be derived from the results shown in Figs. S1C and D, where a cube-root energy dispersion away from the third-order EPs located at (k x = 0.4701, k y = 0.6341) is observed experimentally.
As shown in Fig. 2E of the main text, the third-order EPs are destroyed by the symmetry-broken perturbation We consider the P-symmetry-protected NH Lieb lattice model [35,46], with eigenenergies 0 and E ± = ± 2(2 + cos k x + cos k y − iϵ cos k x + iϵ cos k y − ϵ 2 ). For ϵ = 0, H P describes the standard nearest neighbor Hermitian model, which possesses a triple energy degeneracy at k x = k y = π. For ϵ ̸ = 0, the triple energy degeneracy splits into four third-order EPs, which are located at (k x = ± arccos(ϵ 2 /2 − 1), k y = ± arccos(ϵ 2 /2 − 1)). In Figs. 3A and B of the main text, we measure the energy dispersion along k y by fixing k x = π for H P with ϵ = 0, and k x = 2.6362 for H P with ϵ = 0.5. As shown in Figs. S2A and B, by fixing k y = π for H P with ϵ = 0 and k y = 2.6362 for H P with ϵ = 0.5, respectively, we sample 11 k x in our experiment. The measured energy spectrum for H P with ϵ = 0 presents a triple degeneracy with the linear dispersion along k x . For H P with ϵ = 0.5, an anomalous square-root scaling away from the third-order EP along k x is also observed. In addition, the energy dispersions near the other third-order EP at (k x = −2.6362, k y = −2.6362) are shown in Figs. S2C and D. The similar square-root energy dispersion is observed in our experiment.
As shown in Fig. 3D of the main text, the third-order EPs are destroyed by the symmetry-broken perturbation

�90°90°1
.0  As demonstrated in the main text, given an approximate eigenstate, our setup can be used to efficiently estimate the corresponding eigenvalue. However, to calculate the eigenstates by using conventional classical methods requires resources. When the size of the system grows, rapidly increasing of the resources becomes an obstacle. Furthermore, finding the eigenstates of a given Hamiltonian is a fundamental problem and has lots of applications. Here we develop a method of search for the eigenstates of the NH Hamiltonian, which can also be applied to the Hermitian Hamiltonian, apparently.
The search proceeds by preparing the input state of |Ψ⟩. The state is then evolved through a unit-time evolutioñ U = e −iH governed by the NH HamiltonianH. By projecting the evolved state to the initial state, we can obtain the normalized probability of finding the output state unchanged as P = | ⟨Ψ|Ũ |Ψ⟩ | 2 / ⟨Ψ|Ũ †Ũ |Ψ⟩. For the evolution operator, it can be rewritten asŨ = i e −iẼi |ψ i ⟩⟨χ i | /⟨χ i |ψ i ⟩. Here ⟨χ i | and |ψ i ⟩ are the left and right eigenstates of H with the corresponding eigenenergy ofẼ i , which satisfy the relations of For the input state of |Ψ⟩, it can be represented as where c i = ⟨χ i |Ψ⟩ and b i = ⟨ψ i |Ψ⟩. Thus, we have Subtracting the numerator from the denominator, we can obtain where |Ψ ⊥ l ⟩ are the states orthogonal to the input state of |Ψ⟩. Thus, we can obtain P ≤ 1. If and only if |Ψ⟩ is an egigenstate ofH, the inequality is saturated to an equality with P = 1. Therefore, P can act as an eigenstate probe, where it equals to the maximum value of 1 if the input state is an eigenstate of the Hamiltonian.
As proof of principle, we achieve the searching task to prepare the initial state into one of the eigenstates ofH P T with (k x = −0.6529, k y = −0.5694) and ϵ = 0.5. The setup is shown in Fig. S3A, where the input qutrit state is prepared by tuning the setting angles of the half-wave-plates (HWPs) of H 1−4 . Then, subjecting the state to the evolution governed by the measured Hamiltonian, the probability P can be measured by subjecting the output state to the reverse process of the initial state preparation. We then have P = N 1 /(N 1 + N 2 + N 3 ), where N i is the number of heralded clicks at detector D i . By continuously tuning the angles of H 1−4 to achieve the maximum P , we can prepare the input state as one of the eigenstate of the Hamiltonian, where the searching task can be interpreted as an optimization problem.
In our experiment, we first extend the PT-symmetry-protected system to the four-band model governed by the Hamiltonian with the symmetry operators P ′ = diag(1, −1, 1, 1) and T ′ = diag(−1, 1, i, 1). Here Γ 0 is the 4 × 4 identity matrix and Γ i (i = 1, 2, · · · , 15) denote the Gell-Mann matrices, that span the Lie algebra of the SU(4) group, Using the representation of the P-symmetry operator as P ′ , the P-symmetric non-Hermitian Hamiltonian is restricted to the form of with the corresponding eigenstates of Thus, we tune two parameters to make bd + f h + βγ = 0 satisfied. Then the second-order EPs with E ± = 0 are generated, which also correspond to the four-fold degeneracy points for the NH system. As shown in Figs. 5C and D of the main text, we choose b = 1 + e ikx − iϵ, d = 1 + e −ikx − iϵ, f = 1 + e iky + iϵ, h = 1 + e −iky + iϵ and β = γ = sin k x . By setting ϵ = 0.5 and tuning k x and k y , the two dispersive bands touch each other to form the second-order EPs at {±2.7275, ±2.7275}, where the four-fold degeneracies emerge. As shown in Fig. S4, to construct the energy spectrum of the symmetry-protected four-band models of H ′ P T and H ′ P in the main text, the basis states are encoded by the hybrid polarization-spatial modes of the single photons as |0⟩ ⇔ |H 1 ⟩, |1⟩ ⇔ |V 1 ⟩, |2⟩ ⇔ |H 2 ⟩, and |3⟩ ⇔ |V 2 ⟩. Here 1, 2 denote the different spatial modes and H (V ) denotes the horizontal (vertical) polarization of the single photons. The state preparation is achieved by passing the photons through a polarizing beam splitter (PBS). An adjustable HWP (H 1 ) combined with a sandwich-type set of HWP (H 2 ) and quarter-wave plates (QWPs) with the setting angles of 45 • are used to control the amplitude and the relative phases of the photon with different polarizations. After passing through a beam displacer (BD), the vertically polarized photons are transmitted and the horizontal polarized photons go through a 3-mm lateral displacement into a neighboring mode. Thus, the photons are prepared into two parallel spatial modes-the transmitted first and lateral second modes. By inserting the wave-plates H 1−6 into the spatial modes accordingly, we prepare the state of the photons in the eigenstates of the measured Hamiltonian.
After the state preparation, the photons are injected into the 50 : 50 non-polarizing beam splitter (NPBS). Thus, extra two paths of the transmitted (t) and reflected (r) modes are introduced. To approach the nonunitary unit-time evolution governed by the mapped NH Hamiltonian on the transmission, we decompose the nonunitary operation as HereH ′ denotes the mapped NH Hamiltonian ofH ′ P T orH ′ P . The operations V ij and W ij are the unitary operations acting on the 2D subspaces of the system with the complementary subspace unchanged, and D 4 is a diagonal matrix. As shown in Fig. S4, the unitary operations V ij and W ij can be realized by combining the two acted modes into one spatial mode with different polarizations and applying a 2 × 2 unitary transformation via the set of wave-plates. The diagonal matrix D 4 can be realized by introducing the mode-selective losses to the corresponding modes, where the intensity of the losses can be controlled by HWPs of H 7−9 .
After recombining the photons in the t-and r-modes, the projective measurements are performed on the polarization of the photons in different spatial modes with the bases of {|±⟩ = (|H⟩ ± |V ⟩) / √ 2, |R⟩ = (|H⟩ − i |V ⟩) / √ 2}. The coincidences for the projective measurements are counted as {N ± i , N R i }. Here i ∈ {t ′ 1 , r ′ 1 , t ′ 2 , r ′ 2 } denotes the transmitted (t ′ ) or reflected (r ′ ) path following the second PBS in the first or second spatial modes of the r-path photons. Accordingly, we can obtain the eigenenergy ofH ′ as the complex phase shift between t and r paths, i.e., where N tot denotes the number of the input photons to achieve the state preparation. The eigenenergy E ′ j of H ′ P T or H ′ P can thus be obtained through the relation of E ′ j = i ln ξ ′ + i ln √ Λ ′ , where Λ ′ = max k |ζ k | and ζ k is the eigenvalue of e −iH ′ P T e −iH ′ † P T or e −iH ′ P e −iH ′ † P .
S5 Generalization to arbitrary models.
Our setup can be efficiently generalized to construct the energy spectrum of arbitrary NH Hamiltonians, by taking advantage of the extendable degrees of freedom of photons and specially designed interferometric network.
To extend the proposed method to achieve the preparation of the initial state, it is convenient to expand the spatial modes by additional BDs and tune the parameters by using the wave-plates. Actually, for an N -dimensional pure states, it can be prepared by using (N − 1)/2 BDs when N is an odd number and N/2 − 1 for an even number. All the proportions of each mode and the relative phases can be tuned by varying the angles of N − 1 HWPs and N − 1 sandwich-type sets of QWP-HWP-QWP.
After passing the initial state through the nonunitary unit-time evolution governed by the mapped NH Hamiltonian, the corresponding eigenenergies will be introduced as the complex phase shift to the evolved state. The decomposition method presented in the main text to achieve the nonunitary dynamics ofŨ can also be generalized to any N × N passive nonunitary operation asŨ N = n=1 n=N m=1 m=n−1 V mn D N i=N i=1 j=N j=i+1 W ij . Thus, the simple building blocks of BDs and wave plates to approach the corresponding decomposed operations provide a recipe for an implementation of arbitraryŨ N , experimentally. The maximum numbers of BDs needed to build all the N -dimensional unitary operators of W ij and V mn are both 2N − 4 when N is an even number and 2N − 3 for an odd number. The number of BDs to achieve D N is N − 1. Thus, The maximum number of BDs to approach this process is only linear increasing with N . Besides, to fully control the parameters in the nonunitary operation, the maximum number of set of wave plates to approach the decomposed unitary operations is (N − 1)N , and N − 1 HWPs to realize D N . At last, by performing interferometric measurements on the photons, the complex eigenenergies for each mode in reciprocal space can be constructed experimentally.