Hysteresis from nonlinear dynamics of Majorana modes in topological Josephson junctions
HHysteresis from Nonlinear Dynamics of Majorana Modes in Topological Josephson Junctions
Jia-Jin Feng, ∗ Zhao Huang, ∗ Zhi Wang,
1, 3, † and Qian Niu School of Physics, Sun Yat-sen University, Guangzhou 510275, China Texas Center for Superconductivity, University of Houston, Houston, Texas 77204, USA Department of Physics, The University of Texas at Austin, Austin, Texas 78712, USA (Dated: November 5, 2018)We reveal that topological Josephson junctions provide a natural platform for the interplay between theJosephson effect and the Landau-Zener effect through a two-level system formed by coupled Majorana modes.We build a quantum resistively shunted (RSJ) junction model by modifying the standard textbook RSJ modelto take account the two-level system from the Majorana modes at the junction. We show that the dynamics ofthe two-level system is governed by a nonlinear Schr¨odinger equation and solve the equations analytically viaa mapping to a classical dynamical problem. This nonlinear dynamics leads to hysteresis in the I-V charac-teristics, which can give a quantitative explanation to recent experiments. We also predict coexistence of twointerference patterns with periods h / e and h / e in topological superconducting quantum interference devices. PACS numbers: 74.50.+r, 03.65.Sq, 85.25.Dq, 74.78.Na
I. INTRODUCTION
The topologically protected degeneracy related to nonlocalnature of Majorana modes is among the core features of topo-logical superconductors . This degeneracy is the founda-tion of fascinating topological qubits and also related tosupersymmetry in condensed matter systems . The situ-ation is interesting as well when the degeneracy is split bycouplings between Majorana modes . In particular for theone-dimensional case , the split energy levels form a typ-ical two-level system since other excitation levels have muchhigher energy .The two-level systems with their energy difference in con-trol have proved extraordinarily fertile for interesting quan-tum phenomena . By coupling two Majorana modes witha Josephson junction as in Fig. 1a, two levels with ener-gies E ∝ ± cos θ / θ the Josephson phaseand the plus/minus signs correspond to states with oppositefermion number parity. Either level can coherently transportone electron through the junction, leading to the fractionalJosephson effect I ∝ ± sin θ / . In realistic systems wherethe two levels are inevitably coupled, the two-level system hasavoided level crossings at θ = ( n + ) π as in Fig. 1b. Energyspectra with such avoided crossings are well known for the ex-istence of the Landau-Zener (LZ) transitions : the two-levelsystem enters a superposition state when the phase differenceis driven by a finite voltage drop across the junction . Thetopological Josephson junction thus hosts a natural platformfor the interplay between the LZ effect and Josephson effect .Since LZ effect has proved its impact on qualitatively chang-ing the dynamics in various systems , novel phenomenastemming from this interplay are expected on the topologicaljunctions.In this work, we study a realistic topological Josephsonjunction as sketched in Fig. 1a, where the supercurrent iscontributed by tunneling in the form of both the single elec-tron and Cooper pair. For a junction with negligible capaci-tance, we build a quantum resistively shunted junction (QRSJ)model by including the two quantum levels into the standard R (a) I (b) (c) E M E J I MajoranaCooper pair RI γ L γ R FIG. 1. (Color online) (a) Schematic of a topological Josephson junc-tion with resistance R driven by an injected current I . The single-electron tunneling through the Majorana modes γ L and γ R , and theCooper-pair tunneling induce Josephson couplings are quantified byenergy scales of E M and E J respectively. (b) Energies of the two-levelsystem defined by the two Majorana modes, with δ coming from thecoupling between γ L , R and the other two Majorana modes at the endsof the wire. The Landau-Zener transition happens at the avoidedenergy crossing with P the transition possibility. (c) Schematic ofequivalent electric circuit for topological Josephson junction. RSJ model. Under current injection, the two-level system canpass the avoided crossing again and again. At each passageit experiences a LZ transition at the near diabatic limit. Theaccumulation of multiple LZ transitions, which couple withthe nonlinear dynamical of Josephson phase, induces a noveldamped quantum oscillation. We cast the quantum model intoa classical model to solve this nontrivial dynamics by exploit-ing the method of averaging, and find that the LZ transitionsare effectively described by a nonlinear Schr¨odinger equation.We use phase-space portrait and the Poincar´e map to analyzethis nonlinear LZ effect, and reveal a separatrix which cate- a r X i v : . [ c ond - m a t . s up r- c on ] N ov gorizes the dynamics into two distinct oscillatory behaviors.Within the separatrix, we obtain an analytically solution forthe damped quantum oscillation, which agrees well with nu-merical simulations. We further show that this damped os-cillation leads to hysteresis in the I-V curves, which gives aquantitative explanation to the recently reported “unexpected”hysteresis in HgTe topological Josephson junctions . Wealso predict, based on our theory, that in a topological su-perconducting quantum interference device (SQUID) two in-terference patterns with periods h / e and h / e can coexist.This phenomenon will be an supporting evidence for Majo-rana modes if verified by future experiments. II. QUANTUM RESISTIVELY SHUNTED JUNCTIONMODEL
The topological Josephson junction sketched in Fig. 1aconsists of two topological superconductors, which couldbe one-dimensional nanowires with spin-orbit couplings ,superconducting quantum spin-Hall edge states , or ferro-magnetic atomic chains . The junction hosts two Majo-rana modes γ L , R with their coupling described by H M = − iE M γ L γ R cos ( θ / ) with E M the maximum coupling energy.By defining a Dirac fermion f = γ L + i γ R , the Hamiltoniandescribes a typical two-level system where the empty state | (cid:105) and occupied state | (cid:105) are the two eigenstates. The corre-sponding energy spectra are E ± = ± E M cos ( θ / ) which crossat θ = ( n + ) π . In finite-size materials, the inevitable over-lapping between γ L , R and the other two edge Majorana modesleads to hybridization of the two states (see Appendix A fordetails), which produces avoided energy crossings. By writ-ing the wave function as | ψ (cid:105) = ψ | (cid:105) + ψ | (cid:105) , the dynamicsis determined by the Schr¨odinger equation i ¯ h dd t (cid:18) ψ ψ (cid:19) = (cid:18) E M cos θ δδ − E M cos θ (cid:19) (cid:18) ψ ψ (cid:19) , (1)with δ the hybridization energy. This equation describes atwo-level system which has an energy spectrum with avoidedcrossings at θ = ( n + ) π , as illustrated in Fig. 1b. When θ is driven through the avoided crossings, the LZ transition be-tween the two levels will change the system from the groundstate to the excited state with a textbook LZ transition proba-bility P = e − πδ / ( ¯ h ˙ θ E M ) .We consider a junction with negligible capacitance, wherethe motion of θ under biased current can be described by theRSJ model , which is the current conservation equationwhere the total current I is transported through the resistiveand Josephson channel with I = V / R + I J as shown schemat-ically in Fig. 1c. The Josephson current I J has two parts: theconventional Cooper-pair channel I = I c1 sin θ , and the par-ity dependent Majorana channel I = I c2 (cid:104) ψ | i γ L γ R | ψ (cid:105) sin ( θ / ) which comes from the phase derivative of H M23 (see Ap-pendix A). By invoking the ac Josephson relation we obtainthe equation explicitly asd θ d t = eR ¯ h (cid:20) I − I c1 sin θ − I c2 (cid:0) | ψ | − | ψ | (cid:1) sin θ (cid:21) , (2) where the quantum average over Majorana operators is ex-pressed with the wave function. This equation brings non-linearity to the Schr¨odinger equation (1), and they togetherconstitute the QRSJ model.One important feature here is that when I is large enoughto make the right hand side of Eq. (2) nonzero, the motionof θ would induce the LZ transitions around θ = ( n + ) π .Different from the conventional LZ effect, the injected cur-rent drives the Josephson phase passing the avoided crossingsagain and again with a large velocity. Each time the LZ transi-tion only induces a small change on the two-component wavefunction. However, the accumulation of many LZ transitionsleads to a nonlinear quantum dynamics of the two-level sys-tem as we will show later. Therefore, the LZ transition is thebuilding brick of the complicated but nontrivial dynamics ofthe two-level system in the topological junction.To observe the effect of these LZ transitions, we first nu-merically integrate Eqs. (1) and (2) with initial conditions ψ = θ =
0, and present the time evolution of the wavefunction in Fig. 2a. We see that the wave function oscillates atthe full time range. Looking carefully, the oscillation ampli-tude begins from a small value with the system mainly stayingat | (cid:105) , and then gradually increases. After passing a criticaltime marked by the red dashed line, the wave function beginsto oscillate between | (cid:105) and | (cid:105) . We will see later that thiscritical time relates to passing the separatrix of an effectiveclassical Hamiltonian. We also notice that the oscillating pe-riod is shorter at the two ends of the time range, and becomeslonger nearby the critical time. Besides the rich oscillatoryfeatures, there is also an obvious damping on the envelope ofthe oscillations, with a characteristic time scale much largerthan the oscillation periods. The damped quantum oscillationis unique and reflects the impact of the nonlinear dynamicsof θ which enters the Schr¨odinger equation of the two-levelsystem. III. NONLINEAR DYNAMICS OF MAJORANATWO-LEVEL SYSTEM
Now we analyze this damped quantum oscillation by map-ping the QRSJ model to a nonlinear classical model, whichenables the usage of sophisticated approaches that have beendeveloped for solving nonlinear classical dynamics . Thetrick is to notice that in the QRSJ model the wave functionis subjected to two restrictions: it must be normalized, andthe global phase is decoupled from the dynamics (see Ap-pendix B for details). Then we can define two real variables:the relative amplitude s = | ψ | − | ψ | and the relative phase φ = arg ψ − arg ψ , which are complete for describing the dy-namics of the two-level system . With this trick, we castthe QRSJ model into a purely classical model and write downthe dynamical equationsd θ d t = eRI (cid:20) − I c1 I sin θ − sI c2 I sin θ (cid:21) , (3a)d s d t = − δ (cid:112) − s sin φ , (3b)d φ d t = E M cos θ + δ s √ − s cos φ , (3c)where we take the unit ¯ h = τ θ = / eRI , τ s = / δ and τ φ = / E M which correspond tothe change of θ , s and φ . We note that these three time scalesare different by orders with τ θ (cid:28) τ φ (cid:28) τ s for the junctionparameters shown in Fig. 2a and generally for I > I c1 + I c2 .For classical nonlinear systems with multiple time scales,the method of averaging is a powerful technique . Theessence is to categorize ”fast” variables and ”slow” variablesby typical time scales, then solve the equations for the fastvariables by treating slow variables as constant parameters.After obtaining the solution, the fast variables are averaged over its time scale and used for solving the equations of the (a) Separatrix (b) -1 -0.5 0 0.5 1-1-0.500.51 -1 -0.5 0 0.5 1-1-0.500.51 (c)
Analytical
Numerical
FIG. 2. (Color online) (a) Evolution of the wave function for the two-level system under constant injected current I / I c2 = .
5, obtained bynumerically solving Eqs. (1) and (2). The bottom inset is a zoom-inview in the marked time window. The analytical solution Eq. (10)provides the ˜ τ s , τ d , the dashed envelope line and the top inset. (b)Phase-space portrait of the classical Hamiltonian H c , with P the el-liptic fixed point, P the hyperbolic fixed point, and red-dashed cir-cle the separatrix. (c) Poincar´e map obtained by numerically solv-ing Eq. (3). Parameters of the junction are taken as I c1 / I c2 = . δ / E M = .
02, and R = h / e . slow variables. With this process, the dynamical equations aredecoupled into averaged equations, which significantly sim-plifies the problem.Now we use the method of averaging to analyze the nonlin-ear dynamics in Eq. (3), where θ is treated as the fast variableand s , φ as slow variables, since τ θ is the smallest time scale.We first consider s unchanged in τ θ and solve Eq. (3a) to ob-tain the time average of cos θ , defined as cos θ ≡ (cid:82) d t cos θ with integration range the time for θ to rotate 4 π .By taking the time derivative on both sides of Eq. (3a), wecan obtain terms containing s and ˙ s . Within τ θ , because s and˙ s both vary slowly, we take them as time independent. Withsome tedious but straightforward computation we obtain (seeAppendix C for details)cos ( θ / ) ≈ α s + β ˙ s , (4)with α = I c1 I c2 / I and β = I c2 τ θ / I from the lowest order Tay-lor expansion of I c1 / I and I c2 / I . Here α s is much larger than β ˙ s , and we refer them as zeroth-order and first-order averag-ing respectively.We begin from the zeroth-order averaging and replace cos θ with cos θ = α s in the Schr¨odinger equation (1), and obtain i ¯ h dd t (cid:20) ψ ψ (cid:21) = (cid:20) E M α ( | ψ | − | ψ | ) δδ − E M α ( | ψ | − | ψ | ) (cid:21)(cid:20) ψ ψ (cid:21) , (5) which becomes a typical nonlinear Schr¨odinger equation dueto the nontrivial diagonal elements . This explicitlyshows that coupling to Josephson phase dynamics brings thenonlinearity into the quantum dynamics of the two-level sys-tem, which is the reason for the rich and unusual dynamicalbehaviors shown in Fig. 2a (see Appendix C for details).Now we interpret this nonlinear quantum dynamics withthe classical model. In Eq. (3c) by replacing cos θ with itsaverage , we obtain,d φ d t = E M α s + δ s √ − s cos φ . (6)Now the system is only described by Eq. (6) and Eq. (3b) with θ integrated out. These two equations are the canonical equa-tions of a classical Hamiltonian (see Appendix C for details), H c = − α E M s + δ (cid:112) − s cos φ , (7)where s and φ are the coordinate and canonical momentum.Let us use the phase space portraits of this effective Hamil-tonian, as shown in Fig. 2b, to understand the oscillatory fea-tures shown in Fig. 2a. There is an elliptic fixed point P at ( s , φ ) = ( , ) , and a hyperbolic fixed point P at ( s , φ ) =( , ± π ) (see Appendix C for details). A separatrix connectsthe hyperbolic fixed point, separating the phase space into twodistinct areas: extended trajectories outside the separatrix andorbiting trajectories around the elliptic fixed point inside theseparatrix.The extended trajectories outside the separatrix in Fig. 2bcorrespond to dynamics before the critical time in Fig. 2a. Formotion along these trajectories, the s stays negative or posi-tive, agreeing with the small oscillations with | ψ | > | ψ | atthe beginning of Fig. 2a. Inside the separatrix, the trajecto-ries become orbital, with s oscillating from negative to pos-itive values. This corresponds to the oscillations in Fig. 2aafter the critical time, where | ψ | and | ψ | have overlappedoscillations. When approaching the separatrix, the period ofthe orbits is enlarged since the period should be divergent atthe separatrix . This corresponds to the observed period en-largement near the critical time in Fig. 2a. From the aboveanalysis, we argue that the system begins from outside of theseparatrix, passing through the separatrix at the critical time,and then orbits inside the separatrix and finally reaches theelliptic fixed point.For clarity we demonstrate the Poincar´e map of the numer-ical results for Eq. (3) in Fig. 2c, which is obtained by record-ing the points on the s − φ plane with θ = n π . The local traceof the Poincar´e map follows the trajectories of the classicalHamiltonian, illustrating that the oscillations shown in Fig.2a can be approximately determined by the classical Hamil-tonian. The global structure of the Poincar´e map, however,demonstrates a spiral-in feature from outside the separatrix tothe elliptic fixed point P . This exhibits the effect of a frictionforce which brings all phase-space trajectories to elliptic fixedpoints. This long-time-scale damping, also shown in Fig. 2a,cannot be obtained based on the zeroth-order averaging.Now we explore the damping feature by including the first-order averaging, and replacing cos θ / P , we find that Eqs. (3b) and (3c) leadto (see Appendix C for details)¨ s + β E M ˙ s + ( δ + α E M δ ) s = , (8)which is nothing but a classical damped harmonic oscillator.It has a standard solution of the form, s = e − t / τ d cos ( π t / ˜ τ s ) , (9)with the damping and oscillating time of τ d = eRI I c2 E M δ , ˜ τ s = πδ (cid:112) + α E M / δ . (10)We plot this analytical solution as an inset of Fig. 2a, and findthat it agrees well with the numerical simulations around theelliptic fixed point. Here we have demonstrated a duality be-tween the nonlinear quantum dynamics in this two-level sys-tem and a classical damped harmonic oscillator which is ex-actly solvable. Therefore, this duality enables us to find an an-alytical solution for the damped quantum oscillations despitethe equations for the nonlinear quantum dynamics is rathercomplicated. In fact, we further show a mapping to a solvableanharmonic damped oscillator (see Appendix C for details),which even correctly describes the dynamics far from the el-liptic fixed point. IV. HYSTERESIS IN I-V CURVES
Now we study the I-V characteristics of the topologicalJosephson junction based on the QRSJ model. We numeri-cally simulate the average voltage upon adiabatic current in-jection, which gradually increases to a large value and then (a) (b) (c) -6 -4 -2 0 2 4 6-0.3-0.2-0.100.10.20.3 (d)
FIG. 3. (Color online) I-V curves in absence of LZ effect for (a) I c2 =
0, (b) I c2 / I c1 = δ =
0. (c) I-V curves in presence ofLZ effect with parameters the same as in Fig. 2a. (d) Comparisonbetween our numerical simulation (solid lines) and the experimentaldata (discrete crosses) taken from Ref. [43]. Junction parametersin simulation are adopted the same as in the experiments with theresistance R = Ω , the capacitance C = π -period current I c1 = µ A, the 4 π -period current I c2 = . µ A, and the decoherencetime chosen as τ = ¯ h / E M . decreases back to zero. As a benchmark, we first show theI-V curve for a trivial junction with I c2 = V = R (cid:113) I − I around the criticalcurrent . We then consider an additional 4 π -period Joseph-son current I = I c2 sin θ which corresponds to the case of lo-cal parity conservation with δ =
0, where LZ effect cannottake place. We solve the Eq. (2) with | ψ | = | ψ | = π -period Josephson current modifiesthe shape of the I-V curve but demonstrates no novel phe-nomenon. For both cases, the voltage which is the velocity ofthe phase difference is fully determined by the applied current,so the quantum dynamics is history independent.However, when δ becomes finite and the LZ transitions be-gin to affect the tunneling current, we find an unambiguoushysteretic I-V curve with two critical currents as shown in Fig.3c: a switching current I sw where the voltage jumps from zeroto finite value and a smaller retrapping current I re for the finitevoltage jumping back zero.The origin of this hysteresis can be understood with thetime evolution of | ψ | and | ψ | discussed in Fig. 2a. Ini-tially for a small injected current below the switch value, thevoltage is zero and the two-level system stays at | (cid:105) with cer-tainty ( | ψ | = I J = I c1 sin θ + I c2 sin θ and the critical currentis given by I sw = ( I c1 ζ + I c2 ) (cid:112) − ζ , (11)with ζ = (cid:113) I / I + / − I c2 / I c1 (see Appendix A for de-tails). On the other hand in the current decreasing stage, be-cause the two-level system has finite probabilities on both lev-els due to the LZ transitions, the Josephson current changes to I J = I c1 sin θ + I c2 ( | ψ | − | ψ | ) sin θ . For this case, the crit-ical current would be smaller, since the two levels with oppo-site parities carry opposite currents and cancel each other. Ifthe cancellation is perfect with | ψ | = | ψ | , the correspondingcritical current is I re = I c1 . (12)Consequently, a hysteresis phenomenon emerges due to theexistence of the Majorana modes. We note that this hystere-sis requires neither local nor global parity conservation and isimmune to various quasiparticle poisoning effects in realisticsetups (see Appendix D, E for details).The hysteresis is solely due to the nonlinear dynamics ofthe two-level system formed by Majorana modes. Thereforewe would expect it to disappear after the topological phasetransition into the trivial superconducting phase. When thesystem approaches the transition point from the topologicalnontrivial side, the spatial spreading of Majorana modes in-creases, which gradually annihilates the hysteresis with twomechanisms. First, the overlapping of two Majorana modeson the same side of the junction increases, which greatly en-larges the coupling energy δ . The eigenstates become stateswith approximately equal weight of | (cid:105) and | (cid:105) . Thereforethe state of this two-level system for the current-increasingand -decreasing process become approximately the same, sothe hysteresis gradually disappears. Second, the weight of thewave function of Majorana mode at the edges become smaller.Correspondingly the tunneling current of the Majorana chan-nel I c2 decreases and so does the hysteresis.From the classical model described by Eq. (3), the hys-teresis is similar to the mechanical hysteresis from the dryfriction since Eq. (3a) is actually a friction equation. Thatis, the particle has different friction forces when it is static andmoving in the direction of θ . This difference comes from thehistory dependent trajectories in the s − φ plane (SeeAppendix B for details), and then feedback to the motion inthe θ direction through the last term in Eq. (3)a. This feed-back effectively induces a difference in the static friction anddynamic friction for the particle; therefore the particle wouldbegin and stop moving at different dragging forces. V. DIRECT COMPARISON WITH EXPERIMENTS
In recent experiments, hysteretic I-V curves have beenreported in a number of overdamped topological Joseph-son junctions, which are unexpected from the conventionalshunted junction theory . We argue that these hysteresisbehaviors possibly come from Majorana modes as we demon-strated from the QRSJ model. In order to prove our argument,we quantitatively compare our theoretical results with experi-mental results. For this reason, we consider the resistively andcapacitively shunted junction model, I = ¯ hC d θ e d t + ¯ h d θ eR d t + I c1 sin θ + I c2 (cid:104) i γ L γ R (cid:105) sin θ , (13)and the master equation for the two-level system (see Ap-pendix E for details),d ρ d t = − i ¯ h [ H , ρ ] + τ L , (14)where ρ is the density matrix of the two-level system, τ is thedecoherence time, and L = | ψ g (cid:105)(cid:104) ψ e | is the standard Lindbladform where | ψ e (cid:105) and | ψ g (cid:105) are the two instantaneous eigen-states of the two-level system. This combination of Eq. (13)and (14) can describe the small but nonzero capacitance andthe decoherence in experiments, however, it is too compli-cate for analytical solution. Here we numerically simulate themodel where the junction parameters are taken from the ex-perimental data , with resistance R = Ω and capacitance C = I c1 = µ A and I c2 = . µ A are extracted from the switching and retrappingcurrent of the experimental I-V curve . The decoherencetime is taken as τ = ¯ h / E M . It is much larger than othertime scales ( τ s , τ θ , τ φ ), which is reasonable because the deco-herence is suppressed by the superconducting gap . Theresult of the simulation is presented in direct comparison withthe experimental data as shown in Fig. 3d. Our theoreticalresults agree well with the experimental data. As far as weknow, the experimental results have no convincing explana-tion so far, and it has never been associated with the topologi-cal nature of the junction. Our results give a reasonable expla-nation for the experimentally reported “unexpected” hystere-sis from the aspect of Majorana modes. VI. INTERFERENCE PATTERN OF A TOPOLOGICALSQUID.
Hysteresis is also expected in a SQUID composed by twosuch junctions as shown in Fig. 4a, where the flux depen-dence of critical currents is a routine measurement . Thesame as for the single topological junction, the I-V curve ofthis SQUID should also be hysteretic. Then we expect twointerference patterns of maximum supercurrent, one for theswitching current and the other for the retrapping current. Theswitching current should contain contributions from both the (a) (b) FIG. 4. (Color online) (a) Schematic setup of a topological SQUIDstructure with four Majorana zero modes. (b) The analytical inter-ference pattern for the switching current (blue solid line) and the re-trapping current (orange solid line), and the numerically results forthe interference patterns of switching current (blue circle) and re-trapping current (orange diamond). Josephson currents are taken as I (cid:48) c1 / I c1 = I (cid:48) c2 / I c2 = .
4. Other Parameters are taken the same as Fig.3c for two identical junctions. conventional and Majorana channel and is thus given by I sw ( Φ ) = max θ (cid:2) I c1 sin θ + I (cid:48) c1 sin ( θ + π ΦΦ ) (15) + I c2 sin θ + I (cid:48) c2 sin ( θ + π ΦΦ ) (cid:3) , where I c1 and I c2 represent the supercurrent for the quasipar-ticle and Majorana channels in one junction, I (cid:48) c1 and I (cid:48) c2 repre-sent the supercurrent for the quasiparticle and Majorana chan-nels in the other junction, Φ is the magnetic flux through theSQUID, and Φ = h / e is the superconducting flux quantum.Here we require the total parity conservation of the coupledMajorana modes. This interference pattern, as shown explic-itly in Fig. 4b, is obviously 2 Φ -periodic, which agrees withprevious studies . On the other hand, the currents fromMajorana channels are almost canceled when considering theretrapping current, which leads to I re ( Φ ) ≈ max θ (cid:2) I c1 sin θ + I (cid:48) c1 sin ( θ + π Φ / Φ )] , (16)which is Φ -periodic as shown in Fig. 4b. I sw and I re can bedirectly obtained by numerically studying the dynamics withthe QRSJ model, where the Hamiltonian for the coupled Ma-jorana modes in the SQUID is H = − i γ γ E u cos ( θ / ) − i γ γ E d cos [( θ + π Φ / Φ ) / ]+ i δ l γ γ + i δ r γ γ , (17)with E u , d and δ l , r the corresponding coupling coefficients. Thenumerical results are shown in Fig. 4b, which agree well withour analytical results.From both the analytical and numerical results, in a topo-logical SQUID we can obtain coexistence of h / e and h / e -periodic interference patterns, which as far as we know isnever seen in any SQUID before. The physical reason behindthis phenomenon is that the Majorana channel contributesonly to switching current but negligibly to retrapping current.This unique interference phenomenon, if experimentally veri-fied, will be an evidence for the existence of Majorana modes. VII. CONCLUSION
In summary, we propose that the Landau-Zener effect ofthe two-level system in a topological Josephson junction canlead to hysteresis in the I-V characteristics. We establish aquantum resistively shunted junction model to study the prob-lem. We demonstrate the nonlinear quantum oscillation in thetwo-level system of the junction, with both numerical simula-tion and analytical methods, and show that the hysteretic I-Vcurves naturally follows from it. We compare our theoreti-cal results with existing experimental results and find themquantitatively in agreement. We predict coexistence of h / e -periodic and h / e -periodic interference patterns which aresubjected to further experimental verifications. ACKNOWLEDGMENTS
The authors are grateful for Pavan Hosur, Stefan Lud-wig, Chin-Sen Ting and Hongqi Xu for helpful discus-sions. This work was supported by Grants Nos. NKRDPC-2017YFA0206203, 2017YFA0303302, 2018YFA0305603,the National Natural Science Foundation of China underGrants Nos. 11774435 and No. 61471401, and China Schol-arship Council under Grants No. 201706385057. Zhao Huangis supported by Robert A. Welch Foundation under GrantNo. E-1146. Qian Niu is supported by DOE (DE-FG03-02ER45958, Division of Materials Science and Engineering),NSF (EFMA-1641101) and Robert A. Welch Foundation (F-1255).
Appendix A: Josephson Hamiltonian and Josephson current
Here we present a derivation for the 2 π -period Joseph-son current I c1 , the 4 π -period Josephson current I c2 , and theHamiltonian of the two-level system H M . In realistic topolog-ical Josephson junctions, usually there are both topologicaland non-topological segments . For example, in the topo-logical superconducting nanowire as sketched in Fig. 1, thewire is topological and the substrate s-wave superconductoris non-topological. The topological segment carries the 4 π -period Josephson current due to Majorana modes while thenon-topological segment carries the 2 π -period Josephson cur-rent.Here we use a phenomenological model to describe aJosephson junction with both topological and non-topologicalsegments. It is a hybrid two-layer system with one layer asa spinless Kitaev chain and the other layer as a trivial s-wavesuperconductor.We first consider the trivial layer which is described by asimple Hamiltonian as H α = − t α ∑ (cid:104) i , j (cid:105) , α , σ c † i , α , σ c j , α , σ − µ α ∑ i , α , σ c † i , α , σ c i , α , σ + ∑ i , α ( ∆ α e i θ α c † i , α , ↑ c † i , α , ↓ + h . c . ) , (A1)where α = L , R represents the left and right sides of the wire, c α , j , σ is the electron annihilation operator on the site j andspin σ = ↑ , ↓ , ∆ α is the superconductor gap, θ α is the super-conducting phase, t α is the nearest neighbor hopping, and µ α is the chemical potential. Here for simplicity we take identicalparameters for the left and right segments, except for the su-perconducting phase θ α which must be different in the pres-ence of a Josephson current. The two superconductors areconnected with a tunneling Hamiltonian, H T = ∑ σ ( T c † L , σ c R , σ + h . c . ) , (A2)where c † L , σ is the electron creation operator at the boundary ofthe left superconductor nearby the junction, T is the tunnelingstrength which is determined by the tunneling barrier of thejunction. In a realistic junction, this can be controlled by anapplied gate voltage. The Josephson current can be calculatedwith the standard Green function technique, where the currentis expressed as, I = eT Im [ ∑ k , p , i ω ℑ † ( k , i ω ) ℑ ( p , i ω )]= I c1 sin θ , (A3)where θ = θ L − θ R is the Josephson phase, ℑ is the off-diagonal Matsubara Green function, and I c1 is given by thecontour integral as, I c1 ≈ e ∆ T ( − µ / t ) ¯ ht . (A4)We note that I c is a square function of T which reflects theCooper-pair tunneling. Higher order contributions in the S-matrix expansion can also be included, however, they shouldbe negligible for the tunneling regime where the tunneling T is small compared with the hopping t .Now we consider topological layer, which can be stud-ied with a spinless p-wave superconducting Hamiltonian pro-posed by Kitaev , H α = N α ∑ j = (cid:104) − t α c † α , j c α , j + + ∆ α e i θ α c α , j c α , j + + h . c . (cid:105) − µ α N α ∑ j = c † α , j c α , j . (A5)In this model, the electron operators can be transformedto Majorana operators γ α , j , A = e i θ α / c α , j + e − i θ α / c † α , j and γ α , j , B = − ie i θ α / c α , j + ie − i θ α / c † α , j . Then the Hamiltoniancan be rewritten in this Majorana representation, H α = ( t + ∆ ) N − ∑ j = i γ α , j , B γ α , j + , A − ( t − ∆ ) N − ∑ j = i γ α , j , A γ α , j + , B − µ α N ∑ j = i γ α , j , A γ α , j , B . (A6)It is well known that this Kitaev model enters the topologicalnon-trivial phase for the parameter regime of | t | > | µ | and ∆ (cid:54) =
0, while the Majorana modes γ L , γ (cid:48) L , γ R , and γ (cid:48) R appearsat the ends of the two segments . Then the low energy (belowsuperconducting energy gap ∆ ) physics of the two segments isdescribed by an effective Hamiltonian, H δ = ∑ α i δ α γ (cid:48) α γ α , (A7)where δ α represents the coupling energy within the left/rightsegment, which is exponentially protected by the length of thewire .The two segments are coupled by the electron tunnelingthrough the barrier, which could be described by a standardtunneling Hamiltonian H T = T c †L , N c R , + T ∗ c †R , c L , N . (A8)For low energy physics, the effective Hamiltonian shouldonly involve the four Majorana modes. Therefore the tunnel-ing Hamiltonian should be projected to these four Majoranamodes with a form of , H M = − iE M γ L γ R cos ( θ / ) , (A9)with E M ≈ T / f = ( γ L + i γ R ) / f = ( γ (cid:48) R + i γ (cid:48) L ) / H = H M + H δ = − E M cos ( θ / )( f †1 f − f f †1 )+ δ L ( f − f †2 )( f + f †1 ) + δ R ( f + f †2 )( f − f †1 ) . (A10)There are natural basis states for this Hamiltonian: | (cid:105) , f †1 f †2 | (cid:105) , f †2 | (cid:105) , and f †1 | (cid:105) , with | (cid:105) the vacuum state for f †1 and f †2 . With these basis states, the total Hamiltonian canbe rewritten in the matrix form as, H = E M cos ( θ / ) δ L + δ R δ L + δ R − E M cos ( θ / ) E M cos ( θ / ) − δ L + δ R − δ L + δ R − E M cos ( θ / ) . (A11)This is a block diagonal matrix, with the left-up and right-down blocks corresponding to the even and odd total parities,respectively. Without losing generality, we take the even totalparity and arrive at the matrix shown in Eq. (1) of the maintext with δ = δ L + δ R .Now let us consider the Josephson current through the Ma-jorana channel. The electron number operator on the right-hand side of the junction is N R = ∑ j c †R , j c R , j , and its timederivative gives the tunneling current, I ( t ) = − e (cid:104) d N R d t (cid:105) = − e (cid:104) ψ ( t ) | i ¯ h [ H , N R ] | ψ ( t ) (cid:105) = ie ¯ h (cid:104) ψ ( t ) | − T c †L , N c R , + T ∗ c †R , c L , N ) | ψ ( t ) (cid:105) , (A12)where | ψ ( t ) (cid:105) is the ground state wave function after includ-ing the tunneling Hamiltonian. The single electron tunnelingthrough Majorana modes is obtained by the zero-order degen-erate perturbation as, I = I c2 sin ( θ / ) (cid:104) ψ ( t ) | i γ L γ R | ψ ( t ) (cid:105) , (A13)with the maximum value I c2 ≈ eE M ¯ h = eT h . (A14)Here we notice that I c2 is linear in T which reflects the phasecoherent single electron tunneling. Comparing with Eq. (A4),we obtain the ratio between the amplitude of the 2 π -periodsupercurrent and the 4 π -period supercurrent I c1 I c2 ≈ ∆ T ( − µ / t ) t . (A15)We notice that it is a linear function of the tunneling strength T . That is, the 4 π -period supercurrent contributed by Majo-rana modes dominants the transport for junctions with hightunneling barriers, while the 2 π -period supercurrent con-tributed by quasiparticles dominants the transport for junc-tions with high transparency.Finally we give a derivation for the switch current shownin Eq. (11). It is the maximum current when s = θ is a free variable. We first calculate the Josephson phase forachieving the maximum current, which is denoted as θ c . It isobtained by taking phase derivative of the Josephson currentdd θ c (cid:18) I c1 sin θ c + I c2 sin θ c (cid:19) = , (A16)which gives ζ ≡ cos θ c = (cid:113) I r + / − I r , (A17)with I r = I c1 / I c . We plug it back to the expression for theJosephson current and obtain I sw = I c1 sin θ c + I c2 sin θ c = (cid:112) − ζ ( I c1 ζ + I c2 ) , which gives the Eq. (11). Appendix B: Casting the two-level system to a classicalHamiltonian
Now we demonstrate how to cast the Schr¨odinger equationfor the two-level system i ¯ h dd t (cid:18) ψ ψ (cid:19) = (cid:18) E M cos θ δδ − E M cos θ (cid:19) (cid:18) ψ ψ (cid:19) (B1)into classical equations, and form a classical dynamical sys-tem by combining with the equation for Josephson phase fromresistively shunted junction model. The wave function of thetwo-level system is ( ψ , ψ ) T ≡ ( | ψ | e i φ , | ψ | e i φ ) T which contains two complex numbers. However, it obeys two con-straints. First, it must be normalized | ψ | + | ψ | =
1; sec-ond, the overall phase of the wave function is decoupled fromthe dynamics of the two-level system. With these constraints,the wave function can actually be described by two real dy-namical variables. One convenient choice is the relative am-plitude s ≡ | ψ | − | ψ | and the relative phase φ = φ − φ .Now we derive the equations for these two real variables outof the Schr¨odinger equation. For this purpose, we explicitlywrite down the amplitude and the phase of the wave functionusing s and φ . The amplitude of the wave function is deter-mined by s with | ψ | = (cid:112) ( − s ) / | ψ | = (cid:112) ( + s ) / φ and the total phase φ T = φ + φ with φ =( φ T − φ ) / φ = ( φ T + φ ) /
2. Then we can transform theSchr¨odinger equation into the form, i ¯ h dd t (cid:32) (cid:113) − s e − i φ / (cid:113) + s e i φ / (cid:33) e i φ T / = (cid:18) E M cos θ δδ − E M cos θ (cid:19) (cid:32) (cid:113) − s e − i φ / (cid:113) + s e i φ / (cid:33) e i φ T / , (B2)We note that we have added a factor of 1 / δ and E m are rescaled to be doublingtheir original value. We reach at two complex equations forthe real variables s , φ , and φ T . The first equation is, i ¯ h ( − (cid:115) ( − s ) ˙ s − i (cid:114) − s φ + i (cid:114) − s φ T )= E M θ (cid:114) − s + δ (cid:114) + s e i φ . (B3)The imaginary part of the equation gives,˙ s = − δ ¯ h (cid:112) − s sin φ , (B4)which is the Eq. (3b), while the real part of the equation gives,˙ φ − ˙ φ T = E M ¯ h cos θ / + δ √ + s ¯ h √ − s cos φ . (B5)Checking the second equation we would have,˙ φ + ˙ φ T = E M ¯ h cos θ / − δ √ − s ¯ h √ + s cos φ . (B6)Combining the Eqs. (B5) and (B6), we obtain the Eq. (3c),˙ φ = E M ¯ h cos θ / + δ s ¯ h √ − s cos φ . (B7)Rearranging the formulas we arrive at Eq. (3). We have twoequations for the two-level system,d s ( t ) d t = − δ ¯ h (cid:113) − s ( t ) sin φ ( t )= − τ s (cid:113) − s ( t ) sin φ ( t ) , (B8)and d φ ( t ) d t = E M ¯ h cos θ ( t ) + s ( t ) δ ¯ h (cid:112) − s ( t ) cos φ ( t )= τ φ cos θ ( t ) + s ( t ) cos φ ( t ) τ s (cid:112) − s ( t ) , (B9)and one equation for the Josephson phase,d θ ( t ) d t = eR ¯ h (cid:20) I − I c1 sin θ ( t ) − I c2 s ( t ) sin θ ( t ) (cid:21) = τ θ (cid:20) − I sin θ ( t ) − I s ( t ) sin θ ( t ) (cid:21) , (B10)where τ s = ¯ h / δ , τ φ = ¯ h / E M , τ θ = ¯ h / eRI , and we redefinetwo dimensionless parameters I ≡ I c / I and I ≡ I c / I formathematical simplicity. We see that φ T is decoupled fromthese three equations. -1 0 1-101 (a) -1 0 1-101 (b) -1 0 1-101 (c) FIG. 5. (Color online) Typical trajectories of the particle in the phasespace for the injected current of (a) I / I c1 = . I / I c1 = . I / I c1 = Within this pure classical model, we first analyze the dy-namical stability of the junction with the injected current I asthe control parameter. As shown in Fig. 5, we numericallyexplore three different injected currents. For a small current I < I c1 , the trajectories for all initial conditions are closed,demonstrating circles in the s − φ plane as seen in Fig. 5a.For an intermediate current I = . I c1 , there are two differ-ent types of trajectories, depending on the initial conditions.The trajectory for large initial s is closed while the trajectoryfor small initial s is not closed, falling to s ≈ I = I c1 , all trajectories are falling to s ≈ s . This historydependence suggests the effect of nonlinearity in the dynami-cal evolution and the falling of | s | indicates existence ofdamping mechanism. In the following, we adopt the methodof averaging to analytically study the Eq. (3). Appendix C: Method of averaging
After casting the QRSJ model into a purely classical model,we obtain a set of classical nonlinear equations. At first sight,the new classical equations for s and φ are no simpler than the Schr¨odinger equation in the original QRSJ model. However,the advantage of this pure classical formalism is the availabil-ity of sophisticated mathematical approaches that have beendeveloped to study nonlinear classical dynamics. From Eq.(3) we have extracted three typical time scales which are dif-ferent by orders with τ θ (cid:28) τ φ (cid:28) τ s . For nonlinear dynamicalsystems with multiple time scales, the method of averagingis a powerful mathematical tool . It was initially devel-oped by Krylov and Bogoliubov to tackle nonlinear oscilla-tion problems such as the study of the Einstein equation forMercury , and from then on the method has been found use-ful in many physical systems involving oscillations . Theessence of the method of averaging is to categorize the dynam-ical variables as ‘fast’ variables and ‘slow’ variables depend-ing on their typical time scales of variation. Then the slowvariables are regarded as almost unchanged within the timescale of the fast variables, and the time dependence of the fastvariables can be solved with the slow variables as fixed param-eters. After obtaining this time dependence, the fast variablesare averaged over time and the averaged values are pluggedback into the dynamical equations for the slow variables. Fi-nally, the time dependence of the slow variables can be solvedwith these averaged values of fast variables as external param-eters.The method of averaging allows us to study the dynamics offast variables and slow variables one by one, which is mucheasier than investigating the full complicate coupled nonlin-ear equations. In the following analysis, we can treat θ asthe fast variable and ( s , φ ) as the slow variables. We will seethat within τ θ , the zeroth-order averaging which uses a time-independent s to replace the function s ( t ) , is enough to givethe high-frequency oscillation shown in Fig. 2a. The first-order averaging, where a time-independent ˙ s is also taken intoaccount within τ θ , is capable of reproducing the damping fea-ture.
1. Time Averaging over Fast Variable
As seen in Eq. (B9), the fast variable θ enters the dynamicsof the slow variables through the function cos θ /
2. Now, wetry to calculate the time average for this function. The wholetime of dynamics can be cut into fractions of the time scale forthe fast variable τ θ . The slow variable s ( t ) should be almostunchanged within each fraction τ θ . As a zeroth-order averag-ing, s ( t ) is treated as time independent in the equation for θ .Therefore Eq. (B10) becomes,d θ ( t ) d t = τ θ (cid:20) − I sin θ ( t ) − I s sin θ ( t ) (cid:21) , (C1)which can be solved alone without considering the Eqs. (B8)and (B9) at the moment. The time evolution of θ can be ob-tained by solving only one differential equation, and after-wards we can make time averaging over the cos θ / θ ≡ T θ (cid:90) T θ d t cos θ ( t ) , (C2)0where T θ is the time for θ to increase 4 π which is at the orderof τ θ . Obviously this is a function of the parameter s . Here wetake a simple approach to evaluate the average without solvingEq. (C1) explicitly. We replace the integration over time withan integration over phase θ ,cos θ = (cid:82) π θ ˙ θ cos θ (cid:82) π θ ˙ θ = (cid:82) π d θ cos θ − I sin θ − I s sin θ (cid:82) π d θ − I sin θ − I s sin θ , (C3)where the solution of θ ( t ) from Eq. (C1) is used implicitly toaccomplish the transformation. We note that the time averagecos θ is nonzero because θ is not linear in time.This expression for the time average cos θ only contains s ,therefore corresponds to the zeroth-order averaging. Now wego to first-order averaging by including the influence of ˙ s . Letus derive this term by taking time derivative to Eq. (B10),d θ ( t ) d t = τ θ ( − I cos θ ( t ) − I s ( t ) cos θ ( t ) ) ∗ ( − I sin θ ( t ) − I s ( t ) sin θ ( t ) ) − τ θ I ˙ s ( t ) sin θ ( t ) . (C4)Within τ θ , since s ( t ) and ˙ s ( t ) vary slowly, we consider bothof them as time independent and mathematically replace thedynamical variables with static parameters s ( t ) ≈ s and ˙ s ( t ) ≈ ˙ s , which is the first order approximation as ˙ s is now takeninto account. We note that the time dependent velocity ˙ s ( t ) reverses the sign under time reversal operation t → − t , whilethe parameter ˙ s stays the same. We thus have a plus/minusambiguity in the replacement ˙ s ( t ) ≈ ± ˙ s . Now we arrive at adynamical equation for θ as, τ θ d θ ( t ) d t = (cid:18) − I cos θ ( t ) − I s cos θ ( t ) (cid:19) ∗ (cid:18) − I sin θ ( t ) − I s sin θ ( t ) (cid:19) ± I τ θ ˙ s sin θ ( t ) = − ∂ V ( s , ˙ s , θ ) ∂ θ , (C5)with the potential function V ( s , ˙ s , θ ) = − (cid:18) − I sin θ ± I s sin θ (cid:19) ± I τ θ ˙ s cos θ . (C6)This resembles a Newtonian equation for a particle with mass τ θ moving under a potential V , therefore obeys a conservationlaw within each τ θ , E = τ θ ˙ θ + V ( s , ˙ s , θ ) , (C7) which gives a solution for ˙ θ as˙ θ = ± τ θ (cid:112) [ E − V ( s , ˙ s , θ )]= ± τ θ (cid:115) E + (cid:18) − I sin θ − I s sin θ (cid:19) ± I τ θ ˙ s cos θ . (C8)This formula is the first-order averaging, and should recoverthe formula for zero-order approximation (Eq. (C1)) when˙ s =
0. This constraint requires the plus sign in front of thesquare root at the right hand side of Eq. (C8) and the energyto be E =
0, which leads to˙ θ = τ θ (cid:115)(cid:18) − I sin θ − I s sin θ (cid:19) ± I τ θ ˙ s cos θ . (C9)With this formula for ˙ θ , we can analytically calculate the timeaverage by transforming the integration over time to the inte-gration over θ ,cos θ = (cid:82) π θ ˙ θ cos θ (cid:82) π θ ˙ θ = (cid:82) π d θ cos θ (cid:113) ( − I sin θ − I s sin θ ) ± I τ θ ˙ s cos θ (cid:82) π d θ (cid:113) ( − I sin θ − I s sin θ ) ± I τ θ ˙ s cos θ . (C10)These two integral expressions Eq. (C3) and Eq. (C10) givethe time average over the fast variable up to zeroth-order andfirst-order averaging. Certainly we can go further to includethe influence of ¨ s , etc . However, we find that the s and ˙ s de-pendence is enough to qualitatively understand the dynamics.We will use these two integral expressions to obtain an ex-plicit function and plug it back into the dynamical equationsfor the slow variables.
2. The Zeroth-order Averaging and the Classical Hamiltonian
Let us first examine the zeroth-order averaging where thetime average is given by Eq. (C3). Now we calculate theintegrals by taking Taylor expansions,11 − I sin θ − I s sin θ = + ( I sin θ + I s sin θ ) + ( I sin θ + I s sin θ ) +( I sin θ + I s sin θ ) + ..., (C11)which gives the lowest order result for the denominator of Eq.(C3), (cid:90) π d θ − I sin θ − I s sin θ ≈ (cid:90) π d θ = π . (C12)1Similarly we have the numerator of Eq. (C3) as, (cid:90) π d θ cos θ − I sin θ − I s sin θ = (cid:90) π d θ cos θ (cid:20) + ( I sin θ + I s sin θ )+( I sin θ + I s sin θ ) + ... (cid:21) . ≈ I I s (cid:90) π d θ cos θ θ θ = π I I s . (C13)For both integrations, we take the lowest order nonzero termin the expansion series. Putting the results for numerator anddenominator together, we obtaincos θ ≈ I I s ≡ α s (C14)where for simplicity we define a parameter α = I I / θ in Eq. (B9) with the time averaged functioncos θ , and obtain the equations soly for s and φ asd s d t = − τ s (cid:112) − s sin φ , (C15)d φ d t = α s τ φ + s cos φτ s √ − s . With θ averaged out, obviously these two equations are selfconsistent equations. In fact they are the canonical equationsof a classical Hamiltonian, H c = − τ φ α s + τ s (cid:112) − s cos φ , (C16)where s and φ are the extended coordinate and the canonicalmomentum. This classical Hamiltonian is the Eq. (7) whichrepresents a classical integrable system, with the evolution ofthe phase-space motions of the Hamiltonian shown in Fig. 2b.Let us examine basic features of this classical Hamiltonian.We first study the fixed points ( s c , φ c ) , which are obtained bytaking the stationary condition of the Hamilton equations,d s d t (cid:12)(cid:12)(cid:12)(cid:12) s c , φ c = − τ s (cid:113) − s sin φ c = φ d t (cid:12)(cid:12)(cid:12)(cid:12) s c , φ c = τ φ E M α s c + s c cos φ c τ s (cid:112) − s = . Considering the fact that δ < E M α , there are three sets of fixedpoints, ( s c , φ c ) = ( , ) , ( , ± π ) , ( ± (cid:113) − τ φ / ( τ s α ) , ± π ) . (C18)The first two sets of fixed points are marked as P , P in the 2b.Then we check the classification of these fixed points, which is described by the Jacobian matrix at the fixed points, J ( s c , φ c ) = (cid:32) ∂ ˙ s ∂ s ∂ ˙ s ∂φ∂ ˙ φ∂ s ∂ ˙ φ∂φ (cid:33) s = s c , φ = φ c = s c sin φ c τ s √ − s c − τ s (cid:112) − s c cos φ c ατ φ + cos φ c τ s √ − s c + s cos φ c τ s ( − s c ) / − s c sin φ c τ s √ − s c . (C19)For the fixed point P at the position ( s c , φ c ) = ( , ) , we havethe Jacobian matrix of J ( , ) = (cid:32) − τ s ατ φ + τ s (cid:33) , (C20)This stability matrix has two imaginary eigenvalues λ , = ± i (cid:112) ( ατ s + τφ ) / τ s τ φ , signifying that P is an elliptic fixedpoint. Similarly, we calculate eigenvalues of the Jacobian ma-trices at the other two fixed points, and find that the P is a hy-perbolic fixed point, while ( s c , φ c ) = ( ± (cid:113) − τ φ / ( τ s α ) , ± π ) are elliptic fixed points. These analytical results agree with theinformation we see on the phase space portrait shown in Fig.2b.In classical dynamics, the phase-space trajectory whichconnects the hyperbolic fixed points is called separatrix. Inthis classical Hamiltonian Eq. (7), a separatrix connects thefixed points P , as shown in Fig. 2b. The separatrix dividesthe phase space into distinct regions, where the evolution ofthe phase-space motion orbits around different elliptic fixedpoints. In the action-angle formalize , the phase space areaenclosed by an orbit defines the action, I ( H c ) = π (cid:90) s ( H c , φ ) d φ , (C21)where the s is a function of φ and the energy H c by reversingEq. (7). This action is an adiabatic invariant and its derivativeon energy gives the period of oscillation T = d I ( H c ) d H c . (C22)Since the hyperbolic fixed points locate at energy saddlepoints, the separatrix has the diverging energy derivative.Therefore, we would expect a slow down of the oscillationif the motion is going near to the separatrix, which is clearlyseen in Fig. 2a.The classical Hamiltonian Eq. (7) obtained by the methodof averaging captures a number of features of the simulationresults. The phase space of the classical Hamiltonian is di-vided into two distinct areas by the separatrix. The orbits out-side the separatrix only have small oscillations in s with itsvalue remains positive or negative, while the orbits inside theseparatrix can oscillate between negative minimums and pos-itive maximums. These two distinct types of orbits agree withthe dynamics of the s shown in Fig. 2a, with the oscillationfirst small and only at the negative value, and later becominglarge and between negative and positive values. Interestingly,right at the transition between these two distinct oscillating2behaviors, we observe obvious enlargement of the period inFig. 2a. This indicates that the system is walking through theseparatrix which has the divergent period. Comparing the Fig.2a and Fig. 2b, it is reasonable to argue that the oscillationbehaviors are well described by the effective classical Hamil-tonian, but the damping of the oscillating amplitude cannot beunderstood yet.One direct method to view the resemblance between theclassical Hamiltonian and the original system is to draw thePoincar´e map for the evolution of the motion obtained by nu-merically solving the QRSJ model. The Poincar´e map is theintersections of a chosen surface in the phase space, called asthe Poincar´e surface of section, and the motion trajectories inthe whole phase space . This approach replaces the integra-tion of equations with the study of mappings, and has shownmuch power in nonlinear dynamics. It is particularly advan-tageous in understanding the qualitative features of the sys-tem. In our present case, we naturally choose the s − φ planewith θ = P for the classical Hamiltonian. However, ifexamined more carefully, the points actually spiral to the fixedpoint P . This is consistent with the large number of pointsaround P , which indicate the convergence of the trajectoriesand the breaking down of the phase-space volume conserva-tion. (a) (b) FIG. 6. (Color online) The time evolution of the wave function simu-lated for the nonlinear Schr¨odinger equation Eq. (5), with initial con-dition of (a) ψ = √ . , ψ = √ .
1, and (b) ψ = √ . , ψ = √ . Finally, it is inspiring to directly looking at the Schr¨odingerequation by replacing cos θ with the averaging cos θ in the Eq.1. We can obtain the Eq. (5), which is a nonlinear Schr¨odingerequation where other nontrivial LZ phenomena have been dis-cussed before . From this equation, we see that the cou-pling between the two-level system and the Josephson phasenaturally induces the nonlinearity to the quantum dynamics,which is the origin of such a rich and unusual dynamics forthe wave function of the two-level system as shown in Fig.2a. The dynamical of this nonlinear Schr¨odinger equation is initial value dependent, as expected from the equivalent classi-cal Hamiltonian. We numerically simulate the time evolutionof the wave function with two typical initial values, and showthe results in Fig. 6. We find that the wave function exhibitsoscillation patterns similar to the oscillations at different timeranges in Fig. 2a. However, the damping of the oscillatingamplitude is missing, which will be explained in the follow-ing section.
3. The first-order averaging and the damped harmonicoscillator
The classical Hamiltonian Eq. (7) helps us to understandthe quantum oscillation of the two-level system. However, itcan not describe the damping of the oscillation as shown inFig.2a. From the Poincar´e map shown in Fig. 2c, we knowthat the damping is towards the elliptic fixed point P . The nat-ural guess is that it comes from an extra friction force whichis proportional to the velocity of the extended coordinate ˙ s .This could be obtained from the expression of the first-orderaveraging Eq. (C10). Similar to the calculation of the zeroth-order averaging, the lowest-order Taylor expansion for the in-tegrated function gives the resultcos θ = α s + β ˙ s + O ( s , ˙ s ) , (C23)where β ˙ s represents the small contribution from the first-orderaveraging with β = ± I τ θ . If we plug this averaging resultback to Eq. (3c), we will find that the second term is linearin ˙ s , thus adding a velocity dependent force beyond the clas-sical Hamiltonian Eq. (7). For classical mechanical systemswhere forces only depend on coordinates, the Liouville theo-rem guarantees the phase-space volume conservation, whichguarantees undamped oscillations. On the other hand, the ex-istence of the velocity dependent force stemming from the β ˙ s term in Eq. (C23) breaks the Liouville theorem and the phase-space volume conservation. This is why the Poincar´e map inFig. 2c shows a phase-space volume compression, leading alltrajectories toward the elliptic fixed point P .Now we explicitly show that this ˙ s dependent term inducesdamping to the oscillation. For this purpose, we need to de-couple equations for s and φ . Taking time derivative on bothsides of Eq. (3b), we obtain¨ s = − δ ¯ h (cid:20) − s sin φ √ − s ˙ s + (cid:112) − s cos φ ˙ φ (cid:21) . (C24)Then we put Eqs. (3b) and (3c) into the right side of thisequation to eliminate the ˙ s and ˙ φ and arrive at¨ s = − δ ¯ h (cid:20) − s sin φ √ − s (cid:18) − δ ¯ h (cid:112) − s sin φ (cid:19) + (cid:112) − s cos φ (cid:18) E M ¯ h cos θ + s δ ¯ h √ − s cos φ (cid:19)(cid:21) = − δ ¯ h (cid:20) δ ¯ h s + E M ¯ h cos θ (cid:112) − s cos φ (cid:21) . (C25)3To eliminate φ from the right side of the equation, we noticethat Eq. (3b) can be transformed with trigonometric identityas ˙ s = δ ¯ h (cid:104) ( − s ) − ( (cid:112) − s cos φ ) (cid:105) , (C26)which can be used to replace the φ depend term and we obtain¨ s = − τ s s ± τ s τ φ cos θ (cid:113) − s − ( τ s ˙ s ) , (C27)where the ambiguity of the plus/minus sign comes from takingthe square root. Now we plug the averaging result Eq. (C23)into the equation, and obtain¨ s + (cid:18) τ s + τ s τ φ α (cid:113) − s − ( τ s ˙ s ) (cid:19) s ≈ − I τ θ τ s τ φ (cid:113) − s − ( τ s ˙ s ) ˙ s , (C28)where the correct plus/minus signs are chosen for obtainingconsistent results with numerical simulations. This second or-der differential equation represents a damped oscillator, wherethe angular frequency and the damping ratio depend on s . Thedamping comes from the right side of the equation which is afriction term proportional to ˙ s .By considering the regime of s , τ s ˙ s (cid:28)
1, we have the ap-proximation (cid:112) − s − ( τ s ˙ s ) ≈
1. Then Eq. (C28) can befurther simplified to¨ s + (cid:18) τ s + τ s τ φ α (cid:19) s + I τ θ τ s τ φ ˙ s = , (C29)which has exactly the same form as a classical damped har-monic oscillator. It can be rewritten to the standard form of,¨ s + ξ ω ˙ s + ω s = , (C30)with an angular frequency of ω = τ s + ατ s τ φ , (C31)and a damping ratio of ξ = I τ θ τ φ (cid:112) + ατ s / τ φ . (C32)This damped harmonic oscillator is underdamped with a smalldamping ratio ξ (cid:28) τ θ (cid:28) τ φ . Finally we arrive at thesolution for s around the elliptic fixed point as, s ( t ) = e − t / τ d cos ( t / ˜ τ s ) , (C33)with τ d = ξ ω = τ s τ φ I τ θ , ˜ τ s = ω (cid:112) − ξ ≈ ω = τ s (cid:112) + ατ s / τ φ , (C34)where we used ξ (cid:28) s . The new time scale τ d is the largest time scale which canbe constructed from the three basic time scales of the system.We compare the analytical solution given by Eq. (C33) withthe numerical results directly from Eq. (3) for several differentsets of parameters, as shown in Fig. 7. We find quantitativeagreement between them when s approaches zero.From above, we find that the quantum dynamics of the two-level system is dual to the classical dynamics of a dampedharmonic oscillator after adopting the method of averaging.This helps us to successfully obtain an analytical solution ofthe quantum dynamics of the two-level system within the sep-aratrix of the effective Hamiltonian. This dual relation givesa new insight in studying quantum two-level systems whennonlinearity is introduced.
4. Anharmonic damped oscillator
The solution of the damped harmonic oscillator obtained inEq. (C33) only becomes accurate when s approaches zero.Here we show an improved approximation which works atlarger s . Based on the numerical results and the analyticalsolution Eq. (C33), we know that the solution is a form ofdamped oscillation. Therefore, we propose an ansatz solutionof the form s = A ( t ) cos ( t / τ (cid:48) s ) , (C35)where A ( t ) is the slow damping amplitude and τ (cid:48) s is the oscil-lating period. With this ansatz solution, we can simplify thesquare root term in Eq. (C28) to (cid:113) − s − ( τ s ˙ s ) = (cid:115) − [ A cos ( t / τ (cid:48) s )] − τ s [ ˙ A cos ( t / τ (cid:48) s ) − A τ (cid:48) s sin ( t / τ (cid:48) s )] ≈ (cid:112) − A , (C36)where in the second line we use the fact that ˙ A / τ s (cid:28) τ s / τ (cid:48) s ∼ s + (cid:18) τ s + τ s τ φ α (cid:112) − A (cid:19) s ≈ − (cid:18) I τ θ τ s τ φ (cid:112) − A (cid:19) ˙ s . (C37)This equation is more precise than the simple damped har-monic oscillator Eq. (C30) since the square root is treatedwith a better approximation than rudely taken as unity. Nowwe try to obtain the analytical solution of this damped anhar-monic oscillator with appropriate approximation. We first cal-culate τ (cid:48) s by treating A as a constant within τ (cid:48) s and ignore thefriction term. These two approximations are valid because A varies much slower than τ (cid:48) s and the friction is ignorable in thetime scale of τ (cid:48) s . Then we can obtain a harmonic oscillatingequation, ¨ s + (cid:18) τ s + τ s τ φ α (cid:112) − A (cid:19) s ≈ , (C38)4 (a) (e) (b) (f) (c) (g) (d) (h) FIG. 7. (Color online) Comparison between numerical results and the analytical solutions shown in Eq. (C33) from the equations of dampedharmonic oscillators. (a) The numerical results with the same parameters as Fig. 2a but taken from the time range of t = [3830,6330] in theoriginal figure. The origin of time is shift to zero for comparison. (b) The numerical results with the same parameters as in (a) except for I = I c = I c , with data taken between the time range [250,10250]. (c) The numerical results with the same parameters as in (a) except for R = h / ( e ) , with data taken between the time range [7000,12000]. (d) The numerical results with the same parameters as in (a) except for δ = . E M , with data taken between the time range [650,1900]. (e-h) The analytic solution shown in Eq. (C33) with parameters taken thesame as in (a-d) respectively. The time origins in each figure have been shifted accordingly to make the initial value comparable to numericalresults. which gives the oscillating period as τ (cid:48) s = τ s (cid:113) + τ s τ φ α √ − A . (C39)Comparing with the oscillating period ˜ τ s obtained from thedamped harmonic oscillator approximation, the new oscillat-ing period τ (cid:48) s depends on the oscillating amplitude A . When A increases, the oscillating period τ (cid:48) s becomes larger. Thisagrees with the numerical results shown in Fig. 2a, and alsoagrees with the analysis based on the classical Hamiltonianwhich states that the oscillating frequency becomes largerwhen approaching the separatrix.Now we calculate the slow damping amplitude A ( t ) byplugging the ansatz solution back to the equation,¨ A cos ( t / τ (cid:48) s ) − A τ (cid:48) s sin ( t / τ (cid:48) s ) (C40) = − (cid:18) I τ θ τ s τ φ (cid:112) − A (cid:19) (cid:18) ˙ A cos ( t / τ (cid:48) s ) − A τ (cid:48) s sin ( t / τ (cid:48) s ) (cid:19) . Noticing the fact that ¨ A ( τ (cid:48) s ) (cid:28) ˙ A τ (cid:48) s (cid:28) A for slow varying A ,we obtain the equation for A asd A d t ≈ − (cid:18) I τ θ τ s τ φ (cid:112) − A (cid:19) A = − τ d A (cid:112) − A , (C41) which has the solution A ( t ) = e − t / τ d + e − t / τ d . (C42)Comparing with the result of damped harmonic oscillator ap-proximation in the previous section, this damping function ismore flat when s becomes large.Putting the expression for A ( t ) and τ (cid:48) s together, we finallyarrive at the analytical solution for s ( t ) inside the separatrix ofthe phase-space, s ( t ) = e − t / τ d + e − t / τ d cos ( t / τ (cid:48) s ) . (C43)Clearly this solution reduces back to Eq. (C33) when s ap-proaches zero. However, it provides a better result for thelarger s regime which captures two more details of the dampedoscillation shown in Fig. 2a. First, the oscillating period islarger when s is larger, which also agrees with our analysisbased on the action-angle formalism. Second, the dampingof the oscillating amplitude is slower at larger s , which isdifferent from the pure exponential decay which has a time-independent decay rate. We show the comparison of the nu-merical results with the analytical result Eq. (C43) in Fig. 8,and find better agreement with numerical results than the sim-ple damped harmonic oscillator approximation.5 (a) (e) (b) (f) (c) (g) (d) (h) FIG. 8. (Color online) Comparison between numerical results and the analytical solutions shown in Eq. (C43). (a-d) are the same as the (a-d)in Fig. 7. (e-h) The analytic solution shown in Eq. (C43) with parameters taken the same as in (a-d) respectively. The time origins in eachfigure are shifted accordingly to make the initial value comparable to numerical results.
5. Krylov-Bogoliubov method of averaging
Finally, we take an alternative method, the Krylov-Bogoliubov averaging method , to calculate the dampingfunction A ( t ) , and show identical results as from the dampedanharmonic oscillator approximation. We examine the Eq.(C28) again and make it dimensionless as,d s d τ + s = − τ s τ φ ( α s + βτ s d s d τ ) (cid:114) − s − ( d s d τ ) , (C44)where we define a dimensionless time as τ = t / τ s . Setting theright hand side of the equation to be zero we obtain,d s d τ + s = . (C45)This equation has the general solution of the form s = A (cid:48) cos ( τ + B ) , (C46)˙ s = − A (cid:48) sin ( τ + B ) . Now we recover the right hand side of the equation, andtake an anartz solution with the same trigonometric functionswhere A (cid:48) and B (cid:48) become time dependent, s = A (cid:48) ( τ ) cos ( τ + B (cid:48) ( τ )) (C47a)˙ s = − A (cid:48) ( τ ) sin ( τ + B (cid:48) ( τ )) . (C47b)In the following, we solve Eq. (C44) with these ansatz func-tions. We first take time derivative to Eq. (C47a) and obtain,˙ s = − A (cid:48) sin ( τ + B (cid:48) ) + ˙ A (cid:48) cos ( τ + B (cid:48) ) − A (cid:48) sin ( τ + B (cid:48) ) ˙ B (cid:48) . (C48) This equation must be equivalent to Eq. (C47b) for a self-consistent ansatz function, and thus we have a constraint equa-tion ˙ A (cid:48) cos ( τ + B (cid:48) ) = A (cid:48) sin ( τ + B (cid:48) ) ˙ B (cid:48) . (C49)We then plug the ansatz functions Eq. (C47) back to the orig-inal equation (C44) and obtain another constraint equation, − ˙ A (cid:48) sin ( τ + B (cid:48) ) − A (cid:48) cos ( τ + B (cid:48) ) ˙ B (cid:48) = E M δ ( α A (cid:48) cos ( τ + B (cid:48) ) + I τ θ A (cid:48) sin ( τ + B (cid:48) )) ∗ (cid:113) − A (cid:48) cos ( τ + B (cid:48) ) − A (cid:48) sin ( τ + B (cid:48) )= E M A (cid:48) δ ( α cos ( τ + B (cid:48) ) + I τ θ sin ( τ + B (cid:48) )) (cid:112) − A (cid:48) . (C50)Combining these two constraint equations (C49) and (C50),we arrive at equations for A (cid:48) ( τ ) and B (cid:48) ( τ ) as,dd τ (cid:18) A (cid:48) B (cid:48) (cid:19) = − E M δ (cid:0) α cos ( τ + B (cid:48) ) + I τ θ sin ( τ + B (cid:48) ) (cid:1) ∗ (cid:112) − A (cid:48) (cid:18) A (cid:48) sin ( τ + B (cid:48) ) cos ( τ + B (cid:48) ) (cid:19) . (C51)We note that no approximation has been made yet. The ansatzfunctions Eq. (C47) together with the constraint equationsEq. (C51) give an exact solution to the Eq. (C44). Nowwe concentrate on the slow varying part of A (cid:48) , denoting as A , which captures the slow damping of the oscillation. Wereplace cos ( τ + B (cid:48) ) and sin ( τ + B (cid:48) ) with their average values6within one period and obtaind A d τ = − E M πδ A (cid:112) − A ∗ (cid:90) π d τ (cid:20) α cos ( τ + B ) + I τ θ τ s sin ( τ + B ) (cid:21) sin ( τ + B )= − I τ θ τ φ A (cid:112) − A . (C52)Transforming back to real time with τ = t / τ s and rearrangingthe parameters, we simplify the equation of A to the formd A d t = − A τ d (cid:112) − A , (C53)which is exactly the same as we obtained from the dampedanharmonic oscillator approximation. Appendix D: Hysteresis with external parity flipping
Here we show that the hysteresis in the I-V curve still exitseven if the total parity of Majorana modes is broken by ex-ternal quantum levels from a single quasiparticle or impurity.For a model study, we consider the simplest case of an extraquantum level with a Hamiltonian of H i = ε d † d , (D1)where ε is the energy of the level which is near zero, and d † isthe creation operator on the level. This quantum level coupleswith one Majorana mode through the tunneling Hamiltonian, H T = T γ L d + T ∗ d † γ L = ( f †1 + f )( T d − T ∗ d † ) , (D2)where T is the tunneling strength. After including this quan-tum level, the Hilbert space is expanded and the total Hamil-tonian is an eight-by-eight matrix. It is also block diagonalwith two four-by-four blocks due to the conservation of thetotal parity. We can take one block by picking the basis statesas, d † | (cid:105) , d † f †1 f †2 | (cid:105) , f †2 | (cid:105) , f †1 | (cid:105) . Then we arrive at aneffective Hamiltonian H = ε + E M cos ( θ / ) δ L + δ R T ∗ δ L + δ R ε − E M cos ( θ / ) T ∗ T E M cos ( θ / ) − δ L + δ R T − δ L + δ R − E M cos ( θ / ) . (D3)The quantum average for the supercurrent through the Majo-rana channel is given by (cid:104) ψ | i γ γ | ψ (cid:105) = | ψ ( t ) | − | ψ ( t ) | + | ψ ( t ) | − | ψ ( t ) | . (D4)We plug the Eqs. (D3) and (D4) into the QRSJ model, and nu-merically obtain the I-V curve of the junction as demonstratedin Fig. 9. Clearly, the hysteresis behavior is insensitive to theparity flipping from the external quantum level. (a) (b) (c) FIG. 9. (Color online) Numerical simulation of the I-V curves withthe energies of the external quantum level as (a) ε =
0, (b) ε / E M = .
5, (c) ε / E M = − .
5. Parameters are taken as T / E M = . δ L / E M = .
005 and δ R / E M = .
015 and other parameters are takenthe same as Fig. 2a.
The reason that the parity flipping does not change thehysteresis is that the Hamiltonians for the odd total Majo-rana parity (the left-up 2x2 block) and the even total Majo-rana parity (the right-down 2x2 block) are qualitatively simi-lar. They both have avoided crossings at the Josephson phase θ = ( n + ) π . Naturally, we would expect that the quantumdynamics within each block is qualitatively the same, present-ing a damped oscillation. The small flipping energy T will notchange this quantum dynamics, therefore will not change thehysteresis behavior. Appendix E: Quasiparticle poisoning
In the topological superconductors, the quasiparticle poi-soning is an important obstacle for many signatures of Ma-jorana modes. The difference between the quasiparticle poi-soning and a simple external quantum level from impurity orquantum dot is that the quasiparticle poisoning comes fromthe thermal equilibrium fermionic environment which bringsdecoherence into the quantum two-level system defined byMajorana modes. This decoherence is fundamental fromthe quantum mechanical point of view, and cannot be sim-ply equivalenced to an enlarged Hilbert space. Then it isa natural question whether the decoherence from the quasi-particle poisoning will destroy the LZ effect induced hys-teresis. We analyze this problem by considering the den-sity matrix ρ ( t ) = ρ ( t ) | (cid:105)(cid:104) | + ρ ( t ) | (cid:105)(cid:104) | + ρ ( t ) | (cid:105)(cid:104) | + ρ ( t ) | (cid:105)(cid:104) | for the two-level system where the decoherencecan be naturally included using the Lindblad form. The dy-namics of the two-level system is then described by a masterequation , d ρ d t = − i ¯ h [ H , ρ ] + ∑ i τ i L i , (E1)where L i are all possible Lindblad forms which describe thedecoherence and τ i are the corresponding decoherence times.For a general two-level system, there are only three possi-ble Lindblad forms L = | ψ e (cid:105)(cid:104) ψ g | , L = | ψ g (cid:105)(cid:104) ψ e | , and L = | ψ e (cid:105)(cid:104) ψ e | − | ψ g (cid:105)(cid:104) ψ g | , where | ψ e (cid:105) and | ψ g (cid:105) are the two instan-taneous eigenstates of the two-level system. When consid-ering the decoherence from the quasiparticle poisoning, only7 (a) (b) (c) (d) FIG. 10. (Color online) Numerical results for the I-V curve withthe decoherence time (a) τ = h / E M , (b) τ = h / E M , (c) τ = . h / E M , and (d) τ = h / E M and τ = . h / E M . Otherparameters are taken the same as Fig. 2a. the relaxation processes described by L and the dephasingprocesses described by L are relevant in the low temperaturelimit.Let us first consider the relaxation processes given by theLindblad L , which involves the coupling between the Majo-rana modes and the quasiparticle states above the supercon-ducting gap. The decoherence time for this process is an ex-ponential function of the superconducting gap ,1 τ = λ Te − ∆ / T , (E2)where λ is a dimensionless factor estimated around 0 . .When the temperature is far below the superconducting gap T (cid:28) ∆ , the relaxation time is exponentially protected by thesuperconducting gap and would be quite long compared withall other time scales in the system. We present the results ofthe I-V curve with two different relaxation times of in Figs.10a and 10b. We see that a reasonable long relaxation timehas little influence on the hysteresis, while an extremely shortrelaxation time reduces the hysteresis but still does not changethe qualitative feature.We then consider the decoherence from the dephasing givenby the Lindblad L . Different from the relaxation, the dephas-ing should have a relatively short dephasing time with τ (cid:28) τ . However, looking at the form of L we see thatthe dephasing only introduces a decoherence in the relativephase of the two eigenstates, leaving the relative amplitudeunchanged. Since only the amplitude of the wave functionenters the dynamical equation for the Josephson phase in theQRSJ model, we would expect that the dephasing has little in-fluence on the hysteresis. We present the I-V curve for a veryshort dephasing time in Fig. 10c, and find that it indeed has no (a) (b) FIG. 11. (Color online) Numerical results of the I-V curves for theunderdamped junctions with (a) I c2 = I c2 = I c1 . The ca-pacitance is taken as C = . e / ¯ hI c1 . Other parameters are the sameas Fig. 2a. influence on the hysteresis behavior. Finally, we show the re-sult with a combination of the relaxation and dephasing in Fig.10d, and find that the hysteresis is robust to the decoherencefrom the quasiparticle poisoning. Appendix F: underdamped junction
The conventional Josephson junctions with negligible ca-pacitance show no hysteresis, making the LZ effect inducedhysteresis a novel phenomenon. However, even in the under-damped junctions where hysteresis is already expected fromthe shunted capacitance, the LZ effect still contribute a signifi-cant feature which might be useful for experimental detection.Here, we demonstrate a comparison between the I-V curvesof conventional and topological junctions in the underdampedregime, where the capacitance is included and the resistivelyshunted junction equation is rewritten as the resistively andcapacitively shunted junction equation. We show the numeri-cal results in Fig. 11. There is a hysteresis in the topologicaltrivial junction as expected from the standard theory, however,the difference between the switching and retrapping current islargely enhanced by the LZ effect induced part. Therefore,it is still a useful signal for detecting the Majorana modes inpossible topological junctions. ∗ These authors contributed equally to this work. † [email protected] A. Y. Kitaev, Phys. Usp. , 131 (2001). H. J. Kwon, K. Sengupta, and V. M. Yakovenko, Eur. Phys. J. B , 349 (2003). A. Kitaev, AIP Conference Proceedings , 22 (2009). X. L. Qi and S. C. Zhang, Rev. Mod. Phys. , 1057 (2011). L. Fu and C. L. Kane, Phys. Rev. Lett. , 096407 (2008). M. Sato, Y. Takahashi, and S. Fujimoto Phys. Rev. Lett. ,020401 (2009). Y. Tanaka, T. Yokoyama, and N. Nagaosa, Phys. Rev. Lett. ,107002 (2009). J. D. Sau, R. M. Lutchyn, S. Tewari, and S. Das Sarma, Phys. Rev.Lett. , 040502 (2010). J. Alicea, Rep. Prog. Phys , 076501 (2012). C. W. J. Beenakker, Annu. Rev. Con. Mat. Phys. , 113 (2013) S. R. Elliott and M. Franz, Rev. Mod. Phys. , 137 (2015). D. Aasen, M. Hell, R. V. Mishmash, A. Higginbotham, J. Danon,M. Leijnse, T. S. Jespersen, J. A. Folk, C. M. Marcus, K. Flens-berg, and J. Alicea, Phys. Rev. X , 031016 (2016). R. Aguado, Riv. Nuovo Cimento , 523 (2017). X. L. Qi, Taylor L. Hughes, S. Raghu, and S. C. Zhang, Phys.Rev. Lett. , 187001 (2009). T. H. Hsieh, G. B. Hal´asz, and T. Grover, Phys. Rev. Lett. ,166802 (2016). Z. Huang, S. Shimasaki, and M. Nitta, Phys. Rev. B ,220504(R) (2017). J. Nilsson, A. R. Akhmerov, and C. W. J. Beenakker, Phys. Rev.Lett. , 120403 (2008). M. Cheng, R. M. Lutchyn, V. Galitski, and S. Das Sarma, Phys.Rev. Lett. , 107001 (2009). T. Mizushima and K. Machida, Phys. Rev. A , 023624 (2010). T. D. Stanescu and S. Tewari, J. Phys.: Condens. Matter ,233201 (2013). S. M. Albrecht, A. P. Higginbotham, M. Madsen, F. Kuemmeth,T. S. Jespersen, J. Nyg˚ard, P. Krogstrup, and C. M. Marcus, Na-ture , 206 (2016). J. Cayao, P. San-Jose, A. M. Black-Schaffer, R. Aguado, and E.Prada, Phys. Rev. B , 205425 (2017). L. Fu and C. L. Kane, Phys. Rev. B , 161408 (2009). R. M. Lutchyn, J. D. Sau, and S. Das Sarma, Phys. Rev. Lett. ,077001 (2010). Y. Oreg, G. Refael, and F. von Oppen, Phys. Rev. Lett. ,177002 (2010). V. Mourik, K. Zuo, S. M. Frolov, S. R. Plissard, E. P. A. M.Bakkers, and L. P. Kouwenhoven, Science , 1003 (2012). M. T. Deng, C. L. Yu, G. Y. Huang, M. Larsson, P. Caro, and H.Q. Xu, Nano Lett. , 6414 (2012). F. Dom´ınguez, F. Hassler, G. Platero Phys. Rev. B , 140503(2012). P. San-Jose, E. Prada, and R. Aguado, Phys. Rev. Lett. ,257001 (2012). L. Allen, and J. H. Eberly (1974),
Optical Resonance and Two-level Atoms (Dover, 1975). L. M. K. Vandersypen and I. L. Chuang, Rev. Mod. Phys. ,1037 (2005). O. Morsch and M. Oberthaler, Rev. Mod. Phys. , 179 (2006). S. N. Shevchenko, S. Ashhab and F. Nori, Phys. Rep. , 1(2010). L. Landau, Phys. Z. Sowjetunion, , 46 (1932); C. Zener, Proc. R.Soc. Landon Ser. A , 369 (1932); E. Majorana, Nuovo Cimento 9, 43(1932). W. C. Huang, Q. F. Liang, D. X. Yao and Z. Wang, Phys. Rev A , 012308 (2015). D. Averin and A. Bardas, Phys. Rev. Lett. , 1831 (1995); L.Y. Gorelik, N. I. Lundin, V. S. Shumeiko, R. I. Shekhter, and M.Jonson, Phys. Rev. Lett. , 2538 (1998). B. Wu and Q. Niu, Phys. Rev. A , 023402 (2000). Y. A. Chen, S. D. Huber, S. Trotzky, I. Bloch, and E. Altman,Nature Physics , 61 (2011). X. J. Liu, K. T. Law, T. K. Ng, and P. A. Lee, Phys. Rev. Lett. , 120402 (2013). F. Forster, G. Petersen, S. Manus, P. Hanggi, D. Schuh, W.Wegscheider, S. Kohler, and S. Ludwig, Phys. Rev. Lett. ,116803 (2014). W. Y. He, S. Z. Zhang, and K. T. Law, Phys. Rev. A , 013606(2016). T. Higuchi, C. Heide, K. Ullmann, H.B. Weber, and P. Hommel-hoff, Nature , 224 (2017). J. B. Oostinga, L. Maier, P. Sch¨uffelgen, D. Knott, C. Ames, C.Br¨une, G. Tkachov, H. Buhmann, and L. W. Molenkamp, Phys.Rev. X , 021007 (2013) J. Wiedenmann, R. S. Deacon, S. Hartinger, O. Herrmann, T. M.Klapwijk, L. Maier, C. Ames, C. Br¨une, C. Gould, A. Oiwa, K.Ishibashi, S. Tarucha, H. Buhmann, L. W. Molenkamp, and E.Bocquillon, Nature Communications , 10303 (2016). V. S. Pribiag, A. J. A. Beukman, F. Qu, M. C. Cassidy, C. Char-pentier, W. Wegscheider, and L. P. Kouwenhoven, Nature Nan-otechnology , 593 (2015). S. Nadj-Perge, I. K. Drozdov, J. Li, H. Chen, S. Jeon, J. Seo, A.H. MacDonald, B. A. Bernevig, and A. Yazdani, Science , 602(2014). M. Tinkham,
Introduction to Superconductivity , (Second Edition,McGraw-Hill Book Co. 1996). J. A. Blackburn, M. Cirillo, N. Grønbech-Jensen, Phys. Rep. ,1 (2016). D. Dragoman and M. Dragoman,
Quantum-Classical Analogies (Springer, 2004). H. Fu, Z. C. Gong, T. H. Mao, C. P. Sun, S. Yi, Y. Li, and G. Y.Cao, Phys. Rev. A , 043855 (2016). J. Liu, L. Fu, B. Y. Ou, S. G. Chen, D. I. Choi, B. Wu, and Q. Niu,Phys. Rev. A , 023404 (2002). J. Liu, B. Wu, and Q. Niu, Phys. Rev. Lett. , 170404 (2003). D. R. Smith,
Singular-Perturbation Theory , (Cambridge Univer-sity Press, 1985). L. Pitaevskii, S. Stringari,
Bose-Einstein Condensation , (Claren-don Press, 2003). W. Dittrich and M. Reuter,
Classical and Quantum Dynamics ,(Springer-Verlag, 1994). Y. Peng, F. Pientka, E. Berg, Y. Oreg, and F. von Oppen, Phys.Rev. B , 085409 (2016). D. Rainis and D. Loss, Phys. Rev. B , 174533 (2012). M. J. Schmidt, D. Rainis, and D. Loss, Phys. Rev. B , 085414(2012). J. Wojewoda, A. Stefanski, M. Wiercigroch, and T. Kapitaniak,Phil. Trans. R. Soc. A , 747 (2008). S. Peotta and M. DiVentra, Phys. Rev. Applied , 034011 (2014). S. N. Shevchenko, Y. V. Pershin, and F. Nori, Phys. Rev. Applied , 014006 (2016). C. Guarcello, P. Solinas, M. D. Ventra, and F. Giazotto, Sci. Rep. , 46736 (2017). J. A. Sanders, F. Verhulst, and J. Murdock,
Averaging Methods inNonlinear Dynamical Systems , (Springer, 2007). B. van Heck, F. Hassler, A. R. Akhmerov, and C. W. J. Beenakker,Phys. Rev. B , 180502 (2011). M. Veldhorst, C. G. Molenaar, C. J. M. Verwijs, H. Hilgenkamp,and A. Brinkman, Phys. Rev. B , 024509 (2012). N. M. Krylov and N. N. Bogolyubov,