Non-Hermitian topological phases and dynamical quantum phase transitions: A generic connection

The dynamical and topological properties of non-Hermitian systems have attracted great attention in recent years. In this work, we establish an intrinsic connection between two classes of intriguing phenomena -- topological phases and dynamical quantum phase transitions (DQPTs) -- in non-Hermitian systems. Focusing on one-dimensional models with chiral symmetry, we find DQPTs following the quench from a trivial to a non-Hermitian topological phase. Moreover, the number of critical momenta and critical time periods of the DQPTs are found to be directly related to the topological invariants of the non-Hermitian system. We further demonstrate our theory in three prototypical non-Hermitian lattice models, the lossy Kitaev chain (LKC), the LKC with next-nearest-neighbor hoppings, and the nonreciprocal Su-Schrieffer-Heeger model. Finally, we present a proposal to experimentally verify the found connection by a nitrogen-vacancy center in diamond.

To date, non-Hermitian topological phases (NHTPs) have been classified and characterized according to their protecting symmetries [13][14][15]. Finding the dynamical signatures of these nonequilibrium topological matter has become an urgent topic for further theoretical and experimental explorations. In the literature, several dynamical probes to the topological invariants of one-and two-dimensional non-Hermitian phases have been proposed, such as the non-Hermitian extension of mean chiral displacements [40][41][42] and dynamical winding numbers [43][44][45][46]. In the meantime, DQPTs (i.e., nonanalytic behaviors of certain observables in time domain [48][49][50][51]) following a quench across the EPs of a non-Hermitian lattice model is investigated in Ref. [52], and the monotonic growth of a dynamical topological order parameter in time is observed if an isolated EP is crossed during the quench [52]. This discovery indicates an underlying relationship between the two notably different nonequilibrium phenomena, NHTPs and DQPTs. However, the more general connection between NHTPs and DQPTs, together with its possible experimental observations have not been revealed.
In this work, we uncover an intrinsic connection between the topological phases and DQPTs in one-dimensional (1D) non-Hermitian systems. In Sec. II, we develop our theoretical framework leading to the establishment of this connection. In Sec. III, we demonstrate the found connection in three different non-Hermitian lattice models, the lossy Kitaev chain (LKC), the LKC with next-nearest-neighbor (NNN) hoppings and pairings, and the nonreciprocal Su-Schrieffer-Heeger (NRSSH) model. In each model, a direct link between the bulk topological invariant of a non-Hermitian phase and the number of critical time and momenta of DQPTs following a quench to the corresponding phase is found. In Sec. IV, we discuss an experimental setup, the nitrogen-vacancy (NV) center in diamond, in which the discovered connection may be tested. We conclude this work and discuss potential future directions in Sec. V.

II. THEORY
In this section, we introduce a generic class of non-Hermitian lattice models and describing the topological characterization of its bulk states in Subsec. II A. In Subsec. II B, we introduce relevant quantities to characterize DQPTs in 1D non-Hermitian systems, and establish their connections with the underlying topological properties of the system.

