Anomalous-order exceptional point and non-Markovian Purcell effect at threshold in one-dimensional continuum systems

For a system consisting of a quantum emitter coupled near threshold (band edge) to a one-dimensional continuum with a van Hove singularity in the density of states, we demonstrate general conditions such that a characteristic triple level convergence occurs directly on the threshold as the coupling $g$ is shut off. For small $g$ values the eigenvalue and norm of each of these states can be expanded in a Puiseux expansion in terms of powers of $g^{2/3}$, which suggests the influence of a third-order exceptional point. However, in the actual $g \rightarrow 0$ limit, only two discrete states in fact coalesce as the system can be reduced to a $2 \times 2$ Jordan block; the third state instead merges with the continuum. Moreover, the decay width of the resonance state involved in this convergence is significantly enhanced compared to the usual Fermi golden rule, which is consistent with the Purcell effect. However, non-Markovian dynamics due to the branch-point effect are also enhanced near the threshold. Applying a perturbative analysis in terms of the Puiseux expansion that takes into account the threshold influence, we show that the combination of these effects results in quantum emitter decay of the unusual form $1 - C t^{3/2}$ on the key timescale during which most of the decay occurs. We then present two conditions that must be satisfied at the threshold for the anomalous exceptional point to occur: the density of states must contain an inverse square-root divergence and the potential must be non-singular. We further show that when the energy of the quantum emitter is detuned from threshold, the anomalous exceptional point splits into three ordinary exceptional points, two of which appear in the complex-extended parameter space. These results provide deeper insight into a well-known problem in spontaneous decay at a photonic band edge.

In general, the exceptional point (EP) represents a defective point in the parameter space of a given Hamiltonian where diagonalization is no longer possible as two or more eigenstates coalesce. This situation is different than an ordinary degeneracy, at which the eigenvalues coincide while the corresponding eigenstates remain orthogonal to one another, and hence never coalesce. In the usual picture, we refer to an exceptional point with N coalescing eigenvalues as an EPN. Then the simplest representation of the Hamiltonian at an EPN would contain an N × N Jordan block 1,47,48 . Meanwhile, in the vicinity of the EPN the eigenvalues (and other experimentally measurable quantities) can be expanded in terms of a Puiseux series of the form E j ( ) =Ē + α j ( −¯ ) 1/N + . . . where¯ is the exceptional point andĒ is the coalesced eigenvalue [1][2][3]14,49 . While most studies have focused on the simplest case involving the coalescence of two eigenstates 5,50,51 , higher-order exceptional points involving three [52][53][54][55][56][57][58] or more 59,60,65 eigenstates have received more attention in the last few years [61][62][63][64] . Proposed applications involving exceptional points include enhanced sensing 57,58,66,67 as well as modified spontaneous emission [68][69][70] and dynamical control 46,63,[71][72][73][74][75][76] .
In this paper, we analyze an exceptional point that generically occurs at the continuum threshold (band edge) when a quantum emitter is coupled to a onedimensional (1-D) system under certain conditions. As a result of the van Hove singularity in the 1-D density of states with characteristic divergence 1/ √ E − E th , in which E th is the continuum threshold, there occurs a nearby level triplet consisting of a resonance, an antiresonance and a bound state. In the vicinity of the threshold, their eigenstates and eigenvalues can be expanded in a Puiseux series in terms of g 2/3 , in which g is the coupling to the quantum emitter. This seems to suggest an EP3 at which all three levels will coalesce in the limit g → 0. In fact, we find that the exceptional point is technically an EP2, despite behaving for physical purposes essentially like an EP3. Further, the decay width Γ of the resonance is enhanced near the threshold (Γ ∼ g 4/3 ) in comparison to the usual case in which arXiv:2104.06929v2 [quant-ph] 12 Jul 2021 Fermi's golden rule could be applied (Γ ∼ g 2 ) 77 . This enhancement can be viewed as an implication of the Purcell effect, which predicts that the decay width of a quantum emitter is enhanced when it is tuned to the natural frequency of a cavity or waveguide [78][79][80] .
After establishing these properties of the exceptional point, we analyze in detail its influence on the survival probability of the excited quantum emitter state. During the critical timescale in which most of the decay occurs, a power-law evolution manifests itself in which the exponent is determined by two factors: one coming from the anomalous exceptional point and the other coming from the continuum threshold (branch-point effect) 46,81 . The dynamics on this timescale can be viewed as a non-Markovian correction to the previously-mentioned Purcell effect. We emphasize that this effect should be rather universal; in particular, the same dynamics occurs in a well-known problem involving spontaneous emission near the edge of a photonic band gap [82][83][84] , but neither the power-law decay nor the exceptional point shaping those dynamics were previously noticed.
The models considered in this paper are written in terms of a microscopic Hamiltonian that includes a structured reservoir to describe the environmental influence on the quantum emitter, which means that there is a background continuum with a well-defined bandwidth and density of states 77,[82][83][84][85][86][87][88][89][90][91][92][93][94][95][96][97] . This is opposed to the coupled mode theory that has often been used to describe exceptional-point phenomena, in which one assumes that the essential physics can be described in terms of a few interacting resonance modes, while the microscopic details of the system are set aside (including the continuum). Although the microscopic description applied here generally requires more effort to analyze the influence of the exceptional point, in some situations the properties of the discrete spectrum and those of the continuum cannot so neatly be disentangled 46,48 .
In Sec. II below we introduce a simple system consisting of a quantum dot and a tight-binding chain in order to explicitly demonstrate the properties of the threshold EP. In Sec. II B we rewrite the model in terms of the quadratic eigenvalue problem, which enables us to describe the physical and mathematical properties of the exceptional point. Then in Sec. III we study the spectrum when the quantum emitter frequency is detuned from the threshold, observing several ordinary EP2s that occur in this vicinity. As the detuning and the coupling are both shut off, these ordinary EP2s merge on the threshold as the previously discussed anomalousorder exceptional point appears. We show the Jordan block structure and other technical properties of the anomalous-order exceptional point in Sec. IV [see also App. A]. Then in Sec. V we consider the influence of the EP on the time evolution of the excited quantum emitter state near the threshold. We demonstrate that a (1−Ct 3/2 )-type power law decay occurs as a result of the influence of the anomalous exceptional point. We discuss our results in Sec. VI and present our argument that two conditions should be satisfied for the anomalous-order EP to appear in a generic 1-D quantum system. We also discuss a potential experiment in circuit quantum electrodynamics (QED).

II. MODEL AND FORMALISM
The primary model that we consider in this paper is a quantum wire superlattice with an attached quantum dot 98 . The quantum wire is modeled by an infinite tightbinding array such that our Hamiltonian takes the form in which c † j is the creation operator at site j of the array, while J is the resonant coupling strength between array elements, which we set to unity as the unit of the energy. Meanwhile d † is the creation operator for the quantum dot excited state with excitation energy d . The dot is coupled to the 0 element of the chain c † 0 through the small coupling parameter g. Note that because there is no particle-particle interaction in this Hamiltonian, we can study the dynamics within the single-particle subspace without approximation; further, the following analysis could equally well apply for an equivalent Bosonic realization of the model such as in waveguide QED in Ref. 99 .

A. Eigenvalues: Siegert boundary condition, effective Hamiltonian and dispersion equation
To analyze the spectrum of our model we next write our physical solution with the outgoing (or Siegert 39 ) boundary conditions Writing the Schrödinger equation x|H|ψ = E x|ψ far from the dot |x| > 1 yields the expression where we have set J = 1. Plugging into this Eq. (2) we find the continuous eigenvalue over the domain k ∈ [−π, π]. This defines the continuum in the range E k ∈ [−2, 2] with the associated density of states given by The divergences occurring at the band edges (continuum threshold) E = ±2 are the van Hove singularities characteristic of 1-D systems [100][101][102] .
Next we aim to obtain the discrete spectrum associated with the quantum dot by solving the Schrödinger equation in its vicinity. To do so we first introduce the projection operator P = |0 0| + |d d| for the central region consisting of the dot |d and the chain site |0 to which it couples. We further introduce the operator Q = 1 − P that represents the system environment. We then project out the Q sector according to the Feshbach method 103,104 to write the Schrödinger equation in the P sector as where we have applied Eq. (2). Here, the effective Hamiltonian H eff (E) is given by see Appendices of Refs. 105,106 . We emphasize that H eff depends on E through Eq. (4). Now, taking the discriminant of Eq. (6), we obtain the dispersion equation for the discrete eigenvalues in the form where Σ(E) has been written specifically in the first Riemann sheet and can be analytically continued to the second sheet by the usual methods. After squaring to remove the root, we obtain an equivalent quartic polynomial condition for the eigenvalues as p(E j ) = 0, in which Let us point out an immediate consequence of Eq. (9). Notice that if we take d = −2, the polynomial becomes p(E) = (E + 2) 3 (E − 2) − g 4 = 0. If g < 1 is small, this results in three eigenvalues clustered near the lower band edge at E = −2. We present these eigenvalues as they appear in the complex energy plane for the case g = 0.5, d = −2 in Fig. 1 (a). We further show the associated value of the complex wave vector k j , obtained from the dispersion E j = −2 cos k j , for each eigenvalue in Fig.  1 (b). Notice that one of the three clustered eigenvalues (E B− ) has Im k B− > 0, for which the wave function form Eq. (2) yields a normalizable solution, implying that E B− represents a bound state appearing in the first Riemann sheet of the complex energy plane. The other two among the clustered eigenvalues are a resonance E R and an anti-resonance E A with complex conjugate eigenvalues and non-normalizable wave functions (these reside in the second sheet). Finally, the fourth eigenvalue E B+ is also a bound state but instead appears near the upper band edge. However, this fourth solution plays little role in the physics in this situation, and hence we can mostly ignore it in what follows.
From the form of the polynomial p(E) = (E + 2) 3 (E − 2) − g 4 = 0 at d = −2, we can also easily see that in the limit g → 0, the three clustered eigenvalues must converge on the lower band edge. This behavior is illustrated diagrammatically in Fig. 2. Further, we can obtain an expansion for the eigenvalues in the vicinity of the band edge (threshold) for small, nonzero g values by substituting an ansatz of the form E B− = −2 + χ α g α + χ β g β + . . . into the dispersion Eq. (9) and solving for the coefficients χ s and the exponents. Doing so, we find an expansion for the bound-state eigenvalue as well as the resonance and anti-resonance eigenvalues The resonance width Γ = 2Im E R has the predicted form Γ ∼ g 4/3 , which is the Markovian aspect of the Purcell effect discussed in Sec . I.

B. Eigenstates: Mapping to the quadratic eigenvalue problem
While the preceding formalism has yielded the energy eigenvalue spectrum for the model in Eq. (1), notice that it is not quite adequate to fully describe the behavior of the three apparently converging eigenstates in the g → 0 limit, because in Eq. (7) the effective Hamiltonian is written as a function of its own eigenvalue, which means that any two differing eigenstates technically belong to the solution space of two different copies of H eff . In other words, Eq. (6) represents a nonlinear eigenvalue problem (See Sec. II. C. of Ref. 46 for a similar discussion). To address this problem, notice that we could equivalently write the effective Hamiltonian Eq. (7) in terms of the variable λ = e ik . This is a convenient choice because we can use E = −2 cos k to write the energy variable appearing on the right-hand side of Eq. (6) as Taken together, this allows us to rewrite Eq. (6) as an equivalent quadratic eigenvalue problem in λ in the form of the coupled equations and − gλ 0|ψ + λ 2 + d λ + 1 d|ψ = 0.  In general, the quadratic eigenvalue problem can be solved after mapping to a generalized linear eigenvalue problem 107,108 . In the present case, we do so by rewriting the paired equations (13) and (14) in the form in which and the eigenvector |Ψ takes the form see Ref. 106 . The four discrete eigenvalues λ j are obtained in the present context from the discriminant of Eq. (15), which yields the quartic f (λ j ) = 0 with This is equivalent to the quartic for the energy eigenvalue p(E j ) = 0 in Eq. (9). The norm of the corresponding eigenstates is fixed by the normalization condition Ψ j |B|Ψ j = 1, which gives Combining this with Eq. (13) we obtain the conditions and With the present formalism in hand, the four eigenstate solutions |Ψ j from Eq. (15) are now in one-to-one correspondence with the four discrete eigenvalues from f (λ j ) = 0. This will enable us to study precisely the behavior of the eigenvectors as we approach the triple degeneracy at the threshold in Sec. IV B and App. A as well as the influence of the degeneracy on the survival probability of the occupied dot state in Sec. V. Before turning to these issues, however, we first examine the eigenvalue spectrum in the vicinity of the threshold in greater detail in the next section.

III. EIGENVALUE SPECTRUM IN THE VICINITY OF THE THRESHOLD
In the previous section, we saw a snapshot of the spectrum for the case that d coincided with the lower band In this section, we relax this condition and allow d to vary in the immediate vicinity of the band edge in order to illustrate the nearby presence of three ordinary exceptional points. We will locate these EPs by applying the method from Ref. 109 . In Fig. 3(a) we plot the real part of the three eigenvalues that appear near the lower band edge for the value The solid curve here indicates the bound state E − . For values of d < −2, this bound state appears shifted slightly below the unperturbed energy d itself. However, for d > −2 this eigenvalue breaks off and gets "stuck," always appearing below the lower band edge. Meanwhile, the two other eigenvalues are shown by a dashed curve. When d appears well below the continuum at d <¯ 0 ≈ −2.05518 . . . these solutions appear as two virtual bound states, which are non-normalizable states with real eigenvalue in the second Riemann sheet 49,81,92,106,110,111 However, as d draws closer to the continuum at d =¯ 0 the virtual bound states coalesce before forming a resonance, anti-resonance pair with complex conjugate eigenvalues. This shows that d =¯ 0 is an (ordinary) EP2.
We note that the above picture for the real part of the eigenvalues [see Fig. 3(a)] is reminiscent of the traditional avoided crossing picture in quantum systems; however, the resonance and anti-resonance actually collide while the bound state experiences an avoided crossing with the other two solutions. We will momentarily obtain deeper insight into this picture.
The corresponding imaginary parts of these eigenvalues are then shown in Fig. 3(b), in which we immediately see that the resonance decay width −ImE R ∼ g 4/3 for d ¯ 0 is enhanced near the threshold d ≈ −2 compared to the case as d becomes embedded deeper inside the continuum where the Fermi golden rule would apply. This indicates that the relaxation process should be enhanced near the continuum threshold, as we would expect from the Purcell effect. However, to get a full picture for the decay dynamics we will also have to take into account the branch-point effect; this is worked out later in Sec.

V.
Now we aim to obtain precise information about the location of the EP2 at which the resonance and antiresonance solutions appear. In the vicinity of any exceptional point we can expand the associated eigenvalues in the characteristic Puiseux expansion 1,50,109 , which in the present case we expect takes the form E j = E + α j ( d −¯ 0 ) 1/2 + . . . withĒ being the coalesced eigenvalue and¯ 0 the location of the exceptional point itself. Notice from this expression that the derivative of the eigenvalue with respect to the parameter d blows up at the exceptional point as We make use of this fact to extract information about the EP in the following quick derivation. First we take a full derivative of Eq. (8) and rearrange to obtain Then letting d →¯ 0 (so that E →Ē) we obtain a condition on the coalesced eigenvalue as This double-cubic polynomial yields three positive solutions (associated with the upper band edge) and three negative solutions (associated with the lower band edge). Maintaining our focus on the lower band edge, we obtain the exact form of the three negative solutions as for n = −1, 0, 1. The solutionĒ 0 corresponds to the previously-discussed exceptional point d =¯ 0 at which the resonance and anti-resonance pair appear [see Fig.  3(a)]. MeanwhileĒ 1 andĒ −1 are two complex-valued solutions of Eq. (24) that can only be seen by extending the problem into the complex d parameter space, as shown in Fig. 4. Notice that whereas it was the resonance and anti-resonance that coalesced at the real-valued EP2¯ 0 corresponding toĒ 0 , it is instead the complex-extension of the bound state that coalesces with one of the other two eigenvalues at the complex-valued EP2s corresponding toĒ 1 andĒ −1 , while the third eigenvalue experiences an avoided crossing. (We note this configuration of three EP2s is also studied from a pure mathematics perspective in a recent work, Ref. 112 .) Also notice from Eq. (25) that for g → 0 the three ordinary EP2s converge to the band edge at E = −2. This suggests that in the g → 0 limit the three ordinary EP2s collectively give rise to the anomalous threshold EP that was previously introduced at the end of Sec. II A, and which is the main topic of this work. From this point forward, we focus our attention on the anomalous EP itself.

IV. PROPERTIES OF THE ANOMALOUS THRESHOLD EXCEPTIONAL POINT
For the remainder of the paper, we focus directly on the properties of the anomalous exceptional point at the lower band edge. Hence, we now assume the parameter condition d = −2 is always satisfied, which corresponds to the vertical red line in Fig. 3.

A. Small g behavior of the eigenstates near threshold
Previously, we obtained Puiseux expansions for the three converging energy eigenvalues in Eqs. (10) and (11), which describe their behavior for small values of g at d = −2. Next we obtain similar expansions for the λ eigenvalues introduced in Sec. II B as well as the corresponding norm of the converging eigenstates. These expressions will prove useful in analyzing the time evolution near the threshold later in Sec. V, while also providing insight into the mathematical properties of the anomalous exceptional point.
First, similar to the process that we used to obtain Eqs. (10) and (11), we can obtain Puiseux expansions for the λ eigenvalues for the bound state and resonance and anti-resonance solutions after applying a generic form of these expansions in the λ polynomial Eq. (18).
As our next step, we can apply the above expressions in Eq. (21) to obtain expansions for the d| component of the norm in the vicinity of the anomalous EP. For the bound state we obtain and for the resonance, anti-resonance pair we find Notice in Eqs. (28) and (29) that the component of the norm d|ψ B−,R,A diverges for all three states in the limit g → 0, which conforms with the expected behavior of an exceptional point 46,48,49 . In particular, it again seems to suggest that all three of these states are converging to the dot state |d in this limit. However, as we will show next, connecting the small g behavior of the system with the actual limit g → 0 is a more subtle issue than what this picture suggests.
B. Jordan block structure in the g → 0 limit We now obtain the explicit form of the eigenvalue problem (15) in the limit g → 0. To simplify our discussion somewhat, we first rewrite the generic form of the eigenvalue equation slightly as which takes the form of an ordinary linear eigenvalue problem of the non-Hermitian matrix B −1 A. In the limit of interest g → 0 with d = −2, the non-Hermitian matrix B −1 A takes the form It is then straightforward to show that this matrix can be transformed to the simplest form which contains a 2 × 2 Jordan block, plus an additional degeneracy (the form of R is reported in App. A). This unambiguously demonstrates that the exceptional point at the threshold is in fact an EP2, despite essentially behaving like an EP3 in the vicinity of the coalescence. More precisely connecting the EP3-like behavior for small g with the actual EP2 in the g → 0 limit can be accomplished by realizing that the three states |Ψ B− , |Ψ R , and |Ψ A no longer really exist individually in this limit. Instead, certain linear combinations can be taken to appropriately connect the spectrum in the two cases. However, since the details are somewhat tedious and not strictly necessary for the proceeding analysis, we leave these details to App. A.

V. NON-MARKOVIAN PURCELL EFFECT AND INFLUENCE OF THE EP
Having established the spectral and mathematical properties of the anomalous exceptional point, we now analyze its influence on the relaxation process of the excited quantum emitter state. A number of works have established that exceptional points can influence spontaneous emission in a variety of circumstances 46,63,[68][69][70][71][72][73][74][75][76] . Below, we demonstrate that the combined influence of the anomalous exceptional point and the threshold itself determine the decay dynamics of the dot during the most consequential timescale.
We assume that the system is initially in the state in which the dot |d is fully occupied. The survival probability P (t) = |A(t)| 2 is given as the square modulus of the survival amplitude A(t) = d|e −iHt |d , which in the context of the present formalism is written as a sum over the contributions from each eigenstate in the form The contour C is shown in Fig. 5(a). Note that this expression is derived for a generic tight-binding model in Ref. 106 and for the present model in Ref. 44 .
In Fig. 6(a) we show the evolution for the case g = 0.02, very near the threshold EP. We note that the dynamics can be broken down into several regimes: the earliest (and very brief) quantum-Zeno timescale t 1/| d | yields the usual short-time parabolic decay of the form P (t) ≈ 1 − g 2 t 2113-116 . This is quickly followed by an intermediate timescale 1 < t 1/g 4/3 , during which most of the decay occurs; it is during this period that the influence of the anomalous exceptional point on the dynamics is most pronounced, as detailed below. A dip then appears in the probability that transitions into a decaying oscillatory phase, which in turn gradually settles asymptotically to an incomplete or fractional decay.
The fractional decay is due to population trapping in the bound state E B− near the lower threshold. Note that the influence of the upper bound state E B+ is negligible in this situation.

A. Intermediate-time scale dynamics
To evaluate the exceptional point influence on the dynamics on the key timescale, we focus on the contributions to the integral in Eq. (33) coming from the three coalescing states j = B−, A, R. Changing the integration variable from λ to E = −λ − 1/λ, we can rewrite the outer contour in Fig. 5(a) as an integration over the branch cut in the energy plane, as in Fig. 5(b). Then after neglecting the contribution from the upper bound state B+ the total survival amplitude can be written in which the branch cut contribution now takes the form and the pole from the lower bound state is evaluated as In Eq. (35) the contour C E surrounds the branch cut in the complex E plane as shown in Fig. 5(b). At this stage, we might consider applying the expansions for the eigenvalues near the EP from Sec. IV to evaluate the branch cut contribution A bc (t) (which would be similar to the approach from Ref. 46 ). However, doing so, we find evaluating the pole at the band edge to be a challenge. Instead, we rewrite the pole appearing in A bc (t) as a second integration in the form Again transforming the integration variable from the E plane to the k plane according to E = −2 cos k allows us to rewrite the inner integration in terms of a Bessel function J 1 (x) as We note that the information about the branch cut is now encoded in the Bessel function J 1 (2t). Making a final integral transformation to the variable t = t − τ one can show that the resonance and anti-resonance contributions to A bc (t) take the form 44,46 for j = {R, A}, while the bound state contribution instead evaluates as Finally, combining Eqs. (40) and (41) with the pole contribution in Eq. (37) we obtain the total survival amplitude [Eq. (34)] in the form in which the contributions from the three states involved with the EP all appear on an equal footing. This enables the following analysis.
We are now in a position to apply our near-threshold expansions to study the influence of the EP on the evolution during the key timescale. Introducing the small parameter β ≡ g 2 /2 1/3 we can treat β 2 t as the smallest non-negligible quantity during the timescale 1 < t 1/β 2 (in which 1/β 2 ∼ g −4/3 ). For the first term in Eq. (42) we can then expand the exponential and apply the expansions near the EP from Eqs. (10), (11), (28) and (29) to find The eigenvalue expansion has yielded in the second term the factor t 2 e 2it , which one would expect in the vicinity of an effective third-order exceptional point. However, this term will ultimately cancel in the final expression for the survival amplitude. Next, we perform a similar expansion for the second term in Eq. (42), and then find it useful to rewrite the resulting integrals over the Bessel function according to We then use the formulae and finally expand the Bessel functions for t > 1 according to to find the approximation for the second term of Eq. (42) as Notice that the first term of this expression will cancel with the second term of Eq. (43), as predicted. Further, note that the second term here in (50) contains the factor t 3/2 e 2it , which can be understood as itself consisting of two factors, one from the exceptional point and the other from the threshold. The former is the factor t 2 e 2it coming from the eigenvalue expansion near the effective thirdorder EP that we have already encountered. The second factor t −1/2 appeared from expanding the integral over the Bessel function in Eqs. (48) and (49), which is the contribution from the band edge, similar to Refs. 46,81 . Putting the two terms Eqs. (43) and (50) together, we obtain our approximation for the time evolution integral Eq. (42) as after again replacing β = g 2 /2 1/3 with the physical parameter g. Finally, the resulting approximation for the survival probability itself is given by This approximation is shown by the dashed curve in Fig.  6 (b), which we see captures the dynamics quite well during the period in which the majority of the decay occurs.

B. Long-time scale dynamics
Next, we evaluate the dynamics in the latter stages of the evolution. As a general rule-of-thumb, we expect that the branch point effect will come to dominate the dynamics as time progresses 46,81,99,[117][118][119][120][121] , which often leads to inverse power law decay as t → ∞. To evaluate this, we return to the expression for the survival amplitude near the threshold A(t) ≈ A bc (t) + A B− (t) from Eq. (34). We evaluate the branch-cut term A bc (t) by deforming the integral contour C E in Eq. (35) by dragging it into the lower half of the complex energy plane as shown in Fig. 5(c). In doing so we pick up a pole at the resonance eigenvalue, which is evaluated 1 near the thresh- ) after applying the expansions Eqs. (29) and (11). Similarly, we can expand for the bound-state pole from Eq. (37) ). The remaining integration contour in Fig. 5(c) extends from the branch points at E = ±2 out to infinity in the lower half plane; this contribution yields the typical inverse power law decay with t −3/2 dependence in the amplitude (see also Refs. 81,99 ). Putting these three pieces together the survival probability in this case follows This is shown as the light blue dotted curve in Fig. 6 (c), which focuses on the late-time evolution of the decay. We note that this approximation picks up around the first minimum, just as the oscillatory dynamics kick in. The three cross-terms in Eq. (53) each give a decaying oscillation with frequency ∼ g 4/3 . However, the resonance contribution with lifetime Γ −1 ∼ g −4/3 [about ∼ 184.2 for g = 0.02 as in Fig. 6 (c)] has mostly died out after the first oscillation, so that the dominant contribution comes from the bound-state/power-law cross term, with oscillation frequency E B− +2 = g 4/3 /2 2/3 . Finally, in the long-time limit, the dot survival probability settles down to the asymptotic occupation probability for the bound state, given simply by P ∞ = lim t→∞ |A B− pole (t)| 2 = 4/9. This is shown by the horizontal black line in Fig. 6 (a).

VI. DISCUSSION
In this work, we have shown the existence of an anomalous-order exceptional point at threshold (band edge) in a simple 1-D continuum model with an attached quantum emitter and examined its influence on the decay dynamics when the emitter is tuned to the threshold energy. Although we have shown this in the specific model in Eq. (1), we claim that it is a much more ubiquitous effect. Let us consider a generic model for a simple open quantum system of the form consisting of the quantum emitter q † coupled to the continuum mode c † k . Note that our original Hamiltonian Eq. (1) can be placed in this form after applying a simple Fourier transform. Here the continuum E k ranges from a lower threshold E th− to an upper threshold E th+ (the latter of which might appear at infinity). We claim that there are two conditions to be placed on this model in order for the anomalous threshold EP to be realized. These conditions are most naturally expressed after writing the self-energy function Σ(z) as in which ρ(E) = ∂k/∂E is the density of states (DOS) function and we have re-written v(E) = v k(E) by inverting E k . (We note the appearance of the so-called reservoir structure function ρ(E) |v(E)| 2 in this expression, as previously defined in the literature 88,89 .) The first of the two conditions is that the density of states must contain the square-root divergence at either threshold E th , as occurred for our specific model in Eq. (5). This is indeed the standard form for the van Hove singularity in 1-D systems [100][101][102] . Then the second condition appears on v(E), which simply requires that this quantity should be non-singular at the threshold E = E th . Assuming these two conditions hold, then as shown in App B the self-energy function will reproduce the square-root divergence, as occurred for our specific model in Eq. (8), and hence the spectrum will contain the triple-level convergence illustrated in Fig. 2 with energy eigenvalues that yield the characteristic g 4/3 Puiseux expansion, similar to Eqs. (10) and (11) in the main text [in the case of the general model, the Puiseux expansion follows directly from Eq. (B8)]. Physically, the quantity g 4/3 appears in the inverse period of the oscillations in Fig. 6. We mention that some quantities, such as the norm in Eqs. (28) and (29), exhibit a slightly more general expansion in terms of g 2/3 . Illustrating this universality, we note that a similar picture for the dynamics shown in Fig. 6 has previously appeared in the literature [82][83][84] in the context of spontaneous emission from atoms near a band edge within photonic band gap (PBG) materials. In Ref. 84 the authors note "the peculiar power of 2 3 being one of the signatures of the unconventional nature of the PBG environment." In the present work, we have revealed that the exceptional point is the underlying mechanism that shapes the dynamics near such a band edge, and that its Puiseux expansion is the origin of the "peculiar power of 2 3 ." (See Refs. 77,93,94,98,128 for other models in which a similar expansion appears.) Indeed, in Refs. 83,86 the authors highlight the existence of a "doublet" consisting of an atom-photon bound state and a resonance when an atom is located near the photonic band edge; however, the presence of an anti-resonance is not noted. We have shown that it is the eigenvalue convergence and (partial) coalescence involving all three eigenstates that underlies the physics in this situation.
While the anomalous-order exceptional point should be a fairly universal feature, it should be considered that in this work we have analyzed the problem only at the level of the Hamiltonian dynamics, which cannot account for processes like quantum jumps. Several quite recent works have revealed that the parametric location and dynamical characteristics of the exceptional points can become modified when treating the problem at the level of the Liouvillian formalism with Lindblad terms that can account for such processes 16,64,76,122,123 . However, recent experiments in circuit quantum electrodynamics have to some extent circumvented this issue through the process of data post-selection, in which trials that result in quantum jumps are eliminated from the final analysis 124,125 . We propose that experimental observation of the features of the anomalous exceptional point might be achieved by a modified version of these experiments in which a superconducting qubit is embedded in a waveguide (instead of a cavity as in the original experiments 124,125 ) with the qubit transition frequency tuned to the lowest waveguide cutoff mode.
We close the paper with the following two comments aimed at placing our results in the context of the wider literature.

A. Comment on the anomalous order of the EP
We emphasize that one peculiar mathematical aspect of the present results is the sharp distinction in the nearthreshold spectral properties of the system in the case of small coupling as opposed to the case in which the coupling actually vanishes. In the former case, the system behaves precisely as if there were an EP3 at the threshold, whereas in the latter case there instead appears an exact EP2. To connect these seemingly incongruent pictures, we showed in App. A that a reorganization of the eigenstates occurs in the limit in which the coupling vanishes such that two linear combinations of the three relevant eigenstates converge on the state of the decoupled quantum emitter, while a third linear combination instead merges with the continuum at the threshold. We mention that this resolution of the problem was partly inspired by Ref. 51 , in which Hashimoto, et al, introduced a novel basis near an exceptional point such that the system eigenvectors can connect continuously with the Jordan block representation directly at the EP, which we referred to in our internal discussions as Hashimoto's representation.
It is interesting to note that the mismatch between the Puiseux expansion and the exact order of the EP in the present case is rather different than what has appeared previously in the literature. First, in the most typical case, in the vicinity of an order N exceptional point the relevant eigenvalues can be expanded in a Puiseux series written in terms of an N th order root. However, in Ref. 54 the authors point out the existence of a special case in which the eigenvalues are instead organized into distinct subgroups (called cycles), which have their own separate Puiseux expansions that are lower order than N . For example, they demonstrate that an EP3 can occur for which two coalescing eigenvalues can be expanded in a square root near the EP3 (order 2 Puiseux expansion), while the other coalescing eigenstate instead takes the form of a Taylor series (which can be viewed as an order 1 Puiseux expansion). Notice then that the order of the two cycles adds up to the order of the EP itself 3 = 2 + 1. We emphasize that here the order of the individual Puiseux expansions is always lower than the order of the EP.
In contrast, in our case not only is there a mismatch between the order of the EP and that of the root appearing in the Puiseux expansion, but the Puiseux expansion is instead higher order than the EP itself. Further, this mismatch is resolved in a fundamentally different way: instead of subgroups of eigenvalues with their own Puiseux expansions, the eigenstates in our case form subgroups with different coalescing behaviors (two linear combinations of eigenstates coalesce on the discrete quantum emitter while a third instead merges with the continuum).
The resolution to the mismatch in the present case points to the reason that a more general type of exceptional point is allowed in our model in the first place.
The key point is that in Ref. 54 the authors only consider finite, non-Hermitian matrices, in contrast to our model, which incorporates an eigenvalue continuum. Hence we conclude that the presence of continuous eigenvalue spectra in a given Hamiltonian permits a wider diversity of exceptional point types.

B. Wider context of near-threshold non-Markovian dynamics in 1-D systems
While here we have analyzed a situation in which the eigenvalues of three states converge directly on the continuum threshold, we note that the dynamical influence of individual states appearing near the threshold has previously been studied in the literature in 1-D systems. In Refs. 126,127 it has been shown that a resonance with a near-threshold eigenvalue can lead to full-time nonexponential decay. Meanwhile, in Refs. 81,89,90,111 it is shown that a virtual bound state (or anti-bound state) merging directly with the continuum threshold can also result in full-time inverse power law decay. More generally, the virtual state appearing in the vicinity of the threshold introduces a timescale that characterizes the power law decay 81,111 . This timescale is inversely proportional to the energy gap between the virtual state eigenvalue and the threshold, and it divides the dynamics into an intermediate-time zone during which the survival amplitude (survival probability) follows 1/t 1/2 (1/t) decay, and a long-time zone in which the amplitude (probability) instead falls off as 1/t 3/2 (1/t 3 ) 81,111 . (We emphasize that a localized bound state near the threshold introduces the same timescale as the virtual bound state, but in this case the decay dynamics are to some extent obscured by the incomplete decay resulting from the bound state itself 81 .) Meanwhile, the dynamics associated with multiple (N ) eigenstates that coalesce at a real eigenvalue that itself appears near the threshold is fundamentally different from, yet also related to the above cases. In Ref. 46 the case of two virtual states coalescing before forming a resonance/anti-resonance pair is considered. In that paper, when the eigenvalue at which the two states coalesce is near the threshold the picture is somewhat similar to that of the lone near-threshold virtual state, except that the the intermediate timescale dynamics is replaced with a decay of the form 1 − Ct 1/2 46 . Meanwhile, in the present work, we have shown that a triplet of resonance, anti-resonance and bound eigenstates converging directly at the threshold results in decay on the intermediate timescale of the form 1 − Ct 3/2 . In either case, the intermediate dynamics are replaced by the typical 1/t 3 decay on the asymptotic timescale, which matches with the single virtual state case. Attempting to generalize from this picture, we propose that the intermediate timescale dynamics for N converging eigenstates with real eigenvalue near the threshold might be expected to contain a term with the time dependence ∼ t N −3/2 . This can be understood as appearing after taking the square modulus of the amplitude containing the terms (1 − const 1 t 1/2 × t N −1 )e −iĒ EP t , in which the factor t N −1 e −iĒ EP t in the second term arises from an N th-order pole appearing due to the N converging states, and the factor 1/t 1/2 results from the influence of the continuum threshold, just as in the case of the single near-threshold virtual state discussed in the previous paragraph. However, for the anomalous-order EP in the present case, the number N seems to be equal to the number of converging levels (N = 3) rather than the technical order of the EP, which suggests that one might have to consider the situation in which the EP occurs directly at the threshold on a case-by-case basis.
The results for the near-threshold dynamics in 2-D and 3-D systems reported in Refs. [95][96][97] suggest that one probably must consider these cases separately, as well. Another interesting case that one could consider would be that of multiple eigenstates converging on a complex eigenvalue near a threshold, which would also likely be modified somehow from the above scenario.
In this appendix, we continue our development from Sec. IV B, illustrating explicitly how to connect the 2 × 2 Jordan block appearing in the g → 0 limit with the EP3like behavior for g = 0.
First we report the explicit from of the so-called rotation matrix R appearing in Eq. 32 and are the generalized eigenstates at the threshold EP. Here |Ψ + in the first line represents the g → 0 limit for the lone bound state that previously appeared above the upper band edge in the case g = 0. For g → 0, this state takes the eigenvalue λ + = −1. Comparing the explicit eigenstate |Ψ + = {−1, 0, 1, 0} T with the generic form Eq. (17), we see that this state is completely decoupled from the dot as d|ψ + = 0. This indicates that |Ψ + has merged with the continuous spectrum in the g → 0 limit. The three remaining states |Ψ d , |Φ d and |Ψ − are those that each share the eigenvalue λ j = 1. This includes |Ψ d = {0, 1, 0, 1} T , which is the uncoupled dot state, and |Φ d = {0, −1, 0, 0} T , which is the associated pseudo-eigenstate. Together, these two states satisfy the Jordan-chain relations Finally, |Ψ − = {1, 0, 1, 0} T again represents a state that is completely decoupled from the dot and has joined the continuous spectrum. While |Ψ + merged with the continuum at the upper band edge with eigenvalue λ + = −1 (E = 2) , |Ψ − is instead merged at the lower band edge with eigenvalue λ − = 1 (E = −2).
However, the observant reader may take some discomfort at this point, noting that it is not immediately clear how the three states |Ψ R , |Ψ A and |Ψ B− involved in the convergence on λ = 1 for small g from Sec. IV connect with the three (generalized) states |Ψ d , |Φ d and |Ψ − appearing in the actual g → 0 limit. In particular, the state |Ψ − cannot be directly connected with the limiting behavior of the lower bound state |Ψ B− in the same way that |Ψ + simply emerged as the limit of the upper bound state. Indeed, Eq. (28) seems to suggest that |Ψ B− should just participate in the eigenvalue coalescence with the other two states.
The resolution to this issue comes in the realization that the states |Ψ d , |Φ d and |Ψ − can only be obtained as the limit of specific linear combinations of the bound state, resonance state and anti-resonance state. This approach borrows conceptually from Hashimoto's representation (introduced in Ref. 51 ), in which an extended Jordan block representation is written such that the eigenvectors away from an exceptional point can connect continuously with those appearing directly at the EP (see also Sec. V of Ref. 48 ).
As our initial step, we note that in the case d = −2 we can use Eq. (14) to write gλ j 0|ψ j = (1 − λ j ) 2 d|ψ j . Combining this with Eqs. (28) and (29), we can rewrite (17) as in which α j = 0 for the bound state (with j = B−), α j = −1 for the resonance (j = R) and α j = 1 for the anti-resonance (j = A). Using this expression, we first obtain the dot state from the most straightforward linear combination To obtain the partner pseudo-dot state, however, we need to be a little more clever. From the following linear combination, we can engineer the cancellation of all the leading-order entries in Eq. (A5), so that we get Note that this pseudo-state is actually different from |Φ d appearing in Eq. A2. In fact, this pseudo-state satisfies the Jordan-chain relation  In Sec. VI we claimed for the generic model given in Eq. (54) that two conditions are required for the triplelevel convergence in Fig. 2 and anomalous exceptional point from the main text to occur. One of these conditions was that the continuum must be 1-D, which yields the characteristic van Hove singularity (E th − E) −1/2 . The second condition was that the potential v(E) must be non-singular at the threshold v(E th ). We motivate these claims in what follows.
We briefly note that we have applied the rotating-wave approximation in writing our generic model in Eq. (54). However, Ref. 77 provides an example of a slightly more general model in which counter-rotating terms are retained and yet the threshold EP still occurs.

Form of self-energy for the generic model
Here we quickly outline the conditions such that the self-energy function from Eq. (55) for the generic model reproduces the square-root divergence in the denominator from the DOS, which in turn yields the anomalous EP. We start by rewriting the self-energy function in terms of an integration in the complex wave vector plane as We assume that Σ(E) is defined first for real E appearing below the continuum threshold E < E th (i.e. a bound state). Afterwards the result can be analytically continued into the complex domain. Next we assume a generic dispersion relation of the form E k ∼ k 2 + E th , which might be exact or merely an approximation near the threshold E th . Then we rewrite Σ(E) in the form in which k ± = ±i √ E th − E. Since the integrand of this expression vanishes like 1/k 2 as k → +i∞, we can close the integration contour in the upper half of the k-plane, so that taking the residue at k + gives Thus the self-energy function obtained by integration over the 1-D continuum has recreated the square-root divergence from the DOS, which verifies our first condition.
For the second condition, notice that if v k+ contained some singularity at k + = 0 (E = E th ), this would have the effect of disturbing the square root in the denominator and thus the anomalous exceptional point would no longer appear (as occurs for the models in Refs. 90,92,111 , for example).
Note that a further term that is analytic at E th might arise depending on the precise form of v k or integration over a second continuum with a different threshold, for example. Hence, we can write the self-energy in Eq. (B3) in a slightly more general form as in which ∆(E) and λ(E) are both analytic in E. (Notice that this is now more general than the model from the main text.) One could reasonably suppose that the presence of the ∆(E) term here might disrupt the occurence of the triple-level convergence associated with the exceptional point for this more general case. We show in App. B 2 that indeed the triple-level convergence still occurs.

Generality of the triple degeneracy at threshold
Here we demonstrate that the form of the self-energy reported in Eq. (B4) indeed yields the triple degeneracy at threshold in the limit of vanishing coupling, as described in the main text and in Fig. 2.
The point spectrum for the generic model is obtained from the dispersion equation E − q = Σ(E), which corresponds to Eq. (8) for the model in the main text. To study the point spectrum in the vicinity of the threshold we specify q = E th in this equation and use the explicit form of Eq. (B4) to write Multiplying through by x ≡ √ E th − E we obtain a cubic equation in x of the form x 3 + g 2 ∆(E)x + g 2 λ(E) = 0, which is solved by . (B7) For small g 1, we find that ξ − ≈ −g 2 λ is of lower order than ξ + ≈ g 4 ∆ 3 /27λ, which results in Squaring, then cubing, we immediately see the familiar form (E − E th ) 3 = −g 4 λ 2 from the main text in Sec. II A, which yields the expected triple-level convergence as g → 0 as well as the Puiseux expansion in terms of g 4/3 . (See also Refs. 77,82,84,91,93,94,98,128 for specific models in which a similar expansion has appeared.)