Non-Hermitian Kibble-Zurek mechanism with tunable complexity in single-photon interferometry

Non-Hermitian descriptions of quantum matter have seen impressive progress recently, with major advances in understanding central aspects such as their topological properties or the physics of exceptional points, the non-Hermitian counterpart of critical points. Here, we use single-photon interferometry to reconstruct the non-Hermitian Kibble-Zurek mechanism and its distinct scaling behavior for exceptional points, by simulating the defect production upon performing slow parameter ramps. Importantly, we are able to realise also higher-order exceptional points, providing experimental access to their theoretically predicted characteristic Kibble-Zurek scaling behaviour. Our work represents a crucial step in increasing the experimental complexity of non-Hermitian quantum time-evolution. It thus also furthers the quest to move the frontier from purely single-particle physics towards increasingly complex settings in the many-body realm.

The foundational axioms of quantum mechanics impose a Hermitian structure on Hamiltonians.However, it is now appreciated that rich and unconventional phenomena can arise in settings where the constraints enforced by such Hermitian structures are absent.This happens rather generically for systems in touch with an environment; experimental instances occurring in photonics [1][2][3][4], cold atoms [5,6], mechanical systems [7][8][9], and electric circuits [10][11][12][13] have revealed rich singleparticle properties induced by non-Hermiticity.In particular, the so-called exceptional points (EPs) [14][15][16][17] provide the non-Hermitian [18] counterpart to the familiar critical points [19] in Hermitian band structures [20,21], leading to critical phenomena unique to non-Hermitian systems.Whereas most previous studies focus on singleparticle properties, we are now faced with the challenge of extending their experimental reach, by realising specific phenomena characteristic of non-unitary dynamics on one hand, and allowing the study of more complex problems on the other, shifting the frontier towards the many-body realm.
In this work, we report a progress along both axes using single-photon interferometric networks as an experimental platform.We provide a framework to achieve highly tunable non-Hermitian band structures allowing us to realise different classes of EPs with varying complexity.Mathematically, at an EP, two (or more) complex eigenvalues and eigenstates coalesce [14].These eigenstates then no longer form a complete basis, which in turn feeds into distinct properties of the EP compared to those at Hermitian critical points.We consider various types of non-Hermitian EPs, culminating in a higher-order EP, the counterpart of the familiar higherorder critical points, which generically appear in multidimensional phase diagrams upon tuning a pair of parameters.
Such a framework enables us to systematically characterise the novel non-unitary dynamics pertaining to different types of EPs, which is specifically visible in the dynamics of defect generation upon passage through the exceptional/critical point in the time domain, as captured by the venerable Kibble-Zurek mechanism [22][23][24][25][26][27].We simulate such ramps in a many-body setting, and probe the characteristic universal scaling behavior of the resulting defect density: Here, v denotes the sweep velocity and ν the correlation length exponent.Remarkably, the defect density n obeys a scaling form even in the non-Hermitian case with, however, one crucial difference.While for conventional critical points d * = d is identical to the spatial dimension d, in the non-Hermitian setting the effective dimension [15] d * = d + z contains a shift by the dynamical critical exponent z.Using single-photon interferometric networks we measure n and observe a power-law dependence with exponents in keeping with the Kibble-Zurek prediction both for the Hermitian and the non-Hermitian cases.
With the ability to tune the band structure complexity, we then also address the Kibble-Zurek mechanism for higher-order EPs and determine their characteristic prop-erties experimentally.Our approach is based on simulating sets of independent modes, mimicking the dynamics of effectively noninteracting quasiparticles upon crossing critical and exceptional points.This technique provides a flexible framework for realising and probing non-unitary dynamics with high tunability and control, with the potential to push the frontier further towards the quantum simulation of many-body effects.

Model, observables and protocol
We consider translationally invariant quantum systems which can be described in terms of independent modes: Here, H p denotes a 2 × 2 matrix represented in terms of Pauli operators σ α (α = x, y, z), and parametrized through the momentum p of the mode and the couplings ∆, Γ.This Hamiltonian exhibits genuine non-Hermitian character due to the complex mass term, with which marks the location of a second-order EP.It is the key goal of this work to experimentally study the dynamical consequences of such EPs [28] and to contrast with those of conventional Hermitian critical points [19,27].
Based on the general protocol that we outline below, the non-unitary dynamics can be further enriched by engineering H p as an enlarged 4 × 4 matrix, enabling us to access a mode structure of increased complexity, and the concomitant unconventional higher-order EPs.
For our protocol, we initialise each mode p in its ground state |Ψ p with Γ = 0.At time t = 0, we then start the non-equilibrium process by linearly increasing Γ ∝ t, driving the system either through a critical or exceptional point.The defects in the final state |Ψ p (τ ) at time t = τ are quantified via The total defect density n = (2π) −1 dpσ z (p, τ ) according to the Kibble-Zurek prediction Eq. ( 1) shows universal behavior, when measured with respect to its equilibrium or steady-state value, n eq (see Methods).For ease of notation, we set the velocity of p-modes and to unity, and the density of states to 1/π.