A. NHTPs
We start with a non-Hermitian Hamiltonian H = H † , which describes particles in a 1D lattice subjecting to gains, losses and/or nonreciprocal effects. Under the periodic boundary condition, we can express the Hamiltonian of the system as H = k Ψ † k H(k)Ψ k , where k ∈ [−π, π) is the quasimomentum, Ψ † k (Ψ k ) is the two-component creation (annihilation) operator in momentum representation, and the Bloch Hamiltonian H(k) takes the general form Here h a,b (k) and g a,b (k) are real-valued functions of the quasimomentum k, i denotes the imaginary unit, and σ a,b are any two of the three Pauli matrices σ x , σ y , and σ z , with {σ a , σ b } = 0 for a = b. We will also denote the 2 × 2 identity matrix as σ 0 .
The non-Hermiticity of H implies that H(k) = H † (k) at a generic quasimomentum k.
The dispersion relation of non-Hermitian Bloch Hamiltonian H(k) is given by It is clear that ±E(k) are in general complex numbers. The spectrum of H(k) becomes gapless at zero energy if E(k) = 0. According to Eq. (2), this is achieved when both the following conditions are satisfied By solving these equations, we could obtain the quasimomentum k 0 at which the spectrum gap closes, and find the boundaries separating different gapped phases in the parameter space, which could also be the boundaries among different bulk topological phases of the system.
To characterize the topological properties of the gapped phases of H(k) (i.e., E(k) = 0 for all k), the usual recipe is to identify the symmetries of the system and construct the relevant topological invariants. From the commutation relation of Pauli matrices [σ a , in the sense that S 2 = σ 0 and SH(k)S = −H(k). The bulk topological phases of a non-Hermitian Bloch Hamiltonian with sublattice symmetry S can usually be characterized by a winding number w, defined as which describes the accumulated change of winding angle φ(k) throughout the first Brillouin zone (BZ). Note that the value of w is always real even though φ(k) is in general complex, as the imaginary part of φ(k) has no winding in the first BZ (see Ref. [43] for a proof). Furthermore, the w as defined in Eq. (5) [40] or the dynamical winding numbers of time-averaged spin textures [43].

B. DQPTs and their relations to NHTPs
DQPTs are characterized by nonanalytic behaviors of system observables as functions of time. They are usually found in the dynamics following a quench across the equilibrium phase transition point of a quantum many-body system (see Ref. [48][49][50][51] for reviews). The central object in the description of DQPTs is the return amplitude where |Ψ is the initial many-particle state (usually taken as the equilibrium ground state of the system before the quench) and U (t) is the evolution operator of the system following a quantum quench (or some other nonequilibrium protocols). Formally, G(t) mimics the dynamical partition function of the post-quench evolution. When G(t) = 0 at a critical time t c , the initial state evolves into its orthogonal state. The rate function of return probability is the number of degrees of freedom of the system) or its time derivatives would then become nonanalytic at t = t c , signifying a DQPT. Accompanying theoretical discoveries [53][54][55], DQPTs have been observed in cold atoms [56][57][58][59], trapped ions [60,61], superconducting qubits [62], nanomechanical and photonic systems [63][64][65].
To relate DQPTs with topological phases in non-Hermitian systems, we focus on a unique class of quench protocol, in which the system is initialized with equal populations but no coherence on the two bands of the non-Hermitian Bloch Hamiltonian H(k) in Eq. (1), i.e., an infinite-temperature initial state k∈BZ ρ 0 , with ρ 0 = σ 0 /2 being the single-particle density matrix. The evolution of ρ 0 at time t > 0 is governed by H(k), and the return amplitude G(k, t), defined as the expectation value of evolution operator U (k, t) = e −iH(k)t over the initial state ρ 0 reads where Eqs. (1) and (2) have been used to reach the second equality. When possible DQPTs happen, we would have cos[E(k)t] = 0, leading to the critical times This seemingly innocent expression yields rather different predictions for Hermitian and non-Hermitian systems. In a Hermitian system, where the dispersion relation E(k) is always real and positive, we would have a set of critical times t n (k) for each quasimomentum k. However, the resulting non-analyticity in the rate function g(t) is simply originated from the oscillatory dynamics of a single Bloch state rather then an actual phase transition, which only happens in thermodynamic limit (N → ∞). On the other hand, when H(k) is non-Hermitian, we have E(k) ∈ C in general, and real critical times t n (k) emerge only at the critical momenta k c where E(k c ) ∈ R. According to Eq. (2), this is equivalent to the fulfillment of the following two conditions: which only yield solutions at isolated values of k. For a critical momentum k c satisfying both the Eqs. (8) and (9), we would have G(k c , t n ) = 0 for all n ∈ Z. In the thermodynamic limit, the rate function of return probability for the many-particle initial state k∈BZ ρ 0 is given by which will have discontinuous first-order time derivatives at all t n (k c ). Note that by taking the limit N → ∞, the distribution of t n (k) on the complex time plane changes from isolated points to a continuous line, whose crossings along the real-time axis correspond to the critical time of genuine DQPTs in the sense of Fisher zeros [48].
The connection between DQPTs and NHTPs in chiral-symmetric 1D systems becomes transparent at this stage. First, we note that the Eq. Putting together, we have uncovered an intrinsic relation between the topological phases and DQPTs in non-Hermitian systems, which not only bridges the gap between these two diverse fields, but also provides a way to probe the NHTPs through nonequilibrium dynamics.
To make the connection more explicit, we will study the DQPTs in a couple of prototypical 1D non-Hermitian lattice models in the following section. Besides the rate function g(t), we will also investigate the real-valued, noncyclic geometric phase of the return amplitude G(k, t) [52], which is defined as where the total phase and the dynamical phase (see Appendix A for more details about these phase factors). The noncyclic geometric phase has been shown to contain important information about DQPTs in both Hermitian [72][73][74][75][76] and non-Hermitian [52] systems. At a given time, the winding number of the geometric phase in the first BZ can be further employed to construct a dynamical topological order parameter (DTOP), which is defined as It takes a quantized jump whenever the evolution of the system passes through a critical time of the DQPT. Note that the range of integration over k depends on the symmetry has the inversion symmetry with respect to k = 0, i.e., Φ G (k, t) = Φ G (−k, t), we can perform the integral over a reduced BZ with k ∈ [0, π] in the evaluation of ν in Eq. (14).

