The Kibble-Zurek mechanism at exceptional points

Exceptional points (EPs) are ubiquitous in non-Hermitian systems, and represent the complex counterpart of critical points. By driving a system through a critical point at finite rate induces defects, described by the Kibble-Zurek mechanism, which finds applications in diverse fields of physics. Here we generalize this to a ramp across an EP. We find that adiabatic time evolution brings the system into an eigenstate of the final non-Hermitian Hamiltonian and demonstrate that for a variety of drives through an EP, the defect density scales as τ−(d + z)ν/(zν + 1) in terms of the usual critical exponents and 1/τ the speed of the drive. Defect production is suppressed compared to the conventional Hermitian case as the defect state can decay back to the ground state close to the EP. We provide a physical picture for the studied dynamics through a mapping onto a Lindblad master equation with an additionally imposed continuous measurement.


Results
The model and observables. We consider Hamiltonians of the form [14][15][16][17] H ¼ which can be decomposed into different momentum sectors labeled by p. For iΓ 2 R; the problem is Hermitian. For Γ 2 R; instead, the above Hamiltonian becomes non-Hermitian with a spectrum given by E ± ðpÞ ¼ ± ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi p 2 þ Δ 2 À Γ 2 p . When Δ > Γ, H p has real eigenvalues for each p. For Γ > Δ on the other hand H p has, in general, complex eigenvalues. At sufficiently large p> ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi Γ 2 À Δ 2 p , however, the spectrum becomes real again. A Hamiltonian is PT-symmetric if it commutes with the combined parity and time reversal operators. Its spectrum is real if PTsymmetry is not spontaneously broken, i.e., the eigenstates also respect PT symmetry. With broken PT-symmetry, the spectrum becomes complex. The analysis of higher-order exceptional points 12 is beyond the scope of the current investigation. Non-Hermitian Hamiltonians of the kind in Eq. (2) can be emulated by optical waveguides 6,10,15 , distributed-feedback structures 18 , microcavities 11 , or electric circuits 19 . Eq. (2) also accounts for the low energy dynamics of the quantum Ising chain in an imaginary transverse field 20,21 or an imaginary mass fermion (i.e., tachyon) system 17 . The last term in Eq. (2) assumes balanced gain and loss without loss of generality: one can shift the diagonal term in the Hamiltonian by any complex value without affecting the results (see Methods).
Hermitian dynamics. We investigate several natural scenarios for time-dependent Δ and Γ. Let us start with the Hermitian Kibble-Zurek mechanism, which is contained in our model by choosing Δ = 0, Γ = −iΔ 0 t/τ. This is the conventional Landau-Zener problem 22,23 starting exactly from the critical point, thus it represents only a half crossing. This yields (see Methods) 〈σ y (τ)〉 = 0, while with W the high energy cutoff. Defect production is effective when the adiabatic condition is violated 24,25 , namely when dln|Γ|/ dt~|Γ|. This gives the transition time τ 1/2 , and the defect density scales inversely with this, in accord with ref. 26 . This also follows from the scaling behaviour of the momentum resolved defect density which is typical for the Hermitian Landau-Zener problem, is a universal scaling function in the near-adiabatic limit, which decays exponentially with x. For a general quantum critical point, the momentum resolved defect density is nðp; τÞ ¼f LZ p z τ zν=ðzνþ1Þ ; ð5Þ and the defect density follows the Kibble-Zurek scaling 1,2 as with d, z, and ν being the spatial dimension, dynamical critical exponent, and the exponent of the correlation length, respectively.
Gapless quench. This is realized for The eigenvalues of the Hamiltonian are always ±|p|, irrespective of the value of Δ 0 , which makes this parameter ramp exactly solvable by plugging these parameters into the non-Hermitian Schrödinger equation (see Methods). The dynamics is nevertheless non-trivial and as Δ and Γ evolve with time, defects are produced in spite of the fact that the instantaneous eigenvalues do not change. Starting from the ground state at t = 0, for p < 0 the wavefunction only picks up a phase factor as expðÀiptÞ½1; 1 T = ffiffi ffi 2 p . On the other hand, the p > 0 ground state at t = 0 evolves to where the second term is generated by the time dependence. For τ → ∞, this expression agrees with the right eigenfunction of the final Hamiltonian (up to normalization). Not only does the instantaneous eigenvalue remain unchanged, i.e., E ± (p) = ±|p|, but also the time evolution is only sensitive to the instantaneous eigenenergies, namely the wavefunction contains exp(±ipt) exponential factors only. This is then used in Eq. (15) to yield 〈σ y (τ)〉 + Δ 0 ln(2W/Δ 0 )/2π~τ −2 and 〈σ z (τ)〉~τ −1 . In this case, the τ → ∞ solution coincides with the instantaneous expectation value after the time evolution and an adiabatic theorem seems to hold. Since the instantaneous spectrum remains unchanged and gapless throughout, the above scaling cannot originate from the usual argumentation of Kibble-Zurek scaling. This is analogous to quenching along a gapless line 27 within the Hermitian realm. In order to appreciate the role of wavefunction normalization in Eq. (15), we have also evaluated it without the denominator: 〈σ z (τ)〉 approaches a constant, while 〈σ y (τ)〉~ln(τ), without any well-defined adiabatic limit for τ → ∞.
PT-symmetric ramp. Now we consider a fully PT-symmetric Kibble-Zurek problem, when the instantaneous spectrum is always real. We choose Δ = Δ 0 , Γ = Δ 0 t/τ such that the time evolution ends exactly at an EP. The spin expectation values are 〈σ y (τ)〉 + Δ 0 ln(W/Δ 0 )/π~τ −2/3 , and as shown in Fig. 1 from the numerics. The wavefunction for t = τ → ∞ agrees with the non-normalized right eigenfunction of the final non-Hermitian Hamiltonian, similarly to the gapless quench.
Since the spectrum ±|p| is linearly gapless at the critical point, this defines z = 1, leaving us with ν = 1/2 for the exponent of the correlation length. Therefore, the Kibble-Zurek scaling of the defect density in one dimension predicts~τ −dν/(zν+1) = τ −1/3 scaling. However, this exponent is different from Eq. (8). We demonstrate that the correct exponent is indeed −2/3 and present a generalized Kibble-Zurek scaling to account for that.
First of all, the numerical data indicate that the momentum resolved defect density, σ z (p, τ), follows the scaling relation with f PT (x) the universal scaling function shown in Fig. 2. Upon integrating this with respect to p by changing variable x = p(τΔ 0 ) 1/3 /Δ 0 , the τ −2/3 scaling of the defect density follows. In Eq. (9), the τ exponent 1/3 originates from the zν/(zν + 1) combination of critical exponents and the p stems from the z = 1 dynamical critical exponents. Therefore, this expression is generalized to an arbitrary critical point for the momentum resolved defect density as which, after performing a d-dimensional momentum integral, gives n~τ −(d+z)ν/(zν+1) . Notice the τ-dependent prefactor of the scaling functions in Eqs. (9) and (10) compared to the Hermitian case Eq. (4). We next provide three complementary explanations for this modified scaling. In an a/diabatic picture, excitations are created by populating the excited state similarly to Hermitian dynamics, but only its component perpendicular to the ground state represents defect production. As we approach the EP with increasing time, we enter into the diabatic regime at the transition timet tr , where adiabatic time evolution breaks down, the dynamics gets frozen and defect production kicks in. The component of the excited state perpendicular to the ground state at this instance has an amplitude sin(θ p ) as the ground and excited states are not orthogonal in general 29 . For Eq. (2), this is evaluated for small momentum states close to the EP, which are the most sensitive to diabatic time evolution, as sinðθ p%0 Þ ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi Δ 2 À Γ 2 p =Δ $ 1=t tr at the adiabatic-diabatic transition: namely the angle becomes proportional to the energy gap. This results in a τ −zν/(zν + 1) suppression factor for the defect density. For the Hermitian case, orthogonality ensures that sin(θ) = 1.
The numerically determined defect density. We plot the defect density from the normalized (blue circles) and unnormalized wavefunction (red squares) for the PT-symmetric ramp, as well as for the full nonhermitian drives (green triangles), measured from its adiabatic value. The black dashed lines denote the τ −1/3 , τ −2/3 , and τ −1 scaling. The cutoff |p| < W = 10Δ 0 and 10Γ 0 , respectively, does not alter the dynamics, with other values yielding similar scaling In a more dynamical picture, defects are created directly in the state perpendicular to the ground state, which decays to the ground state with a rate 1=t tr reducing the Hermitian Kibble-Zurek scaling by the probability to remain in the perpendicular state, 1=t tr . At an EP, there is only a single eigenstate and any perpendicular component decays 30 . Close to an EP, the state perpendicular to the ground state initially decays towards the ground state, which is followed by revival and periodic oscillation with frequency E + (p). An unnormalized state perpendicular to the ground state evolves as |⊥(t)〉 = 2icot(θ p )sin (E + (p)t)|GS〉 + exp(iE + (p)t)|⊥〉 due to a time independent non-Hermitian Hamiltonian, where θ p is the angle between the ground state and the excited state, 〈GS|⊥〉 = 0 and θ p vanishes at EP. For t~π/2E + (p), it develops a large component parallel to the ground state with length cot(θ p ). However, when the driving rate, ∂ t E + (p)/E + (p) is larger than the revival frequency, the system does not have enough time for revival and only the initial decay is probed. This gives the very same condition as the a/diabatic transition and the decay time ist tr .
Mathematically, the prefactor in Eq. (10) arises because the norm of the wavefunction also changes due to non-unitary time evolution. This follows from 31 d〈Ψ p (t)|Ψ p (t)〉/dt = 2Γ(t)〈Ψ p (t)|σ z | Ψ p (t)〉. Due to the slow time evolution, states with large momentum p have large energy, and are hardly affected by the time-dependent term and the corresponding wavefunction norm hardly changes. The low energy and small momentum states are the most influenced by the non-Hermitian and non-adiabatic time evolution. At short times t ( τ ð Þ, both the small matrix element and the Γ(t) prefactor block the growth of the wavefunction norm, but at a distancet tr from the critical point, adiabaticity breaks down. Afterwards, diabatic time evolution takes place, and the norm of the wavefunction gets enhanced bŷ t tr $ τ 1=3 , thus suppressing the defect density.
Full non-Hermitian drive. A full non-Hermitian drive is realized for Δ = 0, Γ = Γ 0 t/τ, which represents the non-Hermitian Kibble-Zurek problem and is equivalent to quenching the imaginary tachyon mass 17 . The instantaneous spectrum contains EPs located at |p| = Γ. By expanding around the EP, the spectrum scales as E ± ðp≳ΓÞ ¼ ± ffiffiffiffiffi 2Γ p ffiffiffiffiffiffiffiffiffiffiffi p À Γ p in the PT-symmetric regime, and as E ± ðp≲ΓÞ ¼ ± i ffiffiffiffiffi 2Γ p ffiffiffiffiffiffiffiffiffiffiffi Γ À p p in the broken PT-symmetry sector. Altogether the dynamical critical exponent is thus z = 1/2, while ν = 1. During the time evolution, all instantaneous eigenvalues are imaginary for Γ(t) > |p| and real for Γ(t) < |p|, separated by an EP. This critical point, which is located at p = Γ 0 t/τ, moves in momentum space during the time evolution at a speed Γ 0 /τ, producing defects in the process. This results in 〈σ y (τ)〉 = 0 and This is depicted in Fig. 1 from the numerical solution of Eq. (14) (see Methods). Here, it is again crucial to properly normalize the wavefunction as in Eq. (15). Without the normalization, the spin expectation value changes exponentially in time due to the imaginary energy eigenvalues.
The numerically obtained adiabatic value for 〈σ z (τ → ∞)〉 is also corroborated from diagonalizing the non-Hermitian Hamiltonian analytically, and using its normalized right eigenfunction: By taking a closer look at near-adiabatic dynamics, we calculate the momentum resolved spin expectation value, σ z (p,τ) numerically. This is illustrated, together with its critical scaling, in Fig. 3.
There is a clear difference in the contribution of states with imaginary as opposed to real instantaneous eigenvalues to this expectation value. This also differs significantly from the Landau-Zener transition probability of the corresponding Hermitian system 3 . From the numerical data, the defect density obeys the scaling function with f nh (x) the corresponding universal scaling function, depicted in Fig. 3 and σ eq z ðpÞ ¼ < . Upon integrating this with respect to p, the 1/τ scaling of the defect density follows.
This allows us to conjecture that for a general EP, the momentum resolved defect density satisfies the same form as Eq. (10), which agrees with Eq. (13) with z = 1/2. Therefore, the induced defect density vanishes as τ −(d + z)ν/(zν + 1) upon traversing the EP adiabatically, similarly to the PTsymmetric ramp.

Discussion
Throughout this work, a non-Hermitian Hamiltonian was used to generate the time evolution. From a more direct physical perspective, such dynamics follow from a Lindblad master equation in combination with a continuous measurement 16,32,33 . More specifically, consider a system, described just by the Hermitian part of our Hamiltonian, coupled to an environment inducing radiative decay in the individual two-level systems for each momentum p with a rate Γ. In terms of a Lindblad equation, this results in quantum jump operators equal to σ − describing an incoherent decay upon emitting a photon. Equivalently, one can map the Lindblad dynamics for the system's density matrix onto a quantum jump trajectory picture, where pure states are evolved upon averaging over trajectories according to the following prescription. A single trajectory is specified by a set of quantum jumps at times, that can be sampled from a given probability distribution. At those times the quantum jump operator of the related Lindblad master equation is applied onto the quantum state. In our case this is the σ − operator and a quantum jump event can therefore be identified with the emission of a photon. In between the quantum jumps the dynamics is solely given by a non-Hermitian Hamiltonian, which in our case is the one in Eq.  Fig. 3 Momentum resolved defect density for the non-hermitian drive. The scaling of the numerically determined momentum resolved defect density, f nh (x) in the near-adiabatic limit is shown for several values of τ for the full non-hermitian drive around the equilibrium EP (2). In order to only select those trajectories without any photon emission event, which is the non-Hermitian evolution targeted in this work, we can further continuously monitor the system and measure the number of emitted photons. In case we only consider those realizations, where no emission event has taken place, we end up with an evolution precisely captured by our non-Hermitian Hamiltonian. Note that the measurement modifies the state of our system, since the absence of emission events effectively gradually forces the system over time towards the ground state of the two-level system as otherwise an emission event becomes too likely 32 . This continuous measurement further ensures that the wavefunction is always properly normalized. While the anti-Hermitian contribution forces the system towards the ground state of the two-level system, the Hermitian part counteracts this tendency by coherently transferring population back into the excited state. At an EP, these processes compete most strongly and generate a non-trivial attractor of the dynamics: this is the single eigenstate of the non-Hermitian system. This provides the additional channel for defect annihilation that we identified in our analysis. Far away from the EP, the Hermitian part of the Hamiltonian dominates for our setups, so that the dynamics is almost fully coherent. Finally, we note that a distinct non-Hermitian Kibble-Zurek scaling describing a different physical realization with different definition of the expectation value and the scalar product was studied in ref. 34 .
Single particle Hamiltonians of the form of Eq. (2) have already been realized in non-conservative classical and quantum systems 6,10,11,35 . The time-dependent control of the non-Hermitian term together with measuring the spin components are possible. By creating several copies of this two-level system, corresponding to distinct p's, the effective many-body dynamics and the non-Hermitian Kibble-Zurek scaling could in principle be detected.
To sum up, the universal features of non-Hermitian dynamics across EPs were investigated. We find that the adiabatic time evolution drives the initial wavefunction to a right eigenstate of the final non-Hermitian Hamiltonian, up to normalization, indicating that an adiabatic theorem probably exists for the systems under consideration. For a near-adiabatic crossing of an EP, defects are produced at a reduced rate, whose density obeys a generalized Kibble-Zurek scaling as τ −(d + z)ν/(zν + 1) . For the future it remains an open question how our results extend also to higher-order exceptional points.

Methods
Non-Hermitian time evolution. For the purpose of this work we consider timedependent parameters Δ(t) and/or Γ(t) yielding a Hamiltonian H(t). Initially, before we start our parameter ramps, we choose the Hamiltonian always to be Hermitian, i.e., Γ = 0, so that the initial condition as the ground state of the Hamiltonian is well-defined. At time t = 0 we start our time-dependent protocol over a time span τ. The time evolution follows from with |Ψ(t)〉 = ⊗ p |Ψ p (t)〉 for a given mode p. In general, the norm of the wavefunction is not conserved when time evolution is driven by a non-Hermitian Hamiltonian, so that 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 16,31,36 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 will quantify the defect production via with α = y, z, and N denoting the number of considered momentum states. In the absence of balanced gain and loss, one can shift the diagonal term in the Hamiltonian by any complex value, which does not affect the results. The reason is that such a shift leaves the expectation values in Eq. (16) invariant, since a simple (timedependent) change of the norm of the wavefunction is cancelled by explicitly using normalized expectation values.
Calculation of the defect density. The total defect density for the four considered settings are evaluated from the momentum resolved defect density after momentum integration. For example, by taking the PT-symmetric ramp, Eq. (9) defines the momentum resolved defect density. The total defect density is σ z ðτÞ ¼ Z dp 2π σ z ðp; τÞ: By changing variable x = p(τΔ 0 ) 1/3 /Δ 0 , this becomes σ z ðτÞ ¼ 1 Data availability The data that support the plots within this paper and other findings of this study are available from the corresponding author on request.