Out-of-equilibrium phase diagram of long-range superconductors

Within the ultimate goal of classifying universality in quantum many-body dynamics, understanding the relation between out-of-equilibrium and equilibrium criticality is a crucial objective. Models with power-law interactions exhibit rich well-understood critical behavior in equilibrium, but the out-of-equilibrium picture has remained incomplete, despite recent experimental progress. We construct the rich dynamical phase diagram of free-fermionic chains with power-law hopping and pairing, and provide analytic and numerical evidence showing a direct connection between nonanalyticities of the return rate and zero crossings of the string order parameter. Our results may explain the experimental observation of so-called \textit{accidental} dynamical vortices, which appear for quenches within the same topological phase of the Haldane model, as reported in [Fl\"aschner \textit{et al.}, Nature Physics \textbf{14}, 265 (2018)]. Our work is readily applicable to modern ultracold-atom experiments, not least because state-of-the-art quantum gas microscopes can now reliably measure the string order parameter, which, as we show, can serve as an indicator of dynamical criticality.


I. INTRODUCTION
Criticality and universality in equilibrium 1-3 are well-established concepts thanks in large part to the renormalization group. 4,5 Universality classes in out-ofequilibrium physics have also been developed for classical systems. 6,7 Nevertheless, criticality and universality in out-of-equilibrium quantum many-body systems are still much less understood concepts that are currently at the forefront of research efforts in condensed matter and coldatom physics.
Among the prominent thrusts of this research effort lies the phenomenon of dynamical phase transitions, 8,9 which fall into two major categories when considering sudden quenches of a quantum many-body system. The first is characterized by a local order parameter, whose long-time behavior determines the dynamical phase of the steady state. [10][11][12][13][14][15][16][17][18][19][20][21][22][23][24][25][26] In principle, this type of dynamical phase transition separates a ferromagnetic from a paramagnetic steady state in the wake of a quench, and this has been recently observed experimentally in trapped-ion setups. 27 However, in cases where the system has no finitetemperature phase transition and, therefore, always ends up in a paramagnetic steady state in the infinite-time limit, this dynamical phase transition can still be defined based on the decay behavior of the order parameter. 24,25,[28][29][30] For small enough quenches starting in an ordered state, the order parameter would decay asymptotically to zero without ever crossing zero. On the other hand, when the quench is large enough, the order parameter exhibits zero crossings as a function of time while forming an envelope that decays to zero in the long-time limit. Very interestingly, however, even when the model has no finitetemperature phase transition, the order parameter can exhibit persistent order at moderate-to-long times in the case of sufficiently small quenches in the presence of longrange interactions, 19 which has been shown to be due to domain-wall binding and confined dynamics. 24,31 The second type of dynamical phase transition is based on nonanalyticities in the so-called Loschmidt return rate 32 (see Sec. III for an overview), which can be construed as a dynamical analog of the thermal free energy with complexified evolution time standing for inverse temperature. Consequently, nonanalyticities in the return rate are construed as dynamical quantum phase transitions (DQPT) that occur at critical evolution times after a quench. In the few years since its introduction, 32 DQPT has witnessed a flurry of research activity from the perspective of both theory 9,33 and experiment. 34,35 Even though the seminal work of Ref. 32 showed in the case of the nearest-neighbor transverse-field Ising chain (TFIC) that nonanalyticities in the return rate occur only upon quenching across the quantum equilibrium critical point, soon thereafter it was indicated that this is neither a necessary nor a sufficient condition in other integrable and nonintegrable models. 36,37 Indeed, it was then firmly established that the dynamical critical point, which separates different phases of DQPT, is in general distinct from the quantum equilibrium critical point and is strongly dependent on the initial condition. 18,21,22 Additionally, this dynamical critical point was shown to coincide with that of the type of dynamical phase transition based on a local order parameter. 18,[20][21][22][23][24][25] Also recently, the theory of DQPT has been extended to Floquet systems [38][39][40] and models with many-body localized phases. 41,42 On the one hand, free-fermionic models have been at the center of the theoretical and even experimental investigations of DQPT. Indeed, the connection of DQPT to the underlying topological phase transition 43 was first illustrated in Bogoliubov-de Gennes Hamiltonians, and one of the first experiments on DQPT implemented the Haldane model on a hexagonal lattice. 34 On the other hand, long-range interactions have proven to give rise to rich dynamical critical behavior, 18 just as they do for equilibrium criticality. 44,45 In particular, the DQPT picture qualitatively changes when long-range interactions are introduced. Domain-wall binding and confined dynamics at sufficiently small transverse-field strengths in long-range quantum Ising chains 24,31 lead to the appearance of anomalous cusps 18,20 in the return rate (see Sec. III), which would be entirely smooth otherwise. However, long-range quantum Ising chains are nonintegrable, and most DQPT studies on them have required advanced numerical treatment, such as the time-dependent variational principle (TDVP) within uniform matrix product states. [46][47][48] Consequently, the interest in free-fermionic models with long-range hopping or pairing has grown in recent years, as this allows the study of DQPT in an exactly solvable setting. Indeed, Dutta and Dutta 49 have recently investigated DQPT in the Kitaev chain with nearest-neighbor hopping and long-range pairing, while Defenu et al. 50 have employed a truncated Jordan-Wigner transformation on the so-called extended long-range Ising chain to map it onto a free-fermionic model of long-range hopping and pairing both scaling as 1/r α , with r the inter-spin distance and α > 1. In this paper, we investigate a Kitaev chain with generalized power-law hopping and pairing profiles and draw a connection between the different dynamical phases of DQPT and the dynamic behavior of the string order parameter (SOP). String order parameters are related to hidden symmetry breaking and distinguish between topologically ordered and trivial equilibrium phases. 51 Comparing the dynamics of the SOP after a quench to that of the return rate can thus provide a natural link between DQPT and the underlying equilibrium topological phase transitions.
The remainder of the paper is organized as follows: In Sec. II, we introduce our model and discuss its equilibrium phase diagram. In Sec. III, we provide a brief overview of the theory of DQPT. In Sec. IV, we derive the string order parameter for our model. An analytic treatment with which we construct the dynamical phase diagram of the Kitaev chain with general power-law hopping and pairing terms is then presented in Sec. V. Our results for the Loschmidt return rates and string order parameter dynamics are then discussed in Sec. VI. We conclude in Sec. VII.

II. MODEL
The model we consider is the generalized Kacnormalized 52 long-range Kitaev chain (LRKC) with powerlaw hopping and pairing profiles and closed boundary conditions. It is described by the Hamiltonian where are the hopping and pairing profiles, respectively, with Kac normalization factors Here and throughout N denotes the system size, α, β > 1 the power-law exponents, h the chemical-potential strength, and c j , c † j the fermionic annihilation and creation operators, respectively, on site j that satisfy the canonical anticommutation relations {c l , c † j } = δ l,j and {c l , c j } = 0. In (2) we have adopted the ring convention of Ref. 53, which translates into r ranging from 1 up to only N/2 − 1 rather than to N in (1). This is done in order to adequately deal with closed boundary conditions and effectively utilize the Fourier transformation into momentum space motivated by the translation invariance of (1). In the limit N → ∞, the choice between periodic or anti-periodic boundaries is arbitrary as they yield the same momentum-space Hamiltonian where the Fourier transforms of J r and ∆ r are respectively, and Li s (z) denotes the polylogarithm of order s and argument z. We recall here that we are working in the thermodynamic limit N → ∞, hence the Riemann Zeta functions ζ arise in these coefficients. The Riemann Zeta functions perform a Kac normalization that fixes the positive quantum equilibrium critical point of the model to h e c = 1. To solve for the spectrum of (5) one employs the Bogoliubov transformation where γ k , γ † k are the fermionic Bogoliubov annihilation and creation operators, respectively, with momentum k obeying the canonical anticommutation relations {γ k , γ † p } = δ k,p and {γ k , γ p } = 0, and The transformation (8) diagonalizes (5) into the form

III. DYNAMICAL QUANTUM PHASE TRANSITIONS
The theory of DQPT relies on an intuitive connection between the thermal partition function and the Loschmidt amplitude 32 where |ψ i is the initial state, H f is the quench Hamiltonian, and is set inconsequentially to unity throughout our paper. The Loschmidt amplitude (12) is a boundary partition function, and its logarithm is therefore proportional to a dynamical free energy: where N is the system size. Nonanalyticities in the thermal free energy as a function of temperature indicate critical temperatures at which thermal phase transitions occur. In (13), complexified time replaces the inverse temperature. Consequently, nonanalyticities of (13) as a function of time indicate critical times at which DQPT occur.
The type of nonanalyticities that occur in the wake of a quench can either be regular or anomalous. 18,20 The former are the cusps that occur for quenches crossing the dynamical critical point, while the latter can occur for arbitrarily small quenches within the ordered phase of the system, and they have been shown to be connected to local spin excitations forming the lowest-lying excitations in the spectrum of the quench Hamiltonian. 24,25,50 The advantage of the LRKC is that it is integrable and also allows for a closed-form expression of its return rate, which we briefly derive here. Let us denote by H i the initial Hamiltonian in whose ground state |ψ i we prepare our system, and by H f the final Hamiltonian with which we quench. Hamiltonian H i(f ) corresponds to (5) with h = h i(f ) , where here it should be clear that we quench the chemical-potential strength. Moreover, Hamiltonians H i(f ) are diagonalised by different Bogolioubov bases, which we denote, respectively, as (η k , η † k ) and (γ k , γ † k ), i.e., The respective ground states are the vacua |0 η and |0 γ . For the calculation of G(t) it is useful to write |0 η in terms of |0 γ by removing all η fermions from |0 γ as Here F is a normalisation factor. Following the BCS-theory 54 approach one can show that where the normalisation is given by F 2 = k>0 (1 + Λ 2 k ) and Λ k = tan(ϕ k ). The angle ϕ k stems from a Bogoliubov transformation between the pre-and post-quench basis with Substituting (16) into (12) we obtain Using the unitarity of e −iH f t , we write the product of the last two exponential terms above as where we have defined f k = 2ω f k . To arrive at the second line we have used that which follows from ω f k = ω f −k . Substituting (20) back into (19), and applying Pauli exclusion to the Taylor expansion of the remaining exponential terms, then yields where we have dropped an inconsequential phase arising from the e −itH f |0 γ term. Taking care of the commutation relations when multiplying out the above product, and writing the remaining terms in normal ordered form, allows us to simplify the expression to It is then straightforward to show that the return rate is Appendix A provides a comparison between the return rate (24) in the limit of α, β → ∞ and its counterpart in the TFIC The latter exactly maps onto the LRKC (5) in this limit. Immediately one sees that (24) can show nonanalytic behavior only in the momentum sectors k c that satisfy θ f kc − θ i kc = π/2 + nπ, with n ∈ Z, or, in other words, the mode k c is in an infinite-temperature state. One can further show that these critical momenta satisfy the implicit relation which in the case of nearest-neighbor couplings (α, β → ∞), reduces to As we shall discuss in great detail later, the relation (26) is satisfied if and only if h i and h f are on opposite sides of the dynamical critical point h d c ≤ h e c . Only quenches where the value of h f leads to at least one real solution of (26) can therefore lead to nonanalytic behavior in (24). When (26) is satisfied, one finds nonanalytic behavior in (24) at the critical times The TFIC can be mapped using a Jordan-Wigner transformation onto the LRKC model (5) in the limit α, β → ∞.
In this case, (26) reduces to (27) and yields a real unique value of k c if and only if h i and h f are on different sides of the quantum equilibrium critical point h e c , which in this case equals h d c . In the case of the general LRKC model, the situation is more subtle in that cusps from the topological phase to just above h e c can give rise to three critical momenta; 49 cf. Sec. V. More intriguingly, for certain regimes of (α, β), which we discuss in detail in Sec. V, doubly critical cusps -which are related to two distinct critical momenta -can arise for quenches within the topologically nontrivial phase. 50 As we elaborate in Sec. VI, these cusps may explain the so-called accidental dynamical vortices observed in the experiment of Ref. 34.

IV. STRING ORDER PARAMETER
The nonlocal form of the SOP makes it a useful quantity to use when studying transitions between equilibrium topological phases. den Nijs and Rommelse first showed that the SOP witnesses the hidden Néel order of the Haldane phase of the spin-1 AKLT model. 55 More recent work has extended the use of SOPs in the analysis of topological phases also to bosonic and fermionic systems. [56][57][58][59][60] Importantly, the SOP can also be experimentally observed in a quantum gas microscope. 61,62 Since the return rate is by definition also nonlocal, it is natural to ask whether the dynamical evolution of the SOP after a quench captures physics similar to that captured by the return rate. The recent work of Budich and Heyl 43 has indeed shown a link between these two quantities for a Kitaev chain with next-to-nearest-neighbor hopping and nearest-neighbor pairing interactions; For a quench from the topologically nontrivial phase into the trivial phase, the authors showed that the times at which zero crossings of the SOP occur are exactly the critical times of the return rate. As a natural extension of this work, we examine in the subsequent sections whether this correspondence also holds in the general Kitaev chain with long-range hopping and pairing (1), thereby extending the link between DQPT and the underlying equilibrium phases to this model. This is especially interesting because we know from quantum Ising chains that the presence of long-range interactions leads to a more subtle picture in the connection between DQPT and the longitudinal magnetization. For instance, for sufficiently small quenches, anomalous cusps appear even when the order parameter never changes sign. Also the regular cusps -which are connected to zero crossings of the order parameter -display slightly different periodicities compared to the cusps of the return rate at short to intermediate evolution times, as shown, for example, in the case of the Lipkin-Meshkov-Glick (LMG) model. 21,22 In this section we show how the SOP can be calculated numerically for this model. 63 We write the SOP defined in Ref. 43 as where φ ± i = c † i ± c i and the time evolution is generated by the post-quench Hamiltonian H f . Our aim is to numerically calculate i.e., the absolute value of the expectation value of (29) with respect to the ground state of the pre-quench Hamiltonian H i . Here we give a brief summary of how this is achieved. First, since (30) is defined with respect to |ψ i , it is useful to write (29) in normal-ordered form via Wick's theorem. A direct consequence of this theorem is that the vacuum (ground state) expectation value of a product of an even number of fermionic operators A i reduces to a sum of elementary contractions 64 Second, we note that the right-hand side of (31) is equivalent to the Pfaffian Pf(M ) of an anti-symmetric 2n × 2n matrix M with entries Since O lm (t) is a product of 2(m − l) = 2d fermionic operators, we then have that M(t) = |Pf(M )|. Consequently, for a given time t ≥ 0, calculating (30) is reduced to three steps: (i) Calculate the four distinct types of elementary con- where x is the largest integer smaller than or equal to x ∈ R.

V. DYNAMICAL PHASE DIAGRAM
Now we present one of the main results of our work, which is an analytic derivation of the dynamical phase diagram further augmented by numerical results.
The presence of two different power-law decays for the hopping and pairing matrices produces a very rich dynamical phase diagram in the LRKC (1). Indeed, while at equilibrium only two topological phases exist with equilibrium critical points h e c = (1, −1 + 2 1−α ) -which reduce to the nearest-neighbor case h e c = (1, −1) in the α, β → ∞ limit -the presence of long-range couplings produces novel dynamical phases that have no equilibrium counterpart. As already established, 33 quenching the chemical potential h across one of the equilibrium critical points produces a dynamical phase transition, which takes the form of regular 18 singularities in the system's return rate. In the TFIC the origin of these singularities can be traced back to the appearance of a perfectly athermal mode whose excitation probability is exactly p k = 1/2. For simplicity and without loss of qualitative generality, we shall henceforth restrict ourselves to nonnegative values of the chemical potential.
Conversely, different situations arise for general α and β where more than a single dynamical critical mode can occur for specific quenches. It has already been noticed that in the case of nearest-neighbor hopping and longrange pairing a triply critical phase can emerge with three excitation modes having p k = 1/2 for quenches across the topological phase transition. 49 More recently, the occurrence of a doubly critical dynamical phase in which two critical modes arise has been demonstrated for the case of long-range pairing and hopping with equal decay rate β = α for quenches within the topologically nontrivial phase, i.e., from h i < h e c to h f ∈ (h d c , h e c ). 50 It is worth noting that the properties of these peculiar dynamical phases can hardly be connected with those of equilibrium critical phenomena. Indeed, long-range couplings alter the universal behavior close to the equilibrium critical points only for α < 2 or β < 2. One can immediately see this by considering the equilibrium dynamical critical exponent z, where φ = min(α, β) and we have used the definition ω k ∝ k z at h = h e c . The product between the dynamical critical exponent z and the correlation length exponent ν always remains the same as in the nearest-neighbor interacting case, i.e., zν = 1. Note that the relevant region for long-range couplings in the LRKC is radically different with respect to the case of interacting field theories. 66 On the other hand, the doubly critical phase also arises in the case of irrelevant long-range couplings φ ≥ 2, demonstrating that the dynamical critical behavior of long-range systems is not as strictly tied to the equilibrium universality as in nearest-neighbor systems. 21 A similar discrepancy in the relevance of long-range couplings, first noticed in Ref. 50 for the β = α case, also exists in the present general case.
After a quench of the chemical potential h in the model (1), quasiparticles in certain momentum sectors k may attain infinite temperature with excitation probability p k = 1/2. The value of the final chemical potential h f at which the momentum k becomes critical after a quench starting at the chemical potential h i is given by 50 which is a mere rearrangement of (27), with J k and ∆ k given in (6) and (7). For simplicity, and without loss of qualitative generality, we shall restrict our analysis to quenches starting at h i = 0, which simplifies (35) to By definition, the equilibrium critical mode k = 0 becomes dynamically critical for quenches at h f (k = 0) = 1 = h e c . There also exist a mode k ∞ < π/2 that can never become critical and gives h f (k ∞ ) → ∞. Therefore, we can distinguish four different forms of h f (k) shown in Fig. 1: (i). The function h f (k) is monotonous and the system has only the traditional singly critical dynamical phase with h d c = h e c = 1 with just a single critical momentum. 32 (ii). The function h f (k) has a single minimum with value h d c < 1 such that for h d c < h f < 1 a doubly critical phase with two critical momenta emerges, whereas when h f > 1 a singly critical phase with a single critical mode appears. In this case, h f (k) has no local maximum.
(iii). The function h f (k) has a local minimum and a local maximum both at k > 0, where both the maximum and the minimum lie in the region h f > 1 such that h max > h min > 1. This gives rise to a singly critical phase for 1 < h f < h min or h f > h max , and a triply critical phase, which has three critical momentum modes, for h min < h f < h max . In this case, the local maximum occurs at a smaller momentum relative to the local minimum.
(iv). The local minimum has a value h min = h d c < 1, while the local maximum always occurs at h max > 1, leading to the emergence of all three dynamically critical phases: The doubly critical phase for h d c < h f < 1, the triply critical phase for 1 < h f < h max , and the singly critical phase for h f > h max . Also in this case, the local maximum appears at a smaller momentum relative to the local minimum.
These forms show that the triply critical phase can only emerge at h f > 1, and that the doubly critical phase only appears for h f < 1. In the case of form (ii) where only a local minimum, but not a local maximum, exists, the doubly critical phase can be detected by the sign of the leading correction to the function h f (k) at k 0, while for form (iv) the Taylor expansion around k = 0 cannot be used to detect the doubly critical phase. In short, when the sign of the leading correction to h f (k) is negative, one can be certain that form (ii) arises, whereas for a positive sign any of the other three forms can occur.
In the following, our Ansatz to uncover most of the properties of the dynamical phase diagram will be to study the function (36) around k = 0. To ease the analysis, we split h f (k) into two contributions Conveniently, the two exponents α and β enter separately into these contributions. This allows for an independent analysis of their influence on the forms (i)-(iv).  (36) can take for quenches from hi = 0. In form (i), we see that there exist only the singly critical phase, where quenches to h f > h e c = 1 generate a unique critical momentum mode that gives rise to critical times in the return rate. In form (ii), we see again the singly critical phase, but also the doubly critical phase, which occurs for quenches within the topologically nontrivial phase and allows for two critical momentum modes. In form (iii), the singly critical phase emerges for quenches just above h e c and also for large quenches. Quenches in between lead to a triply critical phase with three critical momentum modes. Finally, form (iv) exhibits the richest dynamical criticality, in that again the doubly critical phase emerges for quenches within the topologically nontrivial phase, the triply critical phase appears for quenches in a region above h e c , but for quenches above that region only the singly critical phase exists.
A. The leading long-range pairing case (β < 2) In the case β < 2 the second contribution is dominated by the nonanalytic power, leading to Consequently, only two possibilities can arise, as follows.
1. The α < (2β − 1) regime In this regime the leading-order expansion for h f (k) is given by H 1 as The coefficient Γ(1 − α) sin(πα/2)/ζ(α) is always negative and therefore the doubly critical phase is always present in this regime. The subleading term is given by H 2 and is always positive. In this regime there is no triply critical phase. Examples of this regime are shown in the top right panel of Fig. 1.

The α > (2β − 1) regime
In this case the roles of α and β are reversed. The leading-order term is given by H 2 and the expansion reads In this regime, the leading term is always positive and thus the existence of the doubly critical phase cannot be directly inferred. The α, β dependence of the subleading correction depends on the value of α.
In such a case, the triply critical phase exists. Examples for this are shown in the bottom right panel of Fig. 1, for α < β + 1 (α = 2.3 and β = 1.5) and for α > β + 1 (α = 2.8 and β = 1.5), where in both cases α > 2β − 1. However, it is also possible that form (i) or (iii) may emerge instead, in which no doubly critical phase exists as both the maximum and minimum of h f (k) occur above unity. As examples, we show in the top left panel of Fig. 1 how for α < β + 1 (α = 2.85 and β = 1.9) and for α > β + 1 (α = 3 and β = 1.9) form (i) can emerge where only a singly critical phase is possible. Additionally, we shown examples of how form (iii) may arise in the bottom left panel of Fig. 1 for α < β + 1 (α = 2.65 and β = 1.7) and for α > β + 1 (α = 2.8 and β = 1.7). Note that all these examples satisfy α > 2β − 1. When β > 2, the second contribution to h f (k) is dominated by the analytic power, which leads to two different regimes as a function of α. It is worth noting that in this region no triply critical phase can ever exist as we shall discuss in the following.

The α < 3 regime
In this case the Taylor expansion for the critical final field reads Since the coefficient is always negative, the doubly critical phase is always present in this regime. This would bring about form (ii), which excludes the emergence of a triply critical dynamical phase. An example is shown in the top right panel of Fig. 1 for α = 1.5 and β = 2.5.

The α > 3 regime
In this regime both of the nonanalytic terms are irrelevant and the leading contribution reads The existence of the doubly critical phase depends on the coefficient For c 1 < 0 the doubly critical phase exists within form (ii), meaning that the triply critical phase cannot exist in this case, and an example of this regime is given in the top right panel of Fig. 1   As we discuss in Sec. VI, of particular interest to us is the doubly critical phase. This is due to recent experimental observations 34 of so-called accidental dynamical vortices that appear in time-resolved Bloch tomography measurements for quenches within the same topological phase. Our results indicate that this is indeed possible even in models with no interactions such as the Haldanelike model realized in Ref. 34. Indeed, the generally accepted explanation -based on results of long-range interacting quantum spin models -for why the dynamical critical point can be smaller than the equilibrium critical point has been that this is most likely due to interactions. [18][19][20][21][22][23][24][25][26] Nevertheless, as indicated in our previous work Ref. 50 and solidified in our current study, the case of h d c < h e c is very prevalent for noninteracting free fermionic systems with long-range hopping and pairing. Fig. 3 shows h d c as a function of the pairing exponent β at fixed values of the hopping exponent α. For small to intermediate values of α this function is not monotonous, but instead goes from very small values at β 1, to a maximum (which is still below unity) at intermediate values of β, and then asymptotically towards a smaller value at β → ∞. Consequently, we see that even for nearestneighbor pairing the doubly critical phase can still arise for small enough α. Indeed, as already discussed in the analysis of (43), in the regime of α < 3 the doubly critical phase arises regardless of the value of β.

VI. RESULTS AND DISCUSSION
Finally, we compare the return rate and the SOP for quenches leading to the three different dynamical phases identified in Sec. V, and we discuss the intimate relation between the periodicities of cusps in the return rates and those of the zero crossings of the SOP.
It has already been shown in the nearest-neighbor case that for quenches of the chemical potential h across the equilibrium phase boundary h e c the kinks appearing in the return rate correspond to zeros of the string order parameter's norm M(t). 43 In Fig. 4 we show the time evolution of the square of the TDVP-computed longitudinal magnetization in the TFIC (25) and that of the norm of the string order parameter in the corresponding Jordan-  1) and the norm of the string order parameter of its Jordan-Wigner mapping for a quench from hi = 0 to h f = 3. The result for the magnetization is simulated using TDVP. [46][47][48] Wigner mapping, which is the LRKC Hamiltonian (1) in the limit of α, β → ∞. When α and β are finite, we can no longer in general find a spin model that can be mapped onto the corresponding LRKC Hamiltonian (1). Consequently, the LRKC can no longer be associated with a local order parameter such as the longitudinal magnetization. Nevertheless, the fact that the SOP norm and the square of the magnetization are the same observable in the nearest-neighbor limit motivates the study of the SOP -which is accessible regardless of the value of α and β -as a viable (nonlocal) order parameter. It further motivates relating its dynamics to that of the return rate, as is done in the case of the local order parameter in spin models.
The picture in the long-range case is similar to its nearest-neighbor counterpart in that large quenches across h e c lead to SOP zeros and return rate cusps that exhibit the same periodicity. This is shown in Fig. 5 where we calculate the dynamical evolution of the system after a quench of the chemical potential to h = h f = 2 starting from the ground state of the h = h i = 0 Hamiltonian for different values of α and β. The return rate and SOP both display regular singularities in the form of cusps and zero crossings, respectively. As in the nearest-neighbor case, 43 these singularities occur simultaneously. It is remarkable that the time evolution of the return rate remains very regular (and is similar to the nearest-neighbor case) in the case of β < α, regardless of the actual β value.
Having identified the doubly critical phase boundary h d c , we can study the dynamical properties of this phase by quenching the system from h i = 0 to the final field h f = 1.075h d c < h e c , where two dynamical critical momenta appear. Not surprisingly, for small α and β, numerical calculations are plagued by considerable finite-size effects. Therefore, we will focus on two particular sets of decay rates (α, β) = (1.75, 1.50) and (1.90, 1.50), both of which i.e., in the singly critical dynamical phase. The SOP is calculated for N = 4000 sites while the return rate r(t) for N = 10000 sites, where these system sizes are chosen such that the corresponding quantity reaches convergence. The second time-derivative (curvature) of the return rate,r(t), is also provided in order to highlight cusps that may not appear clearly in r(t). As can be seen, the cusps of r(t) and the zeros of M(t) exhibit the same periodicity.  . Return rate and SOP norm for a quench from hi = 0 into the the doubly critical phase with h f = 1.075h d c < h e c , for two choices of (α, β). In both panels the top quadrant shows the SOP norm and the middle quadrant the return rate. The bottom quadrant shows the curvaturer(t) of the return rate in order to identify the critical times, which are further indicated by the vertical dashed lines. The critical times exhibit two distinct periodicities, distinguished by the black and gray dashed lines, respectively, arising from the two distinct critical momenta of this quench. The system sizes were chosen as in Fig. 5.
are inside the relevant region where the effect of longrange couplings is expected to be more drastic and the doubly critical phase is larger; cf. Fig. 3. Indeed, as per the analysis of Sec. V, in both these cases β < 2 and α < 2β − 1, meaning that for quenches just below the equilibrium critical point a doubly critical phase emerges, while for quenches above it, only the singly critical phase can appear; cf. top right panel of Fig. 1. The dynamics of the return rate and the SOP for these two parameter sets are shown in Fig. 6. Even though not as clear cut as in the case of the singly critical phase, we see that the cusps in the return rate and the zero crossings of the SOP nevertheless share similar periodicities. Moreover, the number of cusps over the time interval we show equals the number of zero crossings the SOP makes. Once again, this indicates a nontrivial connection between DQPT and the dynamics of the SOP, much the same way the local order parameter in spin chains has been shown to be intimately connected to regular cusps in the return rate. 18,[20][21][22][23][24][25] We note here how both the cusps of the return rate and the zero crossings of the SOP show two periodicities due to these quenches having two emerging critical momenta. These periodicities are indicated by two different fonts of dashed lines, each corresponding to the critical times of a given periodicity.
It is interesting here to relate the doubly critical phase we observe in our simulations to recent experimental findings. In the experiment of Ref. 34, a topological Haldane-like model is implemented using spin-polarized fermions in a driven optical lattice. Initially suppressed with a large energy offset, the tunneling strength between the two sublattices is then suddenly quenched, where near-resonant driving reestablishes tunneling in the final Floquet Hamiltonian. A time-resolved Bloch state tomography is then carried out for the pseudo-spin-1/2 in momentum space arising from the sublattice degree of freedom. This system gives rise to the same physics of spin chains in the limit of no interactions. Dynamical vortices -which indicate that the time-evolved wavefunction is orthogonal to the initial state, and are therefore representative of cusps in the Loschmidt return rate -are directly observed when quenching between two phases of different Chern numbers. However, accidental dynamical vortices are also observed for quenches close to the above phase transition, but not crossing it. Indeed, the right panels of Fig. 1 show two forms of h f (k) described in (36) that lead to the dynamical phase diagram hosting a doubly critical phase in which nonanalytic behavior in the return rate arises for quenches within the topologically nontrivial phase in the LRKC. This could very well be the underlying phenomenon behind the accidental dynamical vortices in Ref. 34. One possible way of ascertaining this would be to study these accidental dynamical vortices as a function of varying the range of the hopping terms of the Floquet system into which they quench. As detailed through our analytic treatment and numerical results in Sec. V, there are certain regimes of the hopping and pairing exponents (α, β) where the doubly critical phase is prominent -such as in forms (ii) and (iv) in Fig. 1 or outright absent -such as in forms (i) and (iii) in Fig. 1. In particular, as proven in Sec. V and illustrated in Fig. 2, when α is sufficiently small (i.e., the hopping is sufficiently long-range), the dynamical critical phase will emerge regardless of the value of β, even when pairing is nearest-neighbor. Therefore, by making α smaller, more of the abovementioned accidental dynamical vortices are expected to appear deeper within the same topological phase and for significantly smaller quench distances.
However, here we have to be careful and note that there are nontrivial differences between our model and that realized in Ref. 34. Not only is their model twodimensional, but their quenches also always start in the topologically trivial phase. In contrast, for the case of the (one-dimensional) LRKC that we consider here, quenches within the topologically trivial phase do not seem to give rise to cusps in the return rate (not shown, but extensively studied by us). Nevertheless, this may be different if the LRKC is extended to two dimensions as in the experimental model. In this case, it may be possible that quenches within the topologically trivial phase, and not just its nontrivial counterpart, will also lead to nonanalytic behavior in the return rate due to higher overall connectivity. We leave this question open for possible future investigation.
Finally, we consider quenches that give rise to the triply critical dynamical phase in the LRKC. Unlike the case of the doubly critical phase, and similar to that of the singly critical phase, the triply critical phase occurs for quenches crossing h e c , albeit it seems to only appear when the quench is right above the equilibrium critical point, h f h e c , and also when the model is in the relevant region of long-range pairing β < 2; cf. Sec. V B. Note that the region of h f to which one must quench in order to observe the triply critical phase can either be immediately above the equilibrium point up to a small h f > 1 (see form (iv) of h f (k) in bottom right panel of Fig. 1), or within a small region slightly above, but separated from, h e c (see form (iii) of h f (k) in bottom left panel of Fig. 1). In Fig. 7 three different fonts of dashed lines indicate the critical times in the return rate arising due to three distinct critical momenta for these quenches. Here, the connection with the zero crossings in the SOP is more subtle than in the case of the singly and doubly critical phases. The fact that the triply critical phase appears in a very small regime of (α, β) makes the investigation of this subtlety more challenging. This is mainly because with h f so close (from above) to the equilibrium critical point, most quenches would require impractically large lattice sizes in order to obtain converged dynamics due to critical fluctuations. Nevertheless, it is still not completely surprising here that the connection is not one-to-one as in the singly critical phase. As briefly discussed in Sec. IV, when longrange interactions are present in quantum spin chains, the connection between cusps in the return rate and zeros in the order parameter becomes less straightforward compared to the nearest-neighbor case. 21,23 This has been clearly demonstrated in the LMG model, which is the fully connected transverse-field Ising model, where especially at short times the periodicity of the cusps in the return rate does not match that of the zero crossings of the order parameter, in stark contrast to the case of the TFIC. At longer times, these two periodicities nevertheless do become equal in the LMG model. In the case of LRKC, the SOP decreases exponentially fast, making it much harder to tell whether this may also be the case here.

VII. CONCLUSION
We have investigated the effect of having both longrange hopping (∝ 1/r α ) and long-range pairing (∝ 1/r β ) on the dynamics of free-fermionic models in the wake of a quench. In particular, we have thoroughly constructed the rich dynamical phase diagram of the long-range Kitaev chain using scaling arguments and extensive numerical simulations. We find three main phases: a singly critical phase known from seminal works in which the cusps of the return rate correspond to a single critical momentum mode; the doubly critical phase in which cusps in the return rate correspond to two distinct critical momenta; and the triply critical phase where cusps arise due to three distinct critical momenta. We identified regimes of (α, β) associated with each of these phases, where the succinct summary would be that small α guarantees the emergence of the doubly critical phase, while small enough β is necessary for the triply critical phase to exist. Both the singly and the triply critical phases can occur for quenches crossing the equilibrium critical point, whereas in contrast the doubly critical phase appears for quenches within the topologically nontrivial phase.
Furthermore, we have related the cusps in the Loschmidt return rate to the zero crossings of the string order parameter. Our results indicate that there is a close connection between the cusps and zero crossings, especially in the case of singly and doubly critical dynamical phases. For large quenches from the topologically nontrivial phase crossing the equilibrium critical point, a singly critical phase arises regardless of the values of α and β. This extends what is already known on DQPT as a probe of the underlying equilibrium physics, but also shows that the connection between DQPT and equilibrium physics is not one-to-one, especially in the case of the doubly critical phase as it arises for quenches within the topologically nontrivial phase, i.e., without having to cross any equilibrium critical point.
Our work shows experimental promise as it may explain the so-called accidental dynamical vortices from time-resolved Bloch tomography measurements in Ref. 34 that appear for quenches close to, but not crossing, the topological phase transition. This parallels our observation of cusps in the return rate for quenches within the same topological phase, giving rise to a doubly critical dynamical phase possessing two distinct periodicities in the return rate cusps. Our analysis suggests that by tuning the hopping range these accidental vortices can be better understood in that for longer hopping ranges more vortices will appear deeper within the same topological phase and for smaller quenches.
Finally, it is worth mentioning that our results suggest that the SOP can offer an indirect way of measuring the nonanalyticities of the Loschmidt return rate. In the work of Refs. 61 and 62, a quantum gas microscope has been used to measure the string order, respectively, in one-dimensional bosonic Mott insulators and the Fermi-Hubbard chain. Given the close connection evident in our simulations between the return rate cusps and the SOP zeros, this would be a viable and much more practical way of observing dynamical critical behavior than actually measuring the return rate itself, which becomes exponentially harder with system size. In Fig. 8 we compare the return rates of the TFIC (25) and the LRKC (1) in the limit α, β → ∞. These two models are identical via an exact Jordan-Wigner mapping.
We see that there is a qualitative difference between the return rates in the case of a quench from the ordered phase (h i < h d c ), with the periodicity of cusps in the LRKC being half that of the TFIC. This is due to the fact that in the fermionic system the ground state is always nondegenerate, while in the case of the TFIC, the ground state in the ferromagnetic phase is always doubly degenerate. The return rates can be brought to quantitative and qualitative agreement by updating the return rate definition to where {|ψ i,0 , |ψ i,1 } are the two degenerate ground states at h i < h e c . On the other hand, the return rates are in full qualitative and quantitative agreement in the case of quenches from the paramagnetic phase. This is because in the latter phase the TFIC has a nondegenerate ground state.