III. MODELS AND RESULTS
In this section, we demonstrate the connection between NHTPs and DQPTs in three typical non-Hermitian 1D lattice models. In each subsection, we introduce the model that will be investigated first and establish its bulk topological phase diagram. After that, we consider the DQPTs in the model following the quench from a trivial phase to different non-Hermitian phases (either trivial or topological), and unveil the relationship between the critical times and momenta of the DQPTs and the topological invariants of the post-quench non-Hermitian system. In the lossy Kitaev chain and its next-nearest-neighbor extension, we observe a one-to-one correspondence between the NHTPs and DQPTs. In the nonreciprocal SSH model, we find that while a topologically nontrivial post-quench system always imply DQPTs following the quench, the reverse may not be true in general, and possible reasons behind such an anomaly will be discussed.

A. The lossy Kitaev chain
We first consider a non-Hermitian variant of the Kitaev chain, which describes a 1D topological superconductor with onsite particle loss. In momentum representation, the Hamiltonian of the model takes the form is the Nambu spinor operator and c † k is the creation operator of an electron with quasimomentum k. The Bloch Hamiltonian H(k) in Nambu basis is given by where Here the real parameters J, ∆ and u denote the nearest-neighbor hopping amplitude, superconducting pairing amplitude and chemical potential. v ∈ R characterizes the strength of onsite particle loss. Following the discussions of Subsec. II A, we see that H(k) possesses the sublattice symmetry S = σ x , i.e., SH(k)S = −H(k). Furthermore, it also has the generalized particle-hole symmetry C = σ x and time-reversal symmetry T = σ 0 , in the sense where performs matrix transpo-sition. H(k) thus belongs to an extension of the symmetry class BDI in the periodic table of non-Hermitian topological phases [14]. In the meantime, H(k) possesses the inversion symmetry P = σ z as PH(k)P −1 = H(−k), which guarantees the correspondence between its bulk topological invariant w (as defined in Eq. (5)) and the number of Majorana edge modes under the open boundary condition [14].
According to Subsec. II A, the complex energy spectrum of LKC takes the form which will become gapless when Combining these equations, we find the gapless quasimomenta for |u| < |J|, and the boundary between different NHTPs as Geometrically, the trajectory of vector forms an ellipse on the h y -h z plane, which is centered at (0, u). When the gapless condition Eq. (21) is satisfied, the spectrum E ± (k) hold a pair of EPs at (±v, 0) on the h y -h z plane, which are passed through by the vector h(k). Whether the two EPs are encircled or not by the trajectory of h(k) when k scans through the first BZ then distinguishes two possible NHTPs. With Eqs. (17) and (21), it is not hard to show that when u 2 /J 2 + v 2 /∆ 2 < 1 (> 1), the two EPs are encircled (not encircled) by the trajectory of h(k). The topological invariant that distinguish these two phases has the form of Eq. (5), where the winding angle φ(k) for the LKC is explicitly given by To link the NHTPs of LKC with the DQPTs, we employ the protocol introduced in Subsec. II B, with the initial state ρ 0 = σ 0 /2 and the dynamics being governed by the Hamiltonian H(k) in Eq. (15). According to Eq. (6), the return amplitude at a later time Eq. (17). From Eqs. (9) and (7), we find the critical momenta and times to be It is clear that when u 2 /J 2 + v 2 /∆ 2 < 1, there are real solutions of t n (k c ) for all n ∈ Z, and the two critical momenta ±k c are coincide with the gapless quasimomenta ±k 0 , yielding On the other hand, there is no critical momenta and t n is always imaginary when u 2 /J 2 + v 2 /∆ 2 > 1, yielding no DQPTs at any real time t. When u 2 /J 2 + v 2 /∆ 2 = 1, which corresponds to a gapless post-quench phase, we will have t n (k c ) → ∞ for any solutions of critical momenta ±k c , and the resulting DQPTs are not observable. For completeness, we numerically compute the return rates and geometric phases with the help of Eqs. (10) and (11) for the cases with and without DQPTs in Fig. 2(a,c) and 2(b,d), respectively. As expected, DQPTs are only observed when the system is quenched to a nontrivial NHTP with the winding number w = 1.
Combining the analysis in this subsection, we obtain a one-to-one correspondence between the NHTPs and DQPTs in the LKC, which is summarized in Table I. This connection not only unifies the NHTPs and DQPTs in the system, but also provides a way to dynamically distinguishing the different NHTPs of LKC and detecting the gapless quasimomenta, as exemplified by Figs. 2(a,c).