Experimental implementation
Here, the rotation operators R(θ j , ϕ j , ϑ j ) (j = 1, 2) are experimentally implemented using a set of sandwichtype wave plates with setting angles ϕ j , θ j , and ϑ j .The polarisation-dependent loss operator L(θ H , θ V ) contributes the non-unitary dynamics, and is experimentally achieved by a combination of beam displacers (BDs) and HWPs with setting angles φ H and φ V .We choose the parameters {θ j , ϕ j , ϑ j , φ H , φ V } such the sequence of three operations reproduces the targeted U p (τ ).Such a scheme enables us to implement any desired non-unitary operator for two-level systems [21].We measure σ z (p, τ ) by recording the relative photon counts in the basis {|H , |V } through a PBS and avalanche photodiodes (APDs), with a typical peak count of 160, 000 photons.In order to estimate the total defect density n, we utilize Gauss-Legendre quadrature which uniquely determines the specific p points once the integration domain is fixed.

Hermitian Kibble-Zurek scaling
We first validate the set-up by reproducing the wellknown Kibble-Zurek scaling for Hermitian quantum critical points with quantitative accuracy.Concretely, we consider ∆ = 0 and Γ = −i∆ 0 t/τ , yielding a mass term linearly increasing with time.Figure 2a shows the results for σ z (p, τ ) in scaled units so as to achieve a data collapse with the predicted exponents: the resulting total defect density n, see Fig. 2c, is consistent with powerlaw behavior over more than one decade.From a fit n − n eq ∼ τ −α we obtain α = 0.494(2), which, within the error bars, agrees with the theoretical prediction α = 0.5 for the underlying equilibrium quantum critical point of Ising universality class with d = z = 1, ν = 1/2.Here, the constant n eq denotes the expected adiabatic value for σ z (p, τ ) in the τ → ∞ limit.

Non-Hermitian Kibble-Zurek scaling
For the central aspect of this work, the non-Hermitian Kibble-Zurek scaling at EPs, we consider two cases.First, a parity-time (PT )-symmetric ramp with ∆ = ∆ 0 and Γ = ∆ 0 t/τ , see Figs. 2b and c, where the energy spectrum E ± (p) of the system remains real throughout the ramp approaching an EP at t = τ .Second, a fully non-Hermitian drive with ∆ = 0, Γ = Γ 0 t/τ , see Fig. 3.
To reduce experimental error for the fully non-Hermitian drive, we sample the mode-resolved defect densities and analyse the scaling behaviour in regions p ≷ Γ 0 separately.Whereas the difference in the integrated densities between the two regions gives the total integrated density as before, the smallness of such a difference would lead to small photon counts and significantly larger error bars.By analysing the two regions separately, as we    show below, a universal scaling behaviour is established, since both regions respect the same scaling law.
For both cases, the data collapse across the varying ramp times τ , characteristic of the Kibble-Zurek mechanism, but now realised for non-Hermitian systems.The resulting scaling functions, however, differ markedly from the Hermitian case in Fig. 2a.Especially, the full non-Hermitian drive provides a manifestly non-Hermitian feature: the double-peak structure in Fig. 3a implies modes p with lower occupation than for the purely adiabatic limit, which is impossible for the Hermitian case.
While the scaling functions appear to have rather unconventional form, the integrated total defect density n robustly exhibits power-law behaviour, as illustrated in Figs.2c and 3b.The associated exponents follow the modified scaling law in Eq. ( 1) with the effective dimension d * .For the PT -symmetric ramp, the fit n − n eq ∼ τ −α yields an exponent α = 0.68(1), consistent with the theoretical prediction α = 2/3 for the underlying EP with critical exponents d = z = 1, ν = 1/2.For the full non-Hermitian drive, the fitted exponents are both α = 0.94 (7) for the two regions with p ≷ Γ 0 , see Fig. 3b.These suggest that the total defect density should obey a power-law scaling in agreement with Eq. ( 1) for d = ν = 1 and z = 1/2.These results demonstrate the high accuracy with which our experiments can probe the dynamics of non-Hermitian systems.
The increasing experimental error bars for longer evo- (σ z (p,τ)-σ eq z (p))�τΓ  lution times in the non-unitary dynamics, Figs.2c and  3b, are mainly due to the non-unitary photon losses: a smaller number of photons implies a larger statistical error.Importantly, the influence of imperfections appears weaker for the total defect density n than for the moderesolved one σ z (p, τ ).This happens because the errors in different modes are statistically independent, and hence suppressed upon integration: the majority of experimental imperfections originate from wave plates and BDs that are independently tuned in different p sectors.

