Joule overheating poisons the fractional ac Josephson effect in topological Josephson junctions

Topological Josephson junctions designed on the surface of a 3D-topological insulator (TI) harbor Majorana bound states (MBS's) among a continuum of conventional Andreev bound states. The distinct feature of these MBS's lies in the $4\pi$-periodicity of their energy-phase relation that yields a fractional ac Josephson effect and a suppression of odd Shapiro steps under $r\!f$ irradiation. Yet, recent experiments showed that a few, or only the first, odd Shapiro steps are missing, casting doubts on the interpretation. Here, we show that Josephson junctions tailored on the large bandgap 3D TI Bi$_2$Se$_3$ exhibit a fractional ac Josephson effect acting on the first Shapiro step only. With a modified resistively shunted junction model, we demonstrate that the resilience of higher order odd Shapiro steps can be accounted for by thermal poisoning driven by Joule overheating. Furthermore, we uncover a residual supercurrent at the nodes between Shapiro lobes, which provides a direct and novel signature of the current carried by the MBS. Our findings showcase the crucial role of thermal effects in topological Josephson junctions and lend support to the Majorana origin of the partial suppression of odd Shapiro steps.

Topological Josephson junctions designed on the surface of a 3D-topological insulator (TI) harbor Majorana bound states (MBS's) among a continuum of conventional Andreev bound states. The distinct feature of these MBS's lies in the 4π-periodicity of their energy-phase relation that yields a fractional ac Josephson effect and a suppression of odd Shapiro steps under rf irradiation. Yet, recent experiments showed that a few, or only the first, odd Shapiro steps are missing, casting doubts on the interpretation. Here, we show that Josephson junctions tailored on the large bandgap 3D TI Bi2Se3 exhibit a fractional ac Josephson effect acting on the first Shapiro step only. With a modified resistively shunted junction model, we demonstrate that the resilience of higher order odd Shapiro steps can be accounted for by thermal poisoning driven by Joule overheating. Furthermore, we uncover a residual supercurrent at the nodes between Shapiro lobes, which provides a direct and novel signature of the current carried by the MBS's. Our findings showcase the crucial role of thermal effects in topological Josephson junctions and lend support to the Majorana origin of the partial suppression of odd Shapiro steps.
Topological superconductivity engineered by coupling superconducting electrodes to topological states of matter has attracted considerable attention due to the prospect of manipulating Majorana states for topological quantum computing [1][2][3].
Another key approach to substantiate the very existence of MBS relies on the fractional ac Josephson effect [13][14][15] that develops in topological Josephson junction [3,15]. Theory predicts that MBS's shall emerge in such junctions as a peculiar, spinless Andreev bound state (ABS). Contrary to the conventional ABS's whose energy level varies 2π-periodically with the phase difference φ between the junction electrodes, the MBS is 4πperiodic and crosses zero-energy for a phase π (see Fig.  1a), yielding a fractional ac Josephson effect at frequency f J /2 = eV /h (e is the electron charge, V the voltage drop across the junction and h the Planck constant), that is, half the Josephson frequency f J [3,15].
Yet, revealing such a 4π-periodic contribution has proven challenging in dc transport experiments due to the presence of often prevailing, conventional ABS's [16][17][18]. Moreover, poisoning processes -stochastic paritychanges of the quasiparticle occupation number -may obscure the MBS contribution by limiting its lifetime [19,20]. Measurement schemes probing at timescales shorter than this lifetime are thus essential. The Shapiro effect comes forth with the combined advantages of a radiofrequency (f rf ) excitation of the phase that can be faster than the poisoning dynamics [15,21,22], and the ease of dc current-voltage (IV 's) characteristics measurements.
The immediate consequence of the fractional ac Josephson effect is an unusual sequence of Shapiro voltage steps ∆V = hf rf e in the IV characteristics, twice of that of conventional Shapiro steps ( hf rf 2e ) [13,[21][22][23][24][25][26], providing direct evidence for the MBS 4π-periodicity. First experiments performed on InSb nanowires [27], on strained HgTe 3D TI [28], and on Bi 1−x Sb x [29] junctions however showed surprises in the sequences of Shapiro steps. In all these cases, only the n = ±1 steps were absent in a given range of rf power and frequency (n is the integer index of the Shapiro steps), an absence which was described as an (incomplete) signature of the fractional a.c. Josephson effect. More recent measurements on 2D HgTe quantum wells showed the absence of odd steps up to n = 9 [30] and Josephson radiation at half the Josephson frequency [31], though without the demonstration of time-reversal symmetry breaking that is required to induce MBS's in quantum spin-Hall edge channels [15,32]. While the latter observations advocate more clearly for the existence of 4π-periodic Andreev modes, the fact that the fractional ac Josephson effect acts only on some odd Shapiro steps depending on the system remains unclear. Whether it provides a signature of the Majorana mode is a central question for identifying topological superconductivity in a variety of implementations.
In this work we report on the observation and understanding of the partial fractional ac Josephson effect in Josephson junctions designed on exfoliated flakes of the 3D topological insulator Bi 2 Se 3 . Our devices exhibit an anomalous sequence of Shapiro steps with (only) the first step absent at low rf power and frequency. The 4π-periodic contribution to the supercurrent is directly identified as a residual supercurrent at the first node of the critical current when the rf power is increased.
To shed light on our findings, we develop a two-channel Resistively Shunted Junction (RSJ) model that includes the quasiparticle overheating induced by Joule effect [33][34][35], and a thermally activated poisoning of the MBS. We show that Joule overheating suppresses the parity lifetime of the MBS and thus terminates the 4π-periodic contribution to any higher index Shapiro steps, accounting for the observed suppression of the first Shapiro step only.

Results
Partial even-odd effect in Bi 2 Se 3 Josephson junctions. Our samples are based on flakes of Bi 2 Se 3 crystals exfoliated with the scotch tape technique. Figure 1c shows a 30 nm thick flake of Bi 2 Se 3 contacted with multiple electrodes of vanadium enabling both magnetotransport and Josephson junction measurements on the same Bi 2 Se 3 crystal (see Appendix for fabrication and measurement details). Analysis of Shubnikov-de-Haas oscillations and Hall effect enable identification of three electronic populations contributing to the sample conductance: Bulk states with a charge carrier density of 4.5 × 10 19 cm −3 , and the top and bottom surface states with densities of 1 × 10 12 and 4 × 10 12 cm −2 respectively (see Supplementary Note 1). All three channels may thus carry supercurrent by proximity effect [36].
We focus here our discussion on the Josephson junction of length L = 125 nm and width W = 2.25 µm (see geometry in Fig. 1b) indicated by the white arrow in Fig. 1c. Below the superconducting transition temperature of the electrodes (T c = 5 K), the proximity effect develops in the TI, leading to a dissipationless supercurrent in the IV 's as shown in Figure 1d. The transition to the resistive state of the junction (R = 7.5 Ω) is hysteretic at 0.05 K with switching and retrapping currents of I sw = 7.3 µA and I r = 5.0 µA, respectively. Such a hysteresis is a common feature of most Josephson junctions made with metallic weak links and results from a quasiparticle overheating in the normal section of the junction [34]. As we will show below, the ensuing quasiparticle overheating is key for understanding the suppression of the n = ±1 Shapiro steps only.
The dc response of the Josephson junction to an rf irradiation is shown in the Shapiro map of Figure 2a that displays the color-coded differential resistance dV /dI versus rf power P rf and dc current I for an rf frequency of 3.5 GHz. At this frequency, well-defined Shapiro steps develop in the IV curves, two of which are shown in Fig. 2c, with voltage steps that match the standard value V n = ±n hf rf 2e = n × 7.2 µV expected for a 2π-periodic current-phase relation. In the Shapiro map, the black areas indicate dV /dI = 0 and hence the position and amplitude of the Shapiro steps in the P rf − I plane. Two features standard for a usual Shapiro map are visible. Energy-phase spectrum of the Andreev bound states for a topological Josephson junction at the surface of a 3D TI. 4π-periodic, spinless, MBS's coexist with conventional 2πperiodic ABS's. The maximum energy EABS is lower than the quasiparticle continuum at ∆ in case of imperfect interface transparency or the presence of a magnetic barrier [3,15,37]. b, Schematic of the Josephson junction geometry showing the superconducting electrodes (S) in orange that contact the top of the Bi2Se3 flake (in blue). c, Optical image of device. Scale bar is 9 µm. d, Current-voltage characteristic of the junction indicated by the arrow in c. Measurements were carried out at 0.05 K.
First, on increasing P rf , the critical current continuously decreases till nearly full suppression at P rf = 9 dBm and then oscillates at higher P rf . Second, the sequence of appearance of the Shapiro steps with P rf is sorted by the Shapiro step index n, and, importantly, starts with the step n = 1.
The central experimental result of this work is displayed in Figure 2b, where we show the rf response of the same junction at a lower frequency of 1 GHz. This Shapiro map exhibits distinct features that markedly differentiate it from the higher frequency map. On increasing P rf , the first Shapiro step n = 1 sets in at high P rf after the steps of higher index for both switching and retrapping currents. This anomaly, sometimes termed even-odd effect [21,22], results in the conspicuous absence of the first Shapiro step in IV 's picked out at low P rf , while steps of higher indexes already appears, see Fig. 2d. Our findings match those recently obtained on InSb nanowires [27], strained HgTe 3D TI [28] and Bi 1−x Sb x alloy [29], which were interpreted as a signature of a 4π-periodic MBS contribution.
A second and new feature emerges at the first minimum of the critical current when the rf power is increased, i.e. at P rf 10.7 dBm. Contrary to a con- ventional Shapiro map where a complete supercurrent suppression is expected at what can be called a resistive node, a small supercurrent I c ∼ 190 nA remains. This appears clearly in the individual IV 's of Fig. 2d, see for instance the blue curve there. Figure 2e displays a similar Shapiro map obtained at a slightly different f rf in the low-frequency regime, but zoomed on the resistive node where the critical current is expected to vanish but does not. We shall see in the following that this residual supercurrent provides a direct signature of the presence of a 4π-periodic mode in the junction.
Determination of the coherent transport regimes. Capturing the ABS spectrum of a Josephson junction on 3D TI's remains difficult as several conduction channels, including bulk and surfaces, intervene. Should all channels carry supercurrent, the nature of charge carriers in them may lead to virtually different regimes of coherent transport that we assess in the following. For bulk carriers, a rough estimate of the mean free path (see Supplementary Note 1) gives a Thouless energy E th = D/L 2 417 µeV, smaller than the superconducting gap of the vanadium electrodes ∆ = 800 µeV. This channel thus belong to the class of long diffusive Josephson junctions [38]. In contrast, for the topological surface states, the spin texture of the Dirac elec-trons stemming from the spin-momentum locking leads to a very strong scattering anisotropy which promotes forward scattering. As a result, the transport time τ tr is expected to be significantly enhanced with respect to the elastic scattering time τ e , with ratio τ tr /τ e up to 60 depending on the disorder source [39]. Recent experiments on Bi 2 Se 3 flakes combining field effect mobility and quantum oscillations assessed a ratio τ tr /τ e 8 [40]. Taking the latter value as a conservative estimate and the surface state elastic mean-free path l e 28 nm of our sample (see Supplementary Note 1) leads to a transport length l tr 225 nm. These considerations suggest that surface transport is ballistic with, importantly, a non-zero probability for straight electronic trajectories impinging both electrodes. This is also consistent with signatures of ballistic transport over 300 nm evidenced in Bi 2 Se 3 nanowires [41,42].
Consequently, we consider the topological surface state channel as ballistic. As such, the relevant energy scale for the ABS's is v F /L = 2.8 meV with v F = 5.4 × 10 5 m/s the Fermi velocity [43]. It is greater than ∆, which shall lead to ABS's in the short ballistic limit. Theory then predicts that a 4π-periodic MBS exists, even in the case where the Fermi level is far from the Dirac point of the surface states, which is here the experimentally relevant regime [37]. This MBS corresponds to ballistic trajectories impinging perpendicularly the superconducting electrodes, all other incidence angles yielding conventional 2π-periodic ABS's [3,37].
Eventually, observability of the 4π-periodicity theoretically implies a strong constraint on the Andreev spectrum: The MBS's must be decoupled from the quasiparticle continuum at φ = 0 and 2π to avoid direct transfer of quasiparticles into or from the continuum. Such a transition would indeed occur from the excited to the ground state every 2π, restoring an effective 2π-periodicity for the MBS. This detrimental effect can be remedied by adding a magnetic layer or magnetic field that break time reversal symmetry, and thus open a gap between the MBS and the quasiparticle continuum [15,21,23,37,44], as sketched in Fig. 1b. In our samples, the vanadium that we use as superconducting electrodes is known to form magnetic dopants in Bi 2 Se 3 and eventually a ferromagnetic phase at large concentration [45]. Given that a smooth ion milling of the Bi 2 Se 3 surface is processed before vanadium deposition, favoring vanadium diffusion into the Bi 2 Se 3 crystal, there is presumably a magnetic layer or local magnetic moments at the superconducting interface as well as magnetic moments on the oxidized vanadium side surfaces of the electrodes. This singular configuration is likely to break time reversal symmetry on the scale of the junction, thus leading to the decoupling of the MBS from the continuum and to the ensuing observability of 4π-periodicity in our Shapiro maps.
Two-channel RSJ model. To understand our experimental findings, we consider a RSJ model comprising a pure Josephson junction in parallel with a shunt resistor R. In the usual scheme of a single Josephson channel with a critical current I c , the key parameter for the phase dynamics is the phase relaxation time τ J = 2eRIc which sets the typical time scale for the phase to adapt to a drive current change [33]. With a rf drive, the RSJ model thus acts on the phase as a low pass filter of cutoff frequency 1/τ J (see Supplementary Note 3). The regime of visibility of Shapiro steps is thus defined by f rf τ J < 1.
To phenomenologically capture the complex dynamics of a topological Josephson junction where a MBS lies within (a majority of) ABS's, we include two different Josephson junctions J1 and J2 in the RSJ model (see Fig. 3a) [24,26]. The first junction J1 stands for the conventional ABS's with a 2π-periodic current-phase relation, and the second one J2 represents the 4π-periodic MBS's. Having both 2π and 4π-periodic contributions in the total supercurrent I s (φ) = I 2π c sin(φ) + I 4π c sin(φ/2) drastically changes the Shapiro steps sequence. The dynamics is now ruled by two different phase relaxation and I 4π c . Consequently, the 4π-periodic contribution will impact the junction dynamics only for drive frequencies f rf < 1/τ 4π J . This can be straightforwardly seen in the Shapiro maps that we obtained by numerically solving the RSJ equation together with the Josephson relation V = 2e < dφ dt >. In the low frequency limit f rf τ 2π J < f rf τ 4π J < 1, all even Shapiro steps develop at rf power P rf lower than their neighboring odd steps, leading to the following appearance sequence |n| = {2, 1, 4, 3, 6, 5...} on increasing P rf (see Fig. 3c). Lowering f rf would enhance this evenodd effect and result ultimately in a quasi-suppression of odd Shapiro steps. Conversely, when the drive frequency is faster than 1/τ 4π J but still lower than 1/τ 2π J , that is, f rf τ 4π J > 1 > f rf τ 2π J , the 4π-periodic contribution is suppressed, restoring the regular sequence of Shapiro steps appearance |n| = {1, 2, 3, 4...} on increasing P rf , as shown in Fig. 3d. In a 2D space (f rf , P rf ), the ranges of existence of the even and odd Shapiro steps are sketched in Fig. 3b. Importantly, the even-odd effect is robust even if I 4π c sounds negligible compared to I 2π c , since the even-odd effect will always be present at low enough frequency as soon as f rf < 1/τ 4π J . This low-frequency observability of the even-odd effect furthermore excludes an explanation for the existence 4π-periodic contribution based on Landau-Zener transitions at a soft gap, which should be otherwise enhanced at high frequencies.
Let us now consider in the same computed maps the P rf -dependence of the Shapiro steps amplitude. In contrast with the standard oscillatory behavior with a complete suppression at the resistive nodes, the even steps, including the supercurrent branch (n = 0), exhibit a non-vanishing amplitude at every two resistive node on increasing P rf , see Fig. 3c and d. This unusual feature has been predicted in a recent paper by Domínguez et al. [26]. It is a direct consequence of the presence of the 4π-periodic channel. It also conspicuously matches the residual supercurrent at the first node of the supercurrent branch in the experimental data of Fig. 2b and e.
We demonstrate below that the amplitude of this residual supercurrent can be quantitatively related to the 4πperiodic critical current I 4π c . Figure 4a displays the computed switching current I 0 of the Shapiro step n = 0, extracted from Fig. 3c as a function of P rf . This plot both highlights the oscillatory behaviour of I 0 with P rf and enables us to identify the residual supercurrent of the first node that we note I k=1 0 , with k the node index. To demonstrate the correlation between I k=1 0 and I 4π c , we performed a numerical study of the dependence of I k=1 0 on the relevant parameters I 4π c and f rf by systematically computing Shapiro maps for different sets of parameters. Figure 4b displays the calculated I k=1 indicating that this residual supercurrent is visible at higher frequencies than the even-odd effect on the step appearance order.
The remarkable consequence of the collapse is that I 4π c is uniquely defined for a fixed set of parameters I k=1 0 , R and f rf . Using a polynomial fit of the data points in Fig. 4b (black dashed line), we can thus directly determine I 4π c from the measured value of the residual supercurrent I k=1 0 at a given frequency f rf (see Supplementary Note 5). With the experimental parameters of Fig. 2b, that is, R = 7.5 Ω, I k=1 0 = 190 nA and f rf = 1 GHz, we obtain I 4π c = 290 nA and thus a ratio I 4π c /I c ∼ 4%. This value of I 4π c is close to the theoretical value of the supercurrent carried by a single mode e∆/ = 195 nA with ∆ = 0.8 meV being the superconducting gap of the vanadium electrodes [46], although the exact value of the induced gap is most likely smaller than ∆ due to nonperfect interface transparency. Furthermore, the resulting values f rf τ 4π J 0.3 and 1 for the 1 GHz and 3.5 GHz Shapiro maps of Fig. 2b and a respectively are in agreement with the presence (absence) of even-odd effect in the data, confirming the consistency of the analysis. Note that our estimate of I 4π c is close those found on strained HgTe [28] and Bi 1−x Sb x [29] systems, which were based on the frequency cut-off criterion.
Joule-induced poisoning of the MBS. The twochannel RSJ model, however, does not account for the experimental absence of the first Shapiro step only. Inclusion of a junction capacitance (geometric or instrinsic [47]) through a two-channel RSCJ model mitigates the suppression of the n ≥ 3 odd steps [48], but is however not relevant for our strongly overdamped Bi 2 Se 3 -based Josephson junctions [49]. Furthermore, in the case of strongly underdamped, hysteresis should also be sizable in the Shapiro steps [50], which is not observed in our experiment. We propose instead that the absence of the first step can be explained by including thermal effects resulting from Joule heating, and their impact on the 4πperiodic mode. Some of us recently showed that electron overheating must be taken into account with the RSJ model to capture Shapiro maps in conventional metallic Josephson junctions [35]. Assuming that quasiparticles in the junction form a thermal distribution with an effective temperature T qp different from the phonon bath temperature T ph , we included the temperature dependence of the critical current, I c (T qp ), and solved the RSJ equation self-consistently together with the heat balance equation P =< I(t)V (t) >= ΣΩ(T 5 qp − T 5 ph ) (Σ is the electron-phonon coupling constant and Ω the volume of the normal part, see Supplementary Note 6 for values) to extract T qp for each dc current. Solving such a thermal RSJ model (tRSJ), we obtain a significant raise of T qp , when a dc voltage drop sets in on the firstly developed Shapiro step [35].
We then conjecture that the ensuing excess of nonequilibrium quasiparticles that we express in terms of an effective quasiparticle temperature poisons the 4πperiodic mode. Note that any non-thermal distribution would have the same consequences. Acting on a single mode, poisoning causes a stochastic paritychange [15,21,22,31,44] that switches in time the quasiparticle occupation from the excited to the ground state, as illustrated in Fig. 1b by the black dotted arrow. Within the two-channel RSJ model, we model this parity-change of the 4π-periodic contribution to the supercurrent as where n sw (t) is a random occupation number of characteristic timescale determined by a switching parity lifetime τ sw (see Supplementary Note 4). The exact microscopic processes that yield such a diabatic event can involve several transitions, including pair breaking, quasiparticle recombination with various rates, and possible coupling to the bulk states. Phenomenologically, but without loosing the generality of the foregoing, we follow the approach of Fu and Kane [15] and consider a thermally activated switching parity lifetime τ sw for the 4π-periodic mode: where ∆ − E M BS is the minimal energy gap separating the 4π-periodic mode to the continuum (see Fig. 1b), τ 0 the characteristic timescale for the activation process, and k B the Boltzmann constant. When the quasiparticle effective temperature T qp is high, the switching time is small τ sw τ J and poisoning suppresses the 4πperiodic Josephson effect, and hence the even-odd effect, by switching frequently the current-phase relation (see Supplementary Note 4). For the opposite limit τ sw τ J , poisoning is irrelevant.
We thus claim that understanding the partial even-odd effect in topological Josephson junctions relies on the interplay between activated poisoning and thermal effects in the two-channel tRSJ model. We show in Fig. 5a-b two Shapiro maps computed for the same rf frequency and junction parameters as those in Fig. 3c. Fig. 5a is the result of the tRSJ model that includes a T -dependence of I 2π c (T ) given by the measured I sw (T ) [35] and realistic parameters for the heat balance equation (see Supplementary Note 6 for an estimate of ΣΩ in our junctions). We assume here that most of the T -dependence of I sw (T ) comes from the 2π-periodic modes. Compared to Fig. 3c, the electron overheating leads to a broadening of the resistive transitions at high P rf between the periodic oscillations of the Shapiro steps. Nevertheless, the full even-odd effect acting on all odd steps remains as in the absence of thermal effect, see Fig. 3c. Figure 5b is the main result of our theoretical analysis. It displays a Shapiro map computed with the same tRSJ than Fig. 5a, but including thermally activated poisoning of the 4π-periodic channel defined by Eq. (1) and (2) which ensues from Joule overheating. The 4π-periodic channel now impacts only the steps n = ±2 by enhancing their amplitude, therefore inverting the appearance order on increasing P rf between step n = 1 and 2. The sequence of appearance of all higher order steps turns out to be regularized due to the suppression of the 4πperiodic contribution by poisoning. This finding, that is, the even-odd effect limited to the first Shapiro step only, is in full agreement with our experiment shown in Fig. 2b and with works on other systems [27][28][29].
The thermally activated poisoning can be captured by inspecting the computed T qp and τ sw for two different P rf 's. Figure 6a displays the IV curves corresponding to the black arrows in Fig. 5b. In Fig. 6b, we show the corresponding T qp versus I, which raises linearly once a dissipative voltage sets in. Accordingly, τ sw is exponentially suppressed, and becomes inferior to τ J on the Shapiro steps n ≥ 2 (Fig. 6c). This explains why the 4πperiodic component acts only on the appearance order of the n = 1 and 2 steps, leaving all other steps unaffected.
Furthermore, inspecting Fig. 5b we see that the residual supercurrent of only the first resistive node I k=1 0 remains in presence of poisoning, as observed in the data of Fig. 2b. This confirms that the residual supercurrent I k=1 0 is a robust feature that provides a new indicator for the presence of 4π-periodic modes, and enables a direct and quantitative determination of the corresponding 4π-periodic critical current as discussed above. Within this approach of thermal poisoning by Joule overheating, it is interesting to compare the power dissipated between the available experiments on different TI systems. For instance, estimates of the dissipated power on the n = 2 Shapiro step give two orders of magnitude difference between different TI materials: P 2 hf rf 2e I ∼ 25 pW for our data and 17 pW for strained HgTe [28]. For InAs nanowires [27] and Bi 1−x Sb x [29], P ∼ 2 pW. For HgTe quantum wells [30], the dissipated power is the smallest P ∼ 0.2 pW. Comparing the strained HgTe to the HgTe quantum wells that share the same electronphonon coupling constant, we estimate a power per unit of volume of ∼ 80 fW.µm −3 and ∼ 2 fW.µm −3 , respectively (∼ 300 fW.µm −3 in our Bi 2 Se 3 samples). Such a small dissipated power density in the experiment on HgTe quantum wells should result in a minimized amount of non-equilibrium quasiparticles and limited poisoning, therefore explaining the observed full suppression of not only the first but of several odd Shapiro steps in this system.

Discussion
Our theoretical approach, combining two Josephson channels in parallel, electronic overheating and quasiparticle poisoning, can be extended to more elaborate situations. For instance, recent theory works predict for the case of the quantum spin Hall regime an 8π-periodicicity due to either interactions [51] or quantum magnetic impurities [52,53]. Although this has not been reported so far in experiments, it would be interesting to study how multiple periodicities mix in the Shapiro response. Given the understanding of the two-channel RSJ model, we expect the 8π-periodicity to enhance every steps of index ±4n and significantly modify the beating pattern at high P rf . Other effects such as the voltage dependence of the phase relaxation time [21] could be included in our model and should enhance the effect of thermal poisoning on the partial suppression of the even-odd effect.
To conclude, our work elucidates the origin of the puzzling suppression of only the first Shapiro step in topological Josephson junctions. In our Bi 2 Se 3 Josephson junctions, this suppression is accompanied by a residual supercurrent that provides a new indicator of the 4π-periodic contribution to the supercurrent. Together, these observations can be captured by a two-channel thermal RSJ model in which Joule overheating activates poisoning of the 4π-periodic mode. The even-odd effect restricted to the first Shapiro step and the residual supercurrent do provide a clear signature of a 4π-periodic mode in the Andreev spectrum, conspicuously pointing to MBS's. Our phenomenological model illustrates a direct consequence of thermal poisoning on Majorana bound states, signaling that dissipation must be scrutinized with attention in dc-biased measurements. Addressing a microscopic description of the enhanced poisoning in such non-equilibrium measurement schemes is a challenging task for theory that should lead to significant progress towards new devices for Majorana physics and possible MBS qubits.

Methods
Bi 2 Se 3 crystals were synthesized by melting growth method with high purity (5N) Bi and Se in an evacuated quartz tube. Crystals were analyzed by Xrays diffraction, and angle-resolved photoemission spectroscopy. Flakes of Bi 2 Se 3 were exfoliated from the bulk crystal on silicon wafer and systematically inspected by atomic force microscopy to ensure crystal quality. V/Au superconducting electrodes were patterned by e-beam lithography and deposited by e-gun evaporation after a soft ion beam etching.
Measurements were performed in a dilution refrigerator equipped with highly filtered dc lines that comprise room temperature feed-through Pi-filters, lossy custommade coaxial cables and capacitors to ground on the sample holder. Radio-frequency are fed through a dedicated coaxial cable ending as an antenna that were adjusted in the vicinity of the devices. Shapiro map measurements were performed with standard lockin amplifier techniques.

I. Transport properties of the Bi2Se3 flakes
We present in this section the magnetotransport properties of sample LC106 discussed in the main text and determine the possible contributions to the conductance of the bulk and surface states. The magnetoresistance measured at 4 K and up to 14 T is shown in Figure S1a. Above 7 T, Shubnikov-de-Haas (SdH) oscillations develop on top of a linear magnetoresistance characteristic of 3D topological insulators [54,55]. The oscillations after background removal are shown in Figure S1c. The oscillatory pattern displays one obvious frequency, but also a beating signaling the presence of additional ones. A Fourier transform shown in Fig. S1d reveals three main peaks, at magnetic frequencies f B1 = 43 T, f B2 = 163 T, f B3 = 400 T (see Fig. S1d), similar to what was observed in previous works [40,56,57]. In Bi 2 Se 3 nanostructures, three electronic populations are expected to contribute to electronic transport: One bulk band and two topological surface states of the lower and upper surfaces [40,[56][57][58], with respective charge carrier densities n b , n SS1 and n SS2 . It is worth noticing that for such strong charge carrier densities no downward band bending, and consequently no trivial charge accumulation 2DEG are expected [57][58][59]. The total carrier density n tot is therefore the sum of these three carrier densities. To identify the magnetic frequencies to the different populations, we use the total carrier density extracted from Hall measurements. Figure S1b displays the transverse resistance measured in sample LC106. A linear fit yields a total carrier density of n Hall tot = 1.5 × 10 14 cm −2 indicating a high electronic doping. From the different possible scenarios for the SdH frequencies repartition, only one is compatible with the Hall data. By assuming either bulk or surface origin for the different magnetic frequencies, we can calculate a total carrier density n SdH tot = n b × d + n SS1 + n SS2 , with d = 30nm the flake thickness. Attributing either f B1 or f B2 to the bulk contribution (and therefore, the other two to TSS) would result in a total density of n SdH tot = 1.8 × 10 13 cm −2 or 4.6 × 10 13 cm −2 , which is much too low compared to the total carrier density n Hall tot . On the contrary, if we associate f B3 with bulk carriers, we obtain a total density n SdH tot = 1.4 × 10 14 cm −2 n Hall tot . We therefore conclude that f B3 is related to the bulk band, leading to a bulk carrier density of n b = 4.5 × 10 19 cm −3 , while f B1 and f B2 correspond to the topological surface states with n SS1 1×10 12 cm −2 and n SS2 4×10 12 cm −2 , respectively. Given that the peak at f B2 gives the largest contribution to the SdH oscillations, we attribute f B2 to the top surface and f B1 to the bottom surface which has most likely a lower mobility than the top one due to its proximity to the SiO 2 substrate.
With the above charge carrier densities, we evaluate the number of transport modes for each conduction channel: N = W k F π 130 and 250 for the bottom and top surface states respectively. For the bulk band, one obtains Moreover, the onset of SdH oscillations provides an estimate of the electronic mobility. Throughout the whole field range, f B2 (corresponding to a TSS) is the main peak of the FFT and the main contribution to the SdH. As a result, the surface mobility estimated from µ SS × B onset = 1 is µ SS 1200 cm 2 /V.s. For topological surface states, the corresponding mean free path is l e = v F µ SS m * /e = 28 nm, with v F = 5.4 × 10 5 m/s the Fermi velocity [43] and m * = √ πn SS2 /v F = 0.076 m e the cyclotron mass at this electronic density (m e the electron mass). This value of mean free path is of the same order of magnitude as reported in previous works [40,60]. As explained in the main text, anisotropic scattering of the topological surface states favors forward scattering, leading to a transport length l tr ≥ 8 l e [39,40]. Ballistic properties of the surface states may then be retained over more than 230 nm, larger than the 125 nm of the junction studied in the main text. We can futhermore estimate the diffusion coefficient of the bulk states assuming that the bulk mean free path is of the order of the surface state mean free path [40]. This leads to D = v bulk We present in this section the Shapiro maps obtained on another sample that displayed a very similar behavior in terms of IVs characteristics and Shapiro maps. The Bi 2 Se 3 flake is 20 nm thick and is contacted with Ti/V/Au (5/70/5 nm) electrodes, see Fig. S2a. The dimensions of the Josephson junction are W = 1.45 µm and L = 140 nm. The normal state resistance is R = 30 Ω and the critical current shown in Fig. S2b is 490 nA. This value is smaller than the one of sample LC106 discussed in the main text. We explain this by the presence of the Ti layer beneath the vanadium which reduces the superconducting proximity effect. The absence of hysteresis at 0.08 K results from this small value of critical current which limits Joule heating so that no thermal bistability develops [34]. Figure S3 displays two Shapiro maps of the junction measured at 1.6 GHz (a) and 4 GHz (b). Whereas the map at 4 GHz shows standard Shapiro steps, an even-odd effect develops in the map at lower frequency. The Shapiro step n = 1 develops after the step n = 2 on increasing P rf , and has a smaller amplitude than the other steps. Furthermore, the resistance at the first node of the critical current, indicated in Fig. S3a by the orange arrow, is less than at all other nodes. This is a reminiscence of the residual supercurrent observed in sample LC106 of the main text and consistent with the reduction of the first Shapiro step amplitude.

III. Frequency response of the RSJ model
We present in this section the frequency response of the RSJ model. The RSJ equation: features a characteristic phase relaxation time τ J = 2eRIc that acts as a low-pass frequency cutoff. For an rf frequency faster than 1/τ J , the phase cannot follow the rf drive thus precluding the phase locking necessary to generate Shapiro steps. To illustrate this low-pass frequency response, we computed three Shapiro maps at different frequencies higher and lower than 1/τ J . The resulting maps are displayed in Figure S4. For f rf τ J << 1 in Fig. S4a, standard Shapiro steps develop, signaling that the phase is well driven by the rf current, leading to phase-locking. For f rf τ J >> 1 in Fig. S4c, the Shapiro steps are nearly absent. We distinguish only the n = 1 step that starts to form at high P rf and develops a small current amplitude. This shows that the phase cannot follow the rf drive above the cutoff frequency and a much higher power is needed to induce phase-locking. Therefore, our simulations shows that the RSJ model behaves as a low-pass filter.

IV. Poisoning within the RSJ model
We implement the parity-change of the 4π-periodic channel by a stochastic sign-change of the current-phase relation: where n sw (t) is a positive integer that counts the number of switching events while solving the time-dependent RSJ equation (A.3) for a given dc current bias I and rf drive. The value n sw at a time t + dt where dt is the time-step of the differential equation solver is given by: where x ∈ [0 1] is a uniformly distributed random number, and I the rounding function to the nearest integer less than or equal to the value argument.  Figure S5a-c displays three computed Shapiro maps obtained for three different τ sw values. We observe in Fig. S5a that when τ sw is of the order of τ J , the poisoning suppresses the Shapiro steps and the IV 's are nearly ohmic (see the line-cut in Fig. S5d).
A Shapiro map computed in the opposite limit τ sw τ J is shown in Fig. S5c together with a line-cut in Fig. S5f. We see that a clear Shapiro steps develop, though with a small resistive slope due to the still non-zero poisoning. Consequently, when the mean switching time is much longer than the phase relaxation time, τ sw τ J , then poisoning becomes inefficient and lets the supercurrent and Shapiro steps nearly unaffected.
In Fig. S5g-i, we display histograms of the parity-lifetime for each τ sw /τ J ratio. Exponential decay fits in red lines confirm the Poisson-distribution of mean τ sw . Note that the apparent noise Fig. S5a-f results from the finite time-span over which the RSJ equation is solved.
These simulations demonstrate how poisoning can affect and ultimately suppress supercurrent and Shapiro steps within the RSJ model. The phase relaxation time τ J is therefore the key parameter to define the sensitivity to the poisoning dynamics.  Parameters of the RSJ model are R = 30 Ω, Ic = 500 nA, and f rf = 2 GHz. The switching time is τsw = 0.2, 0.5 and 2 ns for (a,d,g), (b,e,h) and (c,f,i) respectively. d-f, V in units of hf rf /2e versus I normalized to Ic extracted at the P rf indicated by the black arrows in a-c. For a large ratio τsw/τ 4π J in a and d, the Shapiro steps are barely visible. g-i, Histograms of the number of parity-switch versus parity lifetime for the respective Shapiro maps a-c. The red curves is a fit with an exponential decay which gives a mean value equal to τsw.

VI. Electron-phonon coupling in Bi2Se3 Josephson junctions
To solve the tRSJ model we estimated the electron-phonon coupling constant Σ Bi2Se3 of Bi 2 Se 3 using the formula for dirty metals and two-dimensional phonons [61] : with u l,t being the longitudinal and transverse sound velocities, v F the Fermi velocity l e the electronic mean free path, d the film thickness, Γ and ζ the gamma and zeta functions and β l a dimensionless parameter that we set to 1 [61]. With u l 2900 m/s and u t = 1700 m/s [62], we estimate the electron-phonon coupling constant to be Σ Bi2Se3 = 27.10 9 W.m −3 K −5 .
The validation of the above estimate can be done by fitting the T -dependence of the critical current hysteresis with the tRSJ model. Using the junction resistance, Σ Bi2Se3 and the measured switching current I sw (T ) for each T as input parameters, we adjusted the volume Ω of the junction entering the heat-balance equation to fit the measured retrapping current I r . Figure S7 displays I sw , I r and the resulting fit of I r for the junction studied in the main text.
With the above value of Σ Bi2Se3 , we obtain a volume equal to the product of the thickness d, width W and the sum of the length L and two halves of the electrode width (400 nm): 2.25 µm × 30 nm × (125 nm + 400 nm) = 3.5 × 10 −2 µm 3 . The excellent agreement between fit and the data together with the resulting volume that is consistent with the junction geometry validate our estimation of Σ Bi2Se3 . Overall, in our tRSJ model, the key parameter that controls the quasiparticle overheating is the product ΩΣ Bi2Se3 , which we deduce from fitting the retrapping current.

VII. Parameters of the RSJ model
Solving the RSJ model with the parameters of the junction presented in the main text (R = 7.5 Ω and I c = 7.3 µA) leads to a dense series of Shapiro steps, making the illustration of the two-channel RSJ model and the impact of thermal poisoning more delicate to capture visually.
Therefore, for the sake of clarity, we choose in Fig. 3 and 5 parameters that provide large Shapiro voltage steps and well developed even-odd effect, namely, I c = 500 nA, R = 30 Ω, and I 4π c /I c = 0.2, and f rf = 2 GHz. For Fig. 3d, f rf = 7 GHz.
Parameters for the thermal RSJ model in Fig. 5 are Σ Bi2Se3 = 27.10 9 W.m −3 K −5 , Ω = 3.5.10 −2 µm 3 (see section VI). For the activated poisoning, given that the exact Andreev spectrum is unknown, we conjecture that E M BS ∆ and thus parametrize ∆ − E M BS = 800 µeV [15]. We thus empirically choose τ 0 = 0.8 × 10 −19 s such that τ sw drops with the computed T qp of the Shapiro map below τ 4π J . Importantly, other choice of ∆ − E M BS with E M BS closer to the continuum would lead to virtually the same effects.