B. The lossy Kitaev chain with next-nearest-neighbor hoppings and pairings
We next consider the LKC with NNN hoppings and pairings, which could possesses NHTPs with larger topological invariants. In the momentum space and Nambu spinor basis, the NNN LKC is described by the Hamiltonian

Condition Geometric
Winding Critical time picture number and momenta No EPs are w = 0 No k c and t n encircled by h(k) No DQPTs Here u is the real part of chemical potential, (J 1 , ∆ 1 ) and (J 2 , ∆ 2 ) are the nearest-neighbor and next-nearest-neighbor hopping and pairing amplitudes, respectively. It is not hard to verify that the H(k) here belongs to the same symmetry class as the LKC, with the same set of time-reversal, particle-hole, sublattice and inversion symmetries. The dispersion relations E ± (k) of H(k) share the same form with Eq. (17), yielding the gapless conditions u + J 1 cos k + J 2 cos 2k = 0.
By solving Eq. (27), we could obtain at most four possible gapless quasimomenta ±k ± 0 as According to Eq. (26), the boundary between different NHTPs is then determined by The explicit expression of the phase boundary in terms of system parameters is tedious, and will be left for numerical calculations. Geometrically, the trajectory of real vector Winding Critical times picture number and momenta No k c and t n by h(k) No DQPTs  (7) as where n ∈ Z. In parallel with the discussions of Subsec. III A, we could summarize the relationship between NHTPs and DQPTs in the NNN LKC model by Table II. Again, we obtain a one-to-one correspondence between these two nonequilibrium phenomena, which also provides us with a way to detect NHTPs with large winding numbers and to locate the phase boundaries between them.  Fig. 4(a), with the two pairs of critical momenta ±k ± c given by Eq. (30) and imaged by the 2π-jumps of geometric phase Φ G (k, t) in Fig. 4(d). In Figs. 4(b,e), the post-quench system is in a NHTP with w = 1, and DQPTs are repeated at only one critical period T (k + c ) = π/ h 2 y (k + c ) − v 2 of g(t) in Fig. 4(b), with 2π-jumps of geometric phase Φ G (k, t) observed at the critical momenta ±k + c in Fig. 4(e). In Figs. 4(c,f), the post-quench system is in a trivial phase with w = 0, and no signatures of DQPTs are observed in the rate function g(t) and geometric phase Φ G (k, t). Putting together, our numerical results confirm the connection between the NHTPs and DQPTs of the NNN LKC model, as summarized in Table II. Furthermore, the results presented here , t 6 (k + c ) from left to right, whose explicit values are obtained from Eqs. (30) and (31). In panels (b,e), the ticks along the horizontal axis are the critical times t n (k + c ) for n = 1, ..., 5 from left to right. An extra amount of 2π-jump in the geometric phase Φ G (k, t) is observed at the corresponding critical momentum when a DQPT happens.
should be directly extendable to non-Hermitian models in the same symmetry class as the LKC, but with even longer-range hopping and pairing amplitudes.
is the creation operator on the two sublattices a and b of the SSH model, and k ∈ [−π, π) is the quasimomentum. The Bloch Hamiltonian H(k) is explicitly given by with Here J 1 ± γ and J 2 are the intracell and intercell hopping amplitudes. A finite γ makes the intracell hopping asymmetric, leading to a non-Hermitian H(k). From now on, we assume The bulk spectrum of H(k) takes the form With Eq. (33), we see that the dispersion is gapless at zero energy when the following two conditions are met which directly yield the phase boundary curves and the gapless quasimomenta k 0 = 0, π. The connection between DQPTs and NHTPs in the NRSSH model can be built as follows.
Choosing the initial state to be ρ 0 = σ 0 /2 as in Subsec. II B, the dynamics of the system at Winding Critical times picture number and momenta No k c and t n No EPs are No DQPTs DQPTs at t 0,π n ∀n ∈ Z, k c = 0, π  (7), i.e., where n ∈ Z. Combining these equations with the gapless conditions in Eqs. (35) and (36), we could immediately identify the relationship between NHTPs and DQPTs in the NRSSH model, as listed in Table III.
From the table, we observe that since there is only a single critical momentum for the NHTPs with w = 1/2, there is also a unique set of critical times (t 0 n or t π n for n ∈ Z) for the DQPTs in this case. This is in contrast with the NHTPs having w = 1, for which both k c = k 0 = 0 and π are the critical momenta, and DQPTs at two different critical time periods T (k c = 0, π) are expected in the post-quench dynamics. In the meantime, we also observe an anomalous case as shown in the last row of Table III. In this case, DQPTs are found when the post-quench system is in a trivial phase with w = 0. Therefore, even though a nontrivial topological phase of the NRSSH model always lead to a unique set of DQPTs following the quench to that phase, the reverse is not true in general. Such a breakdown of the one-to-one correspondence between the DQPTs and NHTPs in the NRSSH model might be due to the absence of inversion symmetry, as compared with the situations in the LKC and its NNN extension. Nevertheless, the most intriguing phase of the NRSSH model, i.e., the one with w = 1/2 can still be distinguished from the other phases through the DQPTs.
Therefore, the connection between NHTPs and DQPTs we discovered can still be used as a powerful tool to probe the details of the NHTPs in the NRSSH model.
For completeness, we present the DQPTs in three typical post-quench phases of the NRSSH model in Fig. 6. In Figs. 6(a,d), the post-quench phase has winding number w = 1, and DQPTs are observed as cusps in the rate function g(t) at two sets of critical times t 0,π n . At each critical time, a 2π-jump in the geometric phase Φ G (k, t) is observed at both the critical momenta k c = 0, π. In Figs. 6(b,e), the post-quench phase has winding number w = 1/2, and DQPTs are found at a unique set of critical time t 0 n , where a 2π-jump in the geometric phase Φ G (k, t) is observed around k c = 0. In Figs. 6(c,f), the post-quench phase is trivial and no DQPTs are found in the post-quench dynamics. Putting together, we found that the NHTPs and DQPTs in the NRSSH model are also two closely related phenomena, and the later can be employed to dynamically probe the properties of the former.

IV. EXPERIMENTAL PROPOSAL
With all the theoretical and numerical results presented above, we now sketch an experimental proposal in which our predicted connection between NHTPs and DQPTs may be verified. Recently, a setup containing an NV center in diamond has been employed to realize the PT-symmetry breaking transition of a non-Hermitian two-level Hamiltonian [29]. The general idea is to dilate a PT-symmetric Hamiltonian into a Hermitian one, and execute the dynamics with the dilated Hamiltonian. Since all the three models discussed in the previous section possess two bulk bands together with the PT-symmetry, the setup proposed in Ref. [29] tends out to be an ideal platform in which the topological invariants and DQPTs of our systems can be detected.
The Hamiltonian we are interested in, as shown in Eq. (1) can be generally expressed as The strategy of Ref. [29] is to dilate H(k) into a Hermitian counterpart with the help of an ancilla qubit.
The dilated Hamiltonian H (k, t) yields the Schrödinger equation where |Ω(k, t) denotes the state of the composite system. With an appropriate postselection scheme, the measurement results can be restricted to a unique outcome for the ancilla qubit [29]. Within the scheme, the composite state |Ω(k, t) takes the form where ω(t) is an appropriate linear operator, and the ancilla qubit basis |± are chosen to be the eigenstates of Pauli matrix σ y In the experiment, a −π/2 pulse is applied following the evolution, and only the measurement results inside the state manifold |Ψ(k, t) |− is post-selected.
The explicit form of dilated Hamiltonian H (k, t) is not unique. A convenient choice realized by the experiment in Ref. [29] is where with the time-dependent operator M (t) ≡ ω † (t)ω(t) + σ 0 . Expanding Λ(k, t) and Γ(k, t) by the Pauli matrices σ x,y,z and σ 0 , we can further express H (k, t) as where the real coefficients A i (k, t) and B i (k, t) for i = 0, 1, 2, 3 can be obtained numerically [29].  [47], as suggested in Ref. [43]. In a very recent experiment, the non-Hermitian topological phases of a nonreciprocal SSH model have been detected in an NV center setup following the universal dilation scheme [30], which confirms the applicability of the experimental proposal.

V. SUMMARY
In this manuscript, we establish a relationship between NHTPs and DQPTs in 1D systems. DQPTs are found when the system is quenched from a trivial to a non-Hermitian In this work, our theory is applied to one-dimensional two-band models with chiral symmetry. Our initial attempts also suggest that the theoretical framework presented here is generalizable to chiral-symmetric multiple-band models [41,78]. However, due to the complexity of multiple-band systems in the study of NHTPs and DQPTs, we expect that our theory would subject to appropriate modifications when it is applied to these systems. This where a, b = x, y, z and a = b. It can be equivalently written as where and σ = (σ a , σ b ). It is clear that n(k) is a unit vector with n(k) · n(k) = 1.
According to Eqs. (6) and (12), the total phase of the return amplitude reads It is clear that Φ(−k, t) = Φ(k, t) once E(−k) = ±E(k). This is clearly the case for the models considered in Subsecs. III A and III B according to the expressions of their bulk spectrum E ± (k). Instead, for the model studied in Subsec. III C, the total phase does not have the parity (inversion) symmetry.
To obtain the dynamical phase, we introduce the biorthogonal representation of non-Hermitian systems. In this representation, the right and left eigenvectors {|ψ s (k) |s = ±} and {|ψ s (k) |s = ±} of H(k) satisfy the eigenvalue equations and The Hamiltonian H(k) can also be expressed in this representation as The time-evolution operators in the spaces of right and left eigenvectors are respectively.
According to the definition of dynamical phase Φ D (k, t) in Eq. (13) of the main text, we With Eqs. (A9) and (A10), we can recast Φ D (k, t) into a more explicit form. The denominator of the integrand tends out to be Furthermore, the numerator in the integrand of Eq. (A11) yields Putting together, we find the dynamical phase to be Referring to the main text, we see that as E(−k) = E(k) for the LKC and NNN LKC models, we also have Φ D (−k, t) = Φ D (k, t) for these two models. Therefore, we conclude that the geometric phases Φ G (k, t) = Φ(k, t) − Φ D (k, t) of the LKC and NNN LKC models in the main text both possess the inversion symmetry, i.e., Φ G (−k, t) = Φ G (k, t). This allows us to confine the range of integration to half of the first BZ, e.g., k ∈ [0, π] for the calculation of dynamical topological order parameters in Eq. (14) for these two models.
Comparatively, for the NRSSH model studied in Subsec. III C of the main text, we have , and the whole first BZ k ∈ [−π, π] should be employed in the calculation of its dynamical topological order parameter. Experimentally, information about the geometric phase may be directly obtained by measuring the complex spectrum dispersion E(k) of the system [29]. we also notice that ν(t) may not take quantized values between certain pairs of critical times (e.g., for t ∈ (t 1 , t 2 )), which might be due to our choice of reduced BZ k ∈ [0, π] in the calculation of ν(t). Nevertheless, the quantized jump of DTOP across each critical time already provides us with essential information about the drastic topological change of the system when undergoing a DQPT. jump of DTOP every time when the system evolves across a critical time, implying the existence of a DQPT. Furthermore, the values of DTOP remain quantized between any pair of the critical times, which is expected as the whole Brillouin zone k ∈ [−π, π] is used in the calculation of ν(t). Besides, we also notice that the value of DTOP changes monotically in time when the post-quench system has a half-quantized winding number w = 1/2, which is consistent with the connection bewteen the exceptional non-Hermitian topology and DQPTs as first observed in Ref. [52].
, t 6 (k + c ) from left to right, whose explicit values are obtained from Eq. (31). In panel (b), the ticks along the horizontal axis are the critical times t n (k + c ) for n = 1, ..., 5 from left to right.