Kibble-Zurek at higher-order exceptional points
The key next step is to increase the complexity of non-Hermitian Hamiltonians, accessing previously unreachable physical properties.Here, we achieve this by enlarging the mode matrix H p to be of 4 × 4 form, providing an experimental access to a higher-order EP.Here, For Γ = 0, the model is Hermitian, while for Γ = ∆, it is PT -symmetric and features a fourth-order EP at p = 0.Here the spectrum scales as ∼ p 1/4 due to detuning from the EP, while the gap scales ∼ √ Γ − ∆ at p = 0, suggesting critical exponents z = 1/4 and zν = 1/2.
The four basis states are now encoded in the polarisations and spatial modes of single photons, and given by {|U H = (1, 0, 0, 0) T , |U V = (0, 1, 0, 0) T , |DH = (0, 0, 1, 0) T , |DV = (0, 0, 0, 1) T }.Here |U and |D represent, respectively, the upper and lower spatial modes of photons (see Fig. 1b).The experimental implementation of the calculated U (4) p , however, is different from that of U p for two-level systems which is realised exactly  by BDs and wave plates.Here, instead, we approximate U (4) p with a series of modules, each consisting of two BDs and a set of wave plates (see Methods for details).
We characterize the final momentum-resolved defect density through M y (p) − M eq y (p), where for small p to a very good approximation, with M y given by M y = i ∂H (4)   p ∂Γ and |ϕ p the eigenstate of H (4) p with the largest imaginary eigenvalue under the condition Γ = ∆.
The measured M y (p) − M eq y (p) for various τ , again exhibits a scaling collapse, see Fig. 4a.The defect production in the presence of a higher-order critical point exhibits a distinct feature, different from the previously studied cases in Figs.2b and 3a.Specifically, M y (p)−M eq y (p) develops a pseudogap-like feature around p = 0 and increases monotonically with momentum.The total defect density n, Fig. 4b, exhibits power-law behaviour fitted to α = 1.80 (5), close the Kibble-Zurek prediction, α = 5/3, while error bars are larger compared to the lower-order EPs.

Discussion
Our work constitutes the first experimental investigation of non-Hermitian Kibble-Zurek scaling due to nearadiabatic passage across exceptional points.The programmability of our interferometric photon network allows us to engineer both Hermitian and non-Hermitian band structures with high tunability.In particular, we demonstrated that we can also generate multi-band models which has enabled us to access higher-order exceptional points and their distinctive properties.As a key consequence, we have identified a platform for realizing non-Hermitian dynamics of increasing complexity as a route towards the many-body realm with higher dimensional time evolution matrices.

Theoretical background
As detailed in Ref. [15], we consider the time evolution of a (non-)Hermitian Hamiltonian H(t).The initial state is the ground state of the starting Hamiltonian, which is always Hermitian.At time t = 0 we start our timedependent protocol, governed by with |Ψ(t) = ⊗ p |Ψ p (t) for a given mode represented by momentum p, and the time dependence of H(t) comes from Γ(t) and ∆(t), with Γ(t = 0) = 0 always.Under the non-unitary time evolution, the norm of the wave function is not conserved, so an additional prescription for performing measurements in such states has to be given.When interpreting such dynamics as a result of dissipation in the framework of a Lindblad master equation with an additional continuous measurement, expectation values of an operator O have to be evaluated as [29][30][31] where the left state, Ψ(t)| is taken as the Hermitian conjugate of the time evolved right state, |Ψ(t) .Since the initial condition at t = 0 is chosen to be the ground state of a Hermitian system, the initial right and left states also satisfy this condition.In the following, we quantify the defect production via with N denoting the number of considered momentum states.
The τ → ∞ limit of the momentum-resolved defect density coincides [15] with calculating the expectation value of σ z for a given p momentum state using the normalized right eigenfunction of the final Hamiltonian, corresponding to the smallest eigenvalue for the PTsymmetric case, or to the complex eigenvalue with the largest imaginary part otherwise.This defines σ eq z (p).Upon integrating this over momenta, we obtain n eq .
The density of states of Eq. ( 2) with ∆ = Γ = 0 is defined as ) with E ± (p) = ±|p| and δ(E) the Dirac delta function.This becomes energy, E, independent and takes the value quoted in the main text.
Experimental realisation of (non-)unitary time evolution for the two-level system In order to simulate the time evolution of these twolevel systems, we encode the basis states in the horizontal and vertical polarisations of a single photon, with |H = (1, 0) T and |V = (0, 1) T .We generate heralded single photons via type-I spontaneous parametric downconversion, with one photon serving as the trigger and the other as the signal.The signal photon is initialized in the ground state |Ψ p (0) of H p with Γ = 0 via a polarising beam splitter (PBS), a quarter-wave plate (QWP) and a half-wave plate (HWP), with p-dependent parameters.We then send the single photon to the interferometric network as illustrated in Fig. 1.
To simulate the non-unitary dynamics, driven by a time-dependent H(t) up to a time τ , we decompose the dynamics into different momentum sectors, and directly implement the time-evolution operator U p (τ ).Specifi-cally, we first numerically calculate U p (τ ) through where t k = (k − 1/2)δt, δt = τ /N , with N ∈ N. We assume H p (t k ) to be time-independent within each k, and take sufficiently large N , such that Eq. ( 8) converges.
As illustrated in Fig. 1, we implement U p (τ ) according to where the rotation operator R(θ j , ϕ j , ϑ j ) (j = 1, 2) is realised using a sandwich-type wave-plate set, including a HWP at the setting angle ϕ j , and two QWPs at θ j and ϑ j , respectively.The polarisation-dependent loss operator L = 0 sin 2φ V sin 2φ H 0 is realised by a combination of two beam displacers (BDs) and two HWPs with setting angles φ H and φ V .The setting angles {θ j , ϕ j , ϑ j , φ H , φ V } are fixed according to the numerically calculated U p (τ ).We note that Eq. ( 9) enables us to implement arbitrary non-unitary operators for a twolevel system with different setting angles [21].
After performing the time evolution, we measure the expectation value of σ z through projective measurements.Specifically, we measure the probability of photons in the basis {|H , |V } through a PBS and avalanche photodiodes (APDs).The outputs are recorded in coincidence with trigger photons.Typical measurements yield a maximum of 160, 000 photon counts.We then construct the momentum-resolved defect density through where N H and N V are the photon counts with horizontal and vertical polarisations, respectively.

Spectrum of H
(4) p with a fourth-order EP The 4x4 Hamiltonian, H p in Eq. ( 4) can be diagonalized analytically, yielding four bands as with α = ±, β = ±.The evolution of the instantaneous spectrum is depicted in Fig. 5.For Γ = ∆, it reduces to which is dominated by the term 32i∆ 3 p under the double square root for small momenta, thus realising an EP4.As advertised above, this Hamiltonian can also be rewritten in terms of the Pauli matrices of two interacting spins, denoted by τ i and σ i with i = x, y, z.Such a minimal many-body model reads Time evolution of the four-level system.
Similar to the case with two-level systems, we numerically calculate the time-evolution operator U   4)  p (t k )δt , for sufficiently large N .
Here |U and |D represent, respectively, the upper and lower spatial modes of photons (see Fig. 1b).The experimental implementation of the calculated U (4) p , however, is different from that of U p for two-level systems which is realised exactly by BDs and wave plates.Instead, we approximate U (4) p with a series of modules, each consisting of two BDs and a set of wave plates.Specifically, each module features two sets of sandwiched wave plates (QWP-HWP-QWP), and a sandwiched set of BDs and wave plates (BD-HWP-BD).The QWP- where φ M is the setting angle of the HWPs between two BDs.
To estimate the deviation of the implemented nonunitary time-evolution operator U with respect to the ideal U where Ũ := U/Tr U U † .The distance varies between 0 and 1, with 0 indicating a perfect implementation of U p .As shown in Fig. 6, the distance is already below 10 −4 for the experimentally relevant parameters when we use three sets of modules.In principle, we can achieve even smaller distances by increasing the number of modules.
In order to test the Kibble-Zurek scaling around EP4 experimentally, we probe M y (p) via projective measurements of the time-evolve state |Ψ p (τ ) .More i=1 o i P i , where P i is the photon count from the corresponding APD.We note that M eq y is measured in a similar fashion, by performing projective measurements directly on the state |ϕ p .
Error analysis.
The difference between the experimental data and theoretical predictions is mostly caused by the inaccuracy of the wave-plate parameters, as well as the dephasing of BDs in the cascaded interferometric network.Since both wave plates and BDs are tuned independently in different p-sectors, experimental errors in different p-sectors are uncorrelated.
For the data shown in Fig. 2a, the unitary time evolution is realised by three HWPs and without any interferometers (hence no BDs).In contrast, the data shown in Figs.2b, 3a and 4a, the non-unitary time evolutions therein require more wave plates as well as BDs.This is the reason for the more apparent experimental imperfections in the case of non-unitary time evolutions.However, since the imperfections in each p are independent (see discussion above), and the data shown in Figs.2c, 3b and 4b are obtained by integrating over p, the differences between the experimental data and theoretical predictions are smaller for the total defect densities.

FIG. 1 .
FIG.1.Experimental setup.As a first step a photon pair is created via spontaneous parametric down-conversion (SPDC).One of the photons serves as a trigger, the other is projected into the polarisation state |H or |V with a polarising beam splitter (PBS) and a half-wave plate (HWP) before entering the interferometric network.a Experimental setup for the nonunitary dynamics of a single-qubit PT -symmetric system.b Experimental setup for the non-unitary dynamics of a four-level system, where the qudit is encoded in polarisational and spatial degrees of freedom.To prepare the initial state, heralded single photons pass through a PBS and a HWP with tailored setting angles and are split by a birefringent calcite beam displacer (BD) into two parallel spatial modes: upper and lower modes.After passing through wave plates inserted into optical paths of the two spatial modes, the photons are prepared in the initial state |Ψp(0) or |ϕp (see the main text for definitions), through p-dependent wave plates.The non-unitary evolution is then realised by three modules involving BDs and wave plates, and projective measurements are performed by a cascaded interferometer and avalanche photodiodes (APDs).

FIG. 2 .
FIG.2.Kibble-Zurek scaling of defect production for parametric ramps in Hermitian and PT -symmetric non-Hermitian systems.Data collapse for the momentum-resolved defect density for the Hermitian case a and the PT -symmetric ramp b.Dots with error bars display the experimental measurements and solid lines refer to numerical simulations.c Total defect density n relative to the adiabatic values neq for the two considered cases as a function of the ramp time τ .Solid lines show the results for the best power-law fit n − neq ∼ τ −α to the experimental data.The fitted exponents α = 0.494(2) for the Hermitian scaling (red), and α = 0.68(1) for the PT -symmetric ramp (blue).

FIG. 3 .
FIG. 3. Kibble-Zurek scaling of defect production for ramps across non-Hermitian EPs.Data collapse for the momentum-resolved defect density for the full non-Hermitian drive a.Total defect density relative to the adiabatic values as a function of the ramp time τ b.As discussed in the main text, we sample the left (p < Γ0) and right (p > Γ0) regions separately.The red solid (black dashed) line in b shows the results for the best power-law fits |n − neq| ∼ τ −α for the left (right) region, with fitted exponents α = 0.94(7) in both regions.

FIG. 5 .
FIG.5.The evolution of the instantaneous spectrum of H(4) .The real and imaginary parts as a function of Γ are shown in solid and dashed curves, respectively.The evolution starts from the Hermitian case with Γ = 0 a, through the PT -broken regions b and c, and ends at Γ = ∆ d, where the EP4 appears are at p = 0.

( 4 )
p (τ ) at the final time τ for each p-sector.More explicitly, we calculate U (

FIG. 6 .
FIG.6.Distance between the implemented nonunitary operators and the ideal ones.Calculated distance d for all sampled p and τ , with different number of modules in the implementation of U (4) p .We plot the distance for all sampled p with τ ∆ = 2 n in the x-axis range (10n − 40, 10n − 30).
specifically, we implement the transformation O = |U H o 1 | + |U V o 2 | + |DH o 3 | + |DV o 4 | through cascaded interferometers consisting of BDs and wave plates.Here |o i (i = 1, • • • , 4) are eigenstates of M y , with M y = i o i |o i o i |, where o i are the eigenvalues of M y .A PBS is then used to map the basis states {|U H , |U V , |DH , |DV } to four distinct spatial modes, where photons are collected by four APDs.It follows that M y (p) = 4