Universal dynamical scaling of long-range topological superconductors

TL;DR

This paper investigates the out-of-equilibrium dynamics of long-range p-wave superconducting wires, revealing how the heat generated during a phase transition scales with the quench rate and identifying regimes with universal behavior.

Contribution

It introduces a scaling law for heat production during quenches in long-range superconductors and uncovers anomalous dynamical universality in certain parameter regimes.

Findings

Heat scales with quench rate as δ^θ, with θ depending on interaction range.

Universal scaling can be linked to equilibrium critical exponents in specific regimes.

Anomalous dynamical universality occurs when hopping dominates over pairing.

Abstract

We study the out-of-equilibrium dynamics of $p$-wave superconducting quantum wires with long-range interactions, when the chemical potential is linearly ramped across the topological phase transition. We show that the heat produced after the quench scales with the quench rate $\delta$ according to the scaling law $\delta^\theta$, where the exponent $\theta$ depends on the power law exponent of the long-range interactions. We identify the parameter regimes where this scaling can be cast in terms of the universal equilibrium critical exponents and can thus be understood within the Kibble-Zurek framework. When the electron hopping decays more slowly in space than pairing, it dominates the equilibrium scaling. Surprisingly, in this regime the dynamical critical behaviour arises only from paring and, thus, exhibits anomalous dynamical universality unrelated to equilibrium scaling. The discrepancy from the expected Kibble-Zurek scenario can be traced back to the presence of multiple universal terms in the equilibrium scaling functions of long-range interacting systems close to a second order critical point.

Equations

$$\displaystyle\hat{H}=-\sum_{i}\Bigl{[}\sum_{r>0}\left(j_{r}\hat{c}^{\dagger}_{i}\hat{c}_{i+r}+\Delta_{r}\hat{c}^{\dagger}_{i}\hat{c}^{\dagger}_{i+r}+{\rm H.c.}\right)+\mu\hat{c}^{\dagger}_{i}\hat{c}_{i}\Bigr{]}+\mathcal{C},$$

$$j_r^{\alpha }$$

$$\mu ={\mu }_c-\delta \cdot t,$$

$$\hat{H}(t)=\sum _k{\hat{\psi }}_k^{†}{\hat{\mathcal{H}}}_k(t){\hat{\psi }}_k,$$

$$j_{\alpha }(k)$$

$${\Delta }_{\beta }(k)$$

$$i\frac{d}{dt}\left(\begin{array}{c}{u}_k\\ v_k\end{array}\right)=\left(\begin{array}{cc}{\epsilon }_{\alpha }(k,t)& -{\Delta }_{\beta }(k)\\ {\Delta }_{\beta }(k)& {\epsilon }_{\alpha }(k,t)\end{array}\right)\left(\begin{array}{c}{u}_k\\ v_k\end{array}\right),$$

$$p_k\simeq \exp \left(-\frac{\pi {\Delta }_{\beta }(k{)}^2}{\delta }\right),$$

$$c=\{\begin{array}{ll}{\left(\frac{\cos (\beta \pi /2)\Gamma (1-\beta )}{\zeta (\beta )}\right)}^2& \text{for}1<\beta <2,\\ {\left(\frac{\zeta (\beta -1)}{\zeta (\beta )}\right)}^2& \text{for}\beta >2.\end{array}$$

$$n_{\mathrm{exc}}=\frac{1}{N}\sum _k⟨{\gamma }_k^{†}{\gamma }_k⟩=\frac{1}{N}\sum _k{p}_k,$$

$$\underset{\delta \to 0}{\lim }{n}_{\mathrm{exc}}\sim {\delta }^{\theta }\mathrm{with}\theta =\{\begin{array}{ll}(2\beta -2{)}^{-1}& \text{for }\beta \le 2,\\ 1/2& \text{for}\beta >2.\end{array}$$

$$p_k=1-|u_k^f{u}_k^{\ast }(t_f)+v_k^f{v}_k^{\ast }(t_f){|}^2,$$

$${\alpha }_k\approx \int ⟨k|{\hat{\partial }}_{\mu }|0⟩e^{\frac{i}{\delta }{\int }^{\mu }\left(E_k({\mu }^{\prime })-E_0({\mu }^{\prime })\right)d{\mu }^{\prime }}d\mu ,$$

$$E_k(\mu )-E_0(\mu )$$

$$⟨k|{\hat{\partial }}_{\mu }|0⟩$$

$$H=-\sum _i\sum _{r>0}\left(j_r{c}_i^{†}{c}_{i+r}+{\Delta }_r{c}_i^{†}{c}_{i+r}^{†}+H.c.\right)-\mu \sum _i\left(c_i^{†}{c}_i-\frac{1}{2}\right),$$

$$j_r^{\alpha }$$

$${\Delta }_r^{\beta }$$

$$N_{\gamma }=2\sum _{r=1}^{N/2}\frac{1}{r^{\gamma }}\to 2\zeta (\gamma ),$$

$$c_k=\frac{1}{\sqrt{N}}{e}^{-i\frac{\pi }{4}}\sum _{j\in \mathbb{Z}}{c}_j{e}^{ikj}$$

$$H=\sum _k\left[(c_k^{†}{c}_k-c_{-k}{c}_{-k}^{†}){\epsilon }_{\alpha }(k)+(c_k^{†}{c}_{-k}^{†}+c_{-k}{c}_k){\Delta }_{\beta }(k)\right],$$

$${\epsilon }_{\alpha }(k)$$

$$j_{\alpha }(k)$$

$${\Delta }_{\beta }(k)$$

$$c_k=u_k{\gamma }_k+v_{-k}^{\ast }{\gamma }_{-k}^{†},$$

$$H=2\sum _k{\omega }_k\left({\gamma }_k^{†}{\gamma }_k-\frac{1}{2}\right),$$

$${\omega }_k=\sqrt{{\epsilon }_{\alpha }(k{)}^2+{\Delta }_{\beta }(k{)}^2}.$$

$$(u_k,v_k)=\left(\cos \frac{{\theta }_k}{2},\sin \frac{{\theta }_k}{2}\right),$$

$$\tan {\theta }_k=\frac{{\Delta }_{\beta }(k)}{{\epsilon }_{\alpha }(k)}.$$

$$j_{\alpha }(k)/J$$

$${\Delta }_{\beta }(k)/d$$

Full text

One of the major challenges of contemporary physics is the identification of quantum phases of matter, which can serve as platforms for quantum computers. In this perspective, topological superconductors Hasan2010 ; Bernevig2013 are promising constituents for quantum devices Nayak2008 ; Terhal2015 ; Kraus2013a ; Mazza2013 , thanks to the presence of gapless Majorana modes, the so-called Majorana zero modes (MZM), which are localised at the chain edges and topologically protected. Since the first theoretical evidence of MZMs in superconducting wires Kitaev2001 , several experimental platforms have revealed consistent signatures of Majorana physics both in one-dimensional Mourik2012 ; Deng2012 ; Das2012 ; Albrecht2016 and two-dimensional Wang2012 ; He2014 ; Sun2016 ; He2017 geometries. More recently, models of $p$-wave superconducting wires with long-range (LR) deformations have shown more robust topological properties Viyuela2015 ; Viyuela:2018fpv , while strong enough LR pairing effects alter the nature of the topological phase Vodola2014 ; Lepori2015 ; Lepori2017 ; Lepori2017add ; Alecce2017 and the spreading of correlations Foss-Feig2015 ; Cevolani2015 ; Vodola2016 . Experimental realisations of LR topological superconductors employ one-dimensional arrays of magnetic impurities on top of a conventional superconducting substrate Nadj-Perge2014 ; Pawlak2016 ; Ruby2017 , leading to the realisation of an effective Kitaev Hamiltonian with both LR pairing and LR hopping Pientka2013 ; Klinovaja2013 ; Pientka2014 ; Neupert2016 . In this context, understanding slow variations (quenches) of control fields in quantum systems is fundamental to adiabatic protocols  Nielsen2000 , since Majorana excitations cannot be realised by sudden manipulations of the system Perfetto2013 . These investigations constitute a fundamental contribution towards the understanding of dynamical scaling for quenches across topological phase transitions Ueda2010 .

In this Letter we characterise the out-of-equilibrium dynamics of a $p$-wave superconducting quantum wire with long-range interactions, whose chemical potential is linearly ramped across the equilibrium critical point. We determine the density of defects produced after the ramp and show that it scales as a power of the quench rate. We then connect the power-law exponent with the equilibrium critical properties and topological features and determine the phase diagram for the dynamics as a function of the decay exponent of the hopping and pairing terms.

We consider spinless electrons hopping across the $N$ sites of a linear chain in presence of $p$-wave pairing. The Hamiltonian reads

[TABLE]

where operators $c^{\dagger}_{i}$ create a fermion at site $i$ and fulfill the anticommutation relations $\{c_{i},c_{j}^{\dagger}\}=\delta_{ij}$. Here, $\mu$ denotes the chemical potential, $\mathcal{C}=N\mu/2$ is an energy offset, $j_{r}$ and $\Delta_{r}$ are the hopping and pairing amplitudes, respectively, and depend on the intersite distance $r$ according to the power laws (reported here for open boundary conditions):

[TABLE]

with the hopping exponent $\alpha>1$, the pairing exponent $\beta>1$, the coefficients $J,d>0$, and $N_{\gamma}=2\sum_{r=1}^{N/2}r^{-\gamma}$ the Kac scaling, which guarantees extensivity of the energy Campa2014 . For sufficiently fast decaying interaction and hopping terms the system possesses two different phases separated by the quantum critical point $\mu_{c}=2J$ Kitaev2001 . In the thermodynamic limit the two topological phases can be distinguished by the bulk topological invariant $w$: For $|\mu|>\mu_{c}$ the ground state is nondegenerate and $w=0$; in the nontrivial phase $|\mu|<\mu_{c}$ the bulk topological invariant $w=1$, and the ground state is doubly degenerate and can host MZMs, see Fig. 1. At finite size $N$ the spectrum is always gapped and for open boundary conditions the MZMs remain localized at the edges of the chain. The presence of the LR pairing and hopping terms in Eq. (21) does not alter this phase diagram nor the values of the bulk topological invariant as long as $\alpha,\beta>1$ Vodola2014 ; Viyuela:2018fpv ; Alecce2017 . Nonetheless, LR connectivity modifies the universal critical behaviour of the model by changing the critical exponents. The resulting phases are displayed in Fig. 2. Note that the equilibrium phase diagram of the long-range Kitaev chain radically differs from the one of the long-range quantum Ising model Defenu2016 ; Defenu2017a .

In the following, we analyse the dynamics during slow variations of the chemical potential across the critical value according to

[TABLE]

where time varies in the interval $[-\mu_{c}/\delta,\mu_{c}/\delta]$, i.e., from the topologically trivial phase $\mu=2\mu_{c}$ deep into the nontrivial phase at $\mu=0$. We note that the time-dependent dynamics have been solved for an Ising model in transverse field Kolodrubetz2012 , which can be mapped to the Kitaev model for $\alpha,\beta\to\infty$ Fradkin1989 . Below we derive an exact solution which is valid for general $\alpha,\beta>1$ and in the thermodynamic limit. This solution allows us to determine the thermodynamic functions after the ramp. For this purpose we rewrite the Hamiltonian (21) using momentum-space operators $\hat{c}_{k}=e^{i\pi/4}\sum_{r\in\mathbb{Z}}c_{r}e^{ikr}/\sqrt{N}$ with $k\in[-\pi,\pi)$. Using the spinor representation $\hat{\psi}_{k}=(\hat{c}_{k},\hat{c}_{-k}^{\dagger})^{T}$, the Hamiltonian is the sum of $2\times 2$ block matrices $\mathcal{H}_{k}\equiv\boldsymbol{h}_{k}(t)\cdot\boldsymbol{\hat{\sigma}}$,

[TABLE]

where $\boldsymbol{\hat{\sigma}}$ is the vector of the Pauli matrices $\hat{\sigma}^{j=1,2,3}$ and $\boldsymbol{h}_{k}(t)=(\Delta_{\beta}(k),0,\varepsilon_{\alpha}(k,t))$ is the pseudo-spin vector. Its elements depend on $\varepsilon_{\alpha}(k,t)=\mu(t)/2-j_{\alpha}(k)$ and on the momentum-space hopping and pairing coefficients:

[TABLE]

where $\operatorname{Li}_{\alpha}(z)$ denotes the polylogarithm and $\zeta(\alpha)$ the Riemann zeta function Abramowitz1964 . The Hamiltonian is diagonalized in terms of the fermionic quasiparticle operators $\hat{\gamma}_{k}(t)=u_{k}(t)\hat{c}_{k}+v_{-k}^{*}(t)\hat{c}_{-k}^{\dagger}$ to obtain \hat{H}=\sum_{k}\omega_{k}(t)\bigl{(}\hat{\gamma}^{\dagger}_{k}(t)\hat{\gamma}_{k}(t)-\frac{1}{2}\bigr{)} with the quasiparticle spectrum $\omega_{k}(t)=2\sqrt{\varepsilon_{\alpha}(k,t)^{2}+\Delta_{\beta}(k)^{2}}$. The pseudo-spin vector $\boldsymbol{h}_{k}(t)$ identifies a direction in the two-dimensional plane of the Hamiltonian space. At a given instant of time, integrating the angle $\theta_{k}={\rm atan}(h_{k}^{1}/h^{3}_{k})={\rm atan}(\Delta_{\beta}(k)/\varepsilon_{\alpha}(k))$ over the Brillouin zone yields the bulk topological invariant $w=\oint d\theta_{k}/(2\pi)$ of the corresponding equilibrium phase.

The dynamics of the Kitaev chain can be exactly described by the Heisenberg equations of motion for the original creation and annihilation operators, ${\rm i}d\hat{c}_{k}/dt=[\hat{c}_{k},\hat{H}]$. These equations can be cast into a matrix evolution for the Bogolyubov coefficients,

[TABLE]

which can be mapped into the Landau-Zener form dziarmaga2005 ; dziarmaga2010 ; Dutta2017 , see App. C. The excitation probability $p_{k}(t)$ can be computed exactly Damski2005a ; Bialonczyk:2018sbn , see App. B. For a slow quench to the final time $\tau_{f}=\mu_{c}/\delta$, the excitation probability is well approximated by the Landau-Zener formula

[TABLE]

which becomes exact for $k\ll\frac{\pi}{2}$ in the slow ramp limit $\delta\to 0$. From Eq. (8) we find population inversion $p_{k}\geq 1/2$ when $|k|<k_{\rm th}$. The threshold value $k_{\rm th}$ is determined analytically from a low-momentum expansion and reads $k_{\rm th}=[\delta\log(2)/(\pi c)]^{\theta}$ with

[TABLE]

It is remarkable that $k_{\mathrm{th}}\to\infty$ for $\beta\to 1^{+}$, corresponding to negative temperatures. In the small $\delta$ limit, this effect is only visible very close to the singular limit $\beta\gtrsim 1$, while for intermediate $\beta$’s the tendency is reversed, see Fig. 3. These results have been numerically verified taking the full $k$ dependence of $\Delta_{\beta}(k)$ into account. We note that stable athermal distributions are generally expected in systems with diverging long-range interactions Kastner2011 , and the case we analyse here seems to be no exception.

Using Eq. (8) we can derive several thermodynamic properties for asymptotically slow drive $\delta\to 0$. We discuss here the excitation density zurek2005 ; chandran2012 ; dziarmaga2010 ; polkovnikov2005 ,

[TABLE]

which we compute in the thermodynamic limit, thus replacing the sum by an integral over the interval $k\in[-\pi,\pi)$. In the $\delta\to 0$ limit the exponent of $p_{k}$ in Eq. (8) diverges and the total contribution to the integral only comes from the saddle point. Expanding around the vanishing effective frequency, we find the scaling law

[TABLE]

At the border $\beta=2$ we find $n_{\rm ex}\propto\sqrt{\delta}/\log\delta$. These scalings are valid irrespectively of the value of $\alpha$, since the defect density solely depends on $\beta$. The corresponding dynamical phase diagram is depicted in Fig. 2. Remarkably, the dynamical phases for $\alpha>2$ correspond to the regions of the equilibrium phase diagram, but in the hopping dominated regime $\alpha<2$ a different universal dynamical scaling arises. Such universal dynamical scaling with $\beta$ cannot be related to the equilibrium critical exponents, which involve $\alpha$, as generally happens in Fermi systems dziarmaga2010 ; polkovnikov2005 ; de_grandi2010 and its appearance can be traced back to the violation of the equilibrium scaling hypotheses due to the LR nature of the interactions.

In order to verify our analytical prediction, we numerically integrated Eq. (7). Initially at $t=t_{i}=-\mu_{c}/\delta$ the system is at equilibrium with Bogolyubov coefficients $u^{i}_{k},v^{i}_{k}=\cos\frac{\theta_{i}(k)}{2},\sin\frac{\theta_{i}(k)}{2}$ where $\tan\theta_{i}(k)=\frac{\Delta_{\beta}(k)}{\varepsilon_{\alpha}(k)}$ and $\mu(t_{i})=2\mu_{c}=4J$. The Bogolyubov coefficients are then evolved numerically according to Eq. (7) for a grid of $k$ points in the interval $[-\pi,\pi)$. Due to non-adiabatic effects arising during the critical stage of the dynamics $t\simeq 0$, the resulting amplitudes at the final time $t_{f}=\mu_{c}/\delta$ differ from the ones of the equilibrium Hamiltonian with $\mu(t_{f})=0$. In order to quantify these deviations, we consider the excitation probability of each state $k$, after the slow ramp, with respect to its equilibrium ground state,

[TABLE]

where $(u^{f}_{k},v^{f}_{k})$ are the equilibrium Bogolyubov amplitudes at $\mu=0$, while $(u_{k}(t_{f}),v_{k}(t_{f}))$ are the ones at the end of the dynamical evolution.

The excitation probability at the end of the slow quench, calculated according to Eq. (12), is shown in Fig. 3 as a function of the momentum $k$ for $\delta=0.5,\,0.05$ in panels (a) and (b), respectively. As the dynamical protocol crosses the $\mu_{c}=2J$ critical point, only low momentum modes $k\approx 0$ become soft during the dynamics. Indeed, Eq. (8) only applies to excitations modes with $k<\pi/2$, as follows from the Landau-Zener mapping, see App. C, while high energy modes $k>\pi/2$ remain adiabatic and their excitation probability is not shown. Numerical points for the Bogolyubov modes excitation probability for $\alpha=(\infty,1.75,1.50,1.25)$ are shown by squares, crosses, circles and triangles respectively (see legend of panel b), while the different values of $\beta=(\infty,1.75,1.5,1.25)$ are reported respectively from top to bottom (gray, green, blue and red). In Fig. 3 we observe almost perfect agreement with the predictions of Eq. (8) in the slow ramp case $\delta=0.05$. Indeed, corrections from finite and slightly asymmetric endpoints $t_{i}$ and $t_{f}$ do not influence the universal behaviour obtained in the $\delta\to 0$ limit at small momenta and can be safely discarded, see App. D.

The result of Eq. (11) contradicts the result found using adiabatic perturbation theory, which produces the Kibble-Zurek relation between the universal slow dynamics and the equilibrium critical exponents $\theta=\nu/(1+z\nu)$ dziarmaga2010 ; polkovnikov2005 ; de_grandi2010 . Here, the critical exponents $z\nu$ and $z$ describe the scaling of the spectrum at the critical point, $\omega_{k=0}\propto|\mu-\mu_{c}|^{z\nu}$ and $\omega_{k\to 0}\propto k^{z}$. In particular, at lowest order in the adiabatic expansion, the excitation probability $p_{k}=|\alpha_{k}|^{2}$ of the Bogolyubov quasi-particle states $|k\rangle=\hat{\gamma}_{k}^{\dagger}|0\rangle$ are given by the squared transition amplitudes induced by the perturbation operator $\hat{\partial}_{\mu}=\partial\hat{H}(\mu)/\partial\mu$ over the Bogolyubov vacuum $|0\rangle$ integrated over the whole dynamical trajectory

[TABLE]

where $E_{k}(\mu)$ is the energy of the state $|k\rangle$. In the $\delta\to 0$ limit, the saddle-point approximation holds and the integral only receives contribution from the vanishing gap region of the trajectory, i.e. the critical point. There, one can employ the universal scaling relations dziarmaga2010 ; polkovnikov2005 ; de_grandi2010

[TABLE]

where $\Delta$ is the minimal gap $\Delta\propto|\mu-\mu_{c}|^{z\nu}$. Inserting Eqs. (14) and (15) into Eq. (13) and making the integration dimensionless, one finds the universal scaling variables $\eta=k\delta^{-\frac{\nu}{1+z\nu}}$ and $\zeta=k^{1/\nu}(\mu-\mu_{c})$. Rephrasing the adiabatic perturbation theory expression for the defect density $n_{\mathrm{exc}}\approx\int dk|\alpha_{k}|^{2}/(2\pi)$ in terms of the universal variables $\eta$ and $\zeta$ immediately leads to the Kibble-Zurek result $\theta=\nu/(1+z\nu)$, see Eq. (13) and Refs. dziarmaga2010 ; polkovnikov2005 ; de_grandi2010 . Since for the p-wave superconducting Hamiltonian in Eq. (21) one has $z\nu=1$ and $z=\phi-1$, where $\phi=\mathrm{min}(\alpha,\beta)$, we can conclude that the scaling exponent $\theta$ in Eq. (11) is inconsistent with the Kibble-Zurek scaling in the region $\alpha<\beta$.

We refer to this unexpected behaviour as anomalous universal scaling. We report its extent in Fig. 3 where the example case of LR hopping $\alpha=1.25$ and short range pairing $\beta=\infty$, well inside the anomalous universal scaling region, is studied. The excitation probability is reported as a function of the universal scaling variable $\eta=k\delta^{-\frac{\nu}{1+z\nu}}$ for several $\delta$ values. Remarkably, for this scaling the curves do not collapse, see Fig. 3. Instead, universality is recovered when one considers the proper dynamical exponent $z_{d}=\nu_{d}^{-1}=1$, for nearest neighbour pairing which is the only responsible for the dynamics Fig. 3. Indeed, perfect collapse of the excitation probabilities for various $\delta$ is observed in terms of the correct scaling variable $k/\delta^{\nu_{d}/(1+z_{d}\nu_{d})}=k/\delta^{1/2}$, see the inset in Fig. 3.

In conclusion, we have demonstrated that long-range coupling terms can lead to a novel scaling behaviour of heat produced by slow quenches in critical topological superconductors. We have shown that this behaviour cannot be understood within the framework of the Kibble-Zurek scaling. In particular, the introduction of long-range hopping terms which decay slower than the pairing couplings modifies the equilibrium critical properties but not the dynamical critical exponent $\theta$. In the traditional case, for dominant pairing $\alpha>\beta$, the MZMs are gapped in the broken phase and, approaching the critical point, the gap closes and couples them into a single Dirac mode at $\mu=\mu_{c}$. The Dirac mode arises at the topological phase transition and dominates the low energy spectrum of the system, being also responsible for the universal slow dynamics Vodola2014 ; Dutta2017 .

Instead, for dominant hopping term $\alpha<\beta$, the critical Dirac mode is not relevant in the low energy spectrum, since the pairing term is not the leading operator in the zero momentum limit, and the equilibrium low energy theory at the critical point does not show any trace of the topological order found in the broken phase ($w=1$). However, the subleading pairing term turns out to be dangerously irrelevant and signatures of the topological order are found in the dynamics, which is always governed by the sub-leading pairing term, which is responsible for the topological transition. It is worth noting that the discrepancy between the traditional scaling argument and the anomalous universal scaling is not related to the inapplicability of the adiabatic perturbation theory expression (13), as it may occur in Bose systems due to diverging occupations Polkovnikov2008 ; Bachmann2017 ; Defenu2018 . Rather, the anomalous universal scaling is the consequence of deviations from the universal scaling hypotheses, see Eq. (14), occurring in LR systems. Similar deviations were already noticed in LR classical systems Flores-Sola2015 ; Flores-Sola2016 , but their consequences appear to be much more striking in the dynamics of quantum systems.

Our results can be straightforwardly generalised to higher dimensional cases. Moreover, we expect them to be generally valid for most of the interacting $p$-wave Hamiltonians Sau2010 ; Jason2010 ; Lutchyn2010 , which reduce to the quadratic form of Eq. (21) in the Bogolyubov approximation. These investigations are of fundamental importance in current technological applications, since slow dynamical manipulations of the Hamiltonian are necessary to realise MZMs Perfetto2013 . Finally, due to the possibility of experimentally measuring both the equilibrium and the dynamical critical scaling, the anomalous universal scaling can be used as a diagnostic for the existence of long-range tails in the hopping matrix and of topological excitations in superconducting systems.

Acknowledgements. ND is grateful to S. Ruffo, A. Trombettoni and G. Gori for useful discussions at the early stages of this work. ND and TE acknowledge financial support by Deutsche Forschungsgemeinschaft (DFG) via Collaborative Research Centre SFB 1225 (ISOQUANT) and under Germany’s Excellence Strategy EXC-2181/1-390900948 (Heidelberg STRUCTURES Excellence Cluster). LD acknowledges financial support from the BIRD2016 project of the University of Padova. GM is grateful for financial support by the DFG Priority Program no. 1929 GiRyd and by the German Ministry of Education and Research (BMBF) via the Quantera project “NAQUAS”. Project NAQUAS has received funding from the QuantERA ERA-NET Cofund in Quantum Technologies implemented within the European Union’s Horizon 2020 Programme.

Appendix A Bogolyubov Transformation

The Bogolyubov transformation which diagonalises the Kitaev hamiltonian is described here in details. Our starting point is the real space Hamiltonian for $N$ spinless fermions, Eq. (21) of the main text, which reads

[TABLE]

where the $c^{\dagger}_{i}$s are Fermionic creation operators which fulfil the anticommutation relations $\{c_{j},c_{\ell}\}=\delta_{j\ell}$. We consider power law couplings for the hopping and pairing terms:

[TABLE]

where $\bar{r}=\min(r,N-r)$ and we have considered periodic boundary conditions. The exponents of the power laws can be different and take values $\alpha>1$ and $\beta>1$, which warrant a well defined ferromagnetic state energy. The normalisation coefficients $N_{\gamma}$ ($\gamma=\alpha,\beta$) garantee that the energy is extensive. They read

[TABLE]

where the expression on the right is exact in the thermodynamic limit and $\zeta(\gamma)$ is the Riemann $\zeta$-function Abramowitz1964 .

Hamiltonian (21) is quadratic and can be explicitly integrated in momentum space via a Bogolyubov transformation Vodola2014 . For this purpose we introduce the Fourier Space transformations of operators $c_{j}$,

[TABLE]

where $k\in[-\pi,\pi)$ and it takes continous values in the thermodynamic limit. Using the Fourier representation (20) in the Kitaev hamiltonian one finds

[TABLE]

where the coefficients are a function of $k$ and read:

[TABLE]

It is convenient to employ a Bogolyubov transformation in order to diagonalise the static Hamiltonian. We choose

[TABLE]

where $u_{k}$, $v_{k}$ are the Bogolyubov coefficients and $\gamma_{k}$ satisfy the fermionic anticommutation relations $\{\gamma_{k},\gamma_{k}^{\prime\dagger}\}=\delta_{k,k^{\prime}}$. With this transformation the Hamiltonian takes the diagonal form

[TABLE]

where the eigenfrequencies read

[TABLE]

The diagonal form is found with the Bogolyubov coefficients

[TABLE]

such that

[TABLE]

This is the solution of the equilibrium model.

A.1 Taylor Expansion of the Polylogarithm

At lowest order in $k$ (namely, for $|k|\ll\pi$) , we expand the $k$-dependent coefficients and obtain the expressions:

[TABLE]

and

[TABLE]

These expressions are valid for all exponents $\alpha>1$, once the analytic continuation of the Reimann $\zeta$-function is considered for the cases $\alpha\leq 2$ and $\beta<2$. For $\beta>2$ and $\alpha>3$ the non analytic terms in Eqs. (27) and (30) become sub-leading with respect to further analytic corrections and they can safely be discarded. Now we have all the necessary information to derive a full phase diagram for the extended Kitaev chain.

Appendix B The scaling of the defect density

According to the solution of the effective LZ problem the excitation probability of each low momentum mode for an infinitely slow ramp is

[TABLE]

and the defect density can be computed integrating the excitation probability along $k$

[TABLE]

in the infinitely slow ramp limit $\delta\to 0$ the above integral has to be computed using the saddle point method. Indeed, the integral remains not negligible only on an infinitesimal neighborhood of the saddle point $k=0$, where the pairing term $\Delta_{\beta}(k)$ vanishes. According to the low momentum expansions reported in the above section for $\beta>2$ one has

[TABLE]

as it shall be for a short-range system. In the long-range regime $\beta<2$ the saddle point approximation is less straightforward due to the divergence of the Hessian in the exponent. Considering the low momentum expansion in this regime the integration reads

[TABLE]

It is convenient to define $\theta=2(\beta-1)^{-1}$ and $c=\pi\left(\cos\left(\frac{\beta\pi}{2}\right)\frac{\Gamma(1-\beta)}{\zeta(\beta)}\right)^{2}$, then we shall consider the transformation

[TABLE]

the integral then reduces to

[TABLE]

as argued in the main text.

At $\beta=2$ the low energy behavior for $\Delta_{\beta}(k)$ acquires logarithmic corrections and it shall then be treated separately. The low momentum limit in this case reads

[TABLE]

leading to the excitation probability

[TABLE]

In this case one shall introduce a more complicated transformation

[TABLE]

where $W(s)$ is the Lambert function. The integral has now been reduced to

[TABLE]

The integration boundary has to be treated with care since the transformation is not univocal. However, one shall consider that we are interested only in a small neighbourhood of $k=0$, where the expansion in Eq. (40) is valid. In this regime, it is sufficient to consider the lower branch of the Lambert function $W_{-1}(s)$, which is real in the interval $s\in[0,-1/e]$, leading to the momentum interval $k\in[0,1/e]$. In the $s\to 0^{-}$ limit $W_{-1}(s)$ obeys the asymptotic expansion Corless1996

[TABLE]

Therefore our integral can be finally approximated with

[TABLE]

the reduction of the integral boundaries to $k\in[0,1/e]$ is valid for $\delta\ll\frac{1}{e^{2}}$ and becomes exact in the $\delta\to 0$ limit. In order to proceed further, it is convenient to introduce the limit representation of the logarithm

[TABLE]

which in turns leads to

[TABLE]

Once latter expression is plugged into the integral one obtains

[TABLE]

where we once again deformed the integration range, since $s\gg\sqrt{\delta}$ contributions to the integral vanish exponentially fast in the $\delta\to 0$ limit. The summation in Eq. (50) has to be considered in the $h\to 0$ limit, where the power law contributions $\delta^{hn/2}$ become all relevant, while the $\frac{\Gamma\left(\frac{3+h\,n}{2}\right)}{1+hn}$ terms can be safely approximated as $\Gamma\left(\frac{3}{2}\right)=\sqrt{\pi}/2$ yielding

[TABLE]

As expected the case $\beta=2$ is exactly in between the pure short-range case $n_{\mathrm{exc}}\propto\sqrt{\delta}$ and the weak long range case $\beta=2-\varepsilon$, where one has $n_{\mathrm{exc}}\propto\delta^{\frac{1+\varepsilon}{2}}$ for $\varepsilon\ll 1$ and, therefore the excitation probability decays faster than in the pure short-range case. The $\varepsilon\to 0$, i.e. $\beta\to 2$, limit is placed exaclty in between, with the density of excitation decaying only logarithmically faster than in the short-range case. We have numerically verified that the introduction of sub-leading $k^{2}$ terms in the expansion of $\Delta_{2}(k)$ as well as the extension of the integration range beyond the region of validity of the low momentum expansion in Eq. (40) do not modify the scaling regime in the $\delta\to 0$ limit.

Appendix C Landau-Zener Problem

As discussed in the main text, one can employ the substitution in Eq. (12) into the dynamical evolution Eq. (11) for the Bogolyubov coefficients. The resulting dynamics takes the celebrated Landau-Zener form

[TABLE]

The Landau-Zener (LZ) evolution, Eq. (52), can be solved exactly using several approaches Zener1932 ; Wittig2005 ; damski2005 . However, this exact solution is rather cumbersome, since it is obtained in terms of Weber functions. Therefore, we will rather rely on an approximate solution, which is capable to correctly reproduce the defect scaling, only sacrificing the exactness of numerical coefficients, unimportant to our scopes. The LZ Hamiltonian can be conveniently written using Pauli’s matrices

[TABLE]

where

[TABLE]

This Hamiltonian is diagonalised analogously to the Bogolyubov transformation described in the first section, one has to introduce the angle $\theta=\mathrm{arctan}\left(\frac{1}{\Omega\tau}\right)$, useful to describe the eigenstates

[TABLE]

with the eigen-energies $\omega_{\pm}=\pm\sqrt{(\Omega\tau)^{2}+1}$.

Appendix D The finite ramp case

At the initial stage of the dynamics the system is exactly in the ground state $|\psi_{\tau=0}\rangle\equiv|-\rangle$, while at every finite time the state is given by the superposition $|\psi_{t}\rangle=\alpha_{-}(\tau)|-\rangle+\alpha_{+}(\tau)|+\rangle$, with $\alpha_{+}(\tau)^{2}+\alpha_{-}(\tau)^{2}=1$. According to adiabatic perturbation theory de_grandi2010 , the excitation amplitude at first order in $\Omega$ is given by the formula

[TABLE]

where

[TABLE]

The overlap element between the adiabatic states can be computed exactly for the LZ model,

[TABLE]

yielding the transition amplitude

[TABLE]

where, without loss of generality, we imposed $\tau_{i}=-\infty$. One can explicitly integrate the phase factor

[TABLE]

It is convenient to rescale the integration variable in Eq. (59) according to $x=\Omega\tau$,

[TABLE]

we are interested in the $\tau_{f}\gg 0$ limit of the latter we shall then separate the integral into the two contributions

[TABLE]

The above expression proves that the finite ramp dynamics is always equivalent to an infinite ramp from $\tau_{i}=-\infty$ to $\tau_{f}=+\infty$ plus a correction, which is equivalent to the exctitation amplitude of a finite ramp not crossing the critical point. The phase factor $g(x)$ has no stationary points on the real line, but it possesses an inflection point at $x=0$. Therefore, the first contribution to Eq. (62) needs to be treated separately. This computation has been already carried on in details in Ref. de_grandi2010 , yielding

[TABLE]

where the numerical coefficient $\pi/3\approx 1.05$ is surprisingly close to the exact value $1$. The second contribution can be transformed into

[TABLE]

Since $\tau_{f}$ is positive $-\tau_{f}$ is negative and the formula describes the excitation amplitude of a ramp ending below the critical point. Therefore, the integration in Eq. (64) does not contain the higher order stationary point $x=0$ and it can be safely integrated using the standard procedure for fast oscillating integrals Dingle1975

[TABLE]

Coming back to the Kitaev chain problem one has $\Omega\equiv\delta/\Delta(k)^{2}$ and $\tau_{f}=(j_{\alpha}(k)-g_{f})\Delta_{\beta}(k)/\delta$. The excitation probability for a single momentum state $k$ is

[TABLE]

Therefore, the defect density for the $p$-wave superconducting Hamiltonian in Eq. (21) after a quench starting at $g_{i}=+\infty$ and ending at $g_{f}>1$ without crossing any critical points is given by

[TABLE]

where, as long as $|g_{f}|>1$ the integral remains always convergent. In the $\delta\to 0$ limit such contribution decays quadratically with $\delta$, in agreement with our expectations for adiabatic dynamics. The latter result proves that finite ramp corrections are always negligible with respect to the non analytic scaling in the defect density generated by the low momenta during the full ramp dynamics.