# Fileset [PhysRevResearch.6.033143.pdf](https://mdr.nims.go.jp/filesets/bbc3c720-ff25-46be-885f-8492c9d93960/download) ## Creator [Máté Kedves](https://orcid.org/0000-0002-2057-4891), Tamás Pápai, [Gergo ̋ Fülöp](https://orcid.org/0000-0003-3736-9488), [Kenji Watanabe](https://orcid.org/0000-0003-3701-8119), [Takashi Taniguchi](https://orcid.org/0000-0002-1467-3105), [Péter Makk](https://orcid.org/0000-0001-7637-4672), Szabolcs Csonka ## Rights [Creative Commons BY Attribution 4.0 International](https://creativecommons.org/licenses/by/4.0/) ## Other metadata [Self-heating effects and switching dynamics in graphene multiterminal Josephson junctions](https://mdr.nims.go.jp/datasets/e6fb092a-90fe-4f0b-ab48-3b74390c6582) ## Fulltext Self-heating effects and switching dynamics in graphene multiterminal Josephson junctionsPHYSICAL REVIEW RESEARCH 6, 033143 (2024)Editors’ SuggestionSelf-heating effects and switching dynamics in graphene multiterminal Josephson junctionsMáté Kedves ,1,2 Tamás Pápai,1,2 Gergő Fülöp ,1,3 Kenji Watanabe ,4 Takashi Taniguchi ,5Péter Makk ,1,2,* and Szabolcs Csonka1,3,61Department of Physics, Institute of Physics, Budapest University of Technology and Economics,Műegyetem rkp. 3, H-1111 Budapest, Hungary2MTA-BME Correlated van der Waals Structures Momentum Research Group, Műegyetem rkp. 3, H-1111 Budapest, Hungary3MTA-BME Superconducting Nanoelectronics Momentum Research Group, Műegyetem rkp. 3, H-1111 Budapest, Hungary4Research Center for Electronic and Optical Materials, National Institute for Materials Science, 1-1 Namiki, Tsukuba 305-0044, Japan5Research Center for Materials Nanoarchitectonics, National Institute for Materials Science, 1-1 Namiki, Tsukuba 305-0044, Japan6Institute of Technical Physics and Materials Science, Center for Energy Research, Konkoly-Thege Miklós út 29-33, 1121 Budapest, Hungary(Received 31 January 2024; accepted 24 April 2024; published 6 August 2024)We experimentally investigate the electronic transport properties of a three-terminal graphene Josephsonjunction. We find that self-heating effects strongly influence the behavior of this multiterminal Josephsonjunction (MTJJ) system. We show that existing simulation methods based on resistively and capacitively shuntedJosephson junction networks can be significantly improved by taking into account these heating effects. Wealso investigate the phase dynamics in our MTJJ by measuring its switching current distribution and findcorrelated switching events in different junctions. We show that the switching dynamics is governed by phasediffusion at low temperatures. Furthermore, we find that self-heating introduces additional damping that resultsin overdamped I−V characteristics when normal and supercurrents coexist in the device.DOI: 10.1103/PhysRevResearch.6.033143I. INTRODUCTIONMultiterminal Josephson junctions (MTJJs) consisting ofa single scattering region connected to multiple supercon-ducting terminals have attracted significant attention in recentyears. Theoretical works showed that MTJJs may enable mul-tiplet supercurrents [1–6], and the Andreev bound state (ABS)spectra of MTJJs can exhibit nontrivial topology and simu-late the band structure of Weyl semimetals [7–23]. Althoughsome of the theoretically proposed key features remain unob-served, recent experimental advances led to the observationof hybridized ABSs [24–27], broken spin degeneracy andground-state parity transitions [28], signatures of quartet su-percurrents [29–33], the Josephson diode effect [34–37], andtopological phase transitions [38], highlighting the versatilityof MTJJ devices.On the other hand, several experimental works found thatthe transport characteristics of MTJJs can be reasonably wellmodeled by a network of resistively and capacitively shuntedJosephson junctions (RCSJs) in which each pair of terminalsis connected by an RCSJ element. This relatively simpleapproach is able to qualitatively capture features of current-biased measurements, such as the coexistence of normal andsupercurrents between different terminals [31,32,39–41] and*Contact author: makk.peter@ttk.bme.huPublished by the American Physical Society under the terms of theCreative Commons Attribution 4.0 International license. Furtherdistribution of this work must maintain attribution to the author(s)and the published article’s title, journal citation, and DOI.multiplet resonances [31,41]. In spite of some agreementbetween simulations and measurements, these models ingeneral fail to quantitatively capture the observations whennormal and supercurrents coexist in the scattering region. Thislack of agreement can be attributed to heating effects dueto the presence of normal currents [39] that influence thesupercurrent flowing in other parts of the device. Furthermore,the observation of more exotic phenomena, such as multipletsupercurrents [31,41] and quantized transconductance [7,8],also rely on the presence of finite voltages between some ofthe terminals that necessarily imply the existence of normalcurrents and heating effects. Due to the large superconductinggap � of the terminals which prevents the outflow of hotelectrons, these heating effects can significantly modify thesuperconducting properties of MTJJs. Moreover, heating ef-fects can have an impact on the switching dynamics of singleJosephson junctions [42,43], which could be enhanced in thecase of MTJJs, due to the complex geometry and the nontrivialcurrent distribution.In this work we experimentally investigate a three-terminalgraphene Josephson junction and compare our current-biasedmeasurements to an RCSJ network model, which enables usto identify the limitations of these models. Next, we presentan improved simulation method, incorporating heating effectsdue to the presence of normal currents, which results in asignificantly better agreement with the measurements. Fur-thermore, we investigate the switching dynamics of our deviceand observe a nontrivial behavior of the switching current dis-tribution (SCD) at low temperatures that is governed by phasediffusion. We find that this behavior is also modified by theheating effects due to normal currents. Finally, we investigatethe charge-carrier-density dependence of the measured and2643-1564/2024/6(3)/033143(13) 033143-1 Published by the American Physical Societyhttps://orcid.org/0000-0002-2057-4891https://orcid.org/0000-0003-3736-9488https://orcid.org/0000-0003-3701-8119https://orcid.org/0000-0002-1467-3105https://orcid.org/0000-0001-7637-4672https://ror.org/02w42ss30https://ror.org/026v1ze26https://ror.org/026v1ze26https://ror.org/03ftngr23https://crossmark.crossref.org/dialog/?doi=10.1103/PhysRevResearch.6.033143&domain=pdf&date_stamp=2024-08-06https://doi.org/10.1103/PhysRevResearch.6.033143https://creativecommons.org/licenses/by/4.0/MÁTÉ KEDVES et al. PHYSICAL REVIEW RESEARCH 6, 033143 (2024)Riφi=i132IRI1I3I2IUIUIRV1V3V2hBNhBNS L GMoRe−2 −1 0 1 2I ( A)R μ−2−1012I(A)UμI*R0.00.51.0dV/dI(k)1UΩ−2 −1 0 1 2I ( A)R μ−2−1012I(A)Uμ0.00.51.0dV/dI(k)3UΩ−2 −1 0 1 2I ( A)R μ−2−1012I(A)Uμ0.00.51.0dV/dI(k)3UΩ−2 −1 0 1 2I ( A)R μ−2−1012I(A)Uμ0.00.51.0dV/dI(k)1UΩ(a) (b) (c)(d) (e) (f )123123FIG. 1. (a) Schematic representation of the multiterminal Josephson junction. Current biases IU and IR are applied via two separate contacts,and the third contact is grounded. Voltages Vi are measured between the three pairs of contacts. (b, c) Differential resistance dV1/dIU anddV3/dIU as a function of the current biases. In panel (b), white arrows illustrate the position of T -dependent measurement of I − V curves[Fig. 3(d)], and colored arrows correspond to bias values where SCD measurements were performed [Fig. 3(b)]. The white star symbol showsthe extended region where a finite voltage develops between all terminals simultaneously. White arrows in panel (c) point to resonant featuresattributed to MAR. (d) RSJ network model of our device. (e), (f) Simulated differential resistance maps analogous to panels (b) and (c),respectively. I∗R corresponds to the single-current-bias value of IR, where all three junctions switch to normal state simultaneously as IU isramped.simulated resistance maps, which gives us further insight intothe possible cooling mechanisms via which the dissipated heatescapes from the device.II. SUPERCURRENT CHARACTERIZATIONOF A MULTITERMINAL DEVICEA. Experimental resultsFigure 1(a) shows the schematic representation of our de-vice. A cross-shaped hBN/graphene/hBN heterostructure isconnected to three MoRe superconducting electrodes (op-tical image is shown in Appendix A). The separation ofneighboring contacts is around 150 nm. The charge-carrierdensity n in graphene can be tuned via the voltage appliedto the doped Si substrate that acts as a global back gate,while a 300-nm-thick SiO2 layer forms the gate dielectric.In our experiments, one of the electrodes is grounded andtwo independent dc current biases, IR and IU , are appliedvia the remaining two contacts, and differential voltages—V1,V2 and V3—between the three different pairs of terminalsare measured. Transport measurements were carried out in aLeiden dilution refrigerator at a base temperature of 40 mK(unless otherwise stated). Figures 1(b) and 1(c) show thedifferential resistance dV1/dIU (dV2/dIU )—obtained from themeasured V1 (V2) voltage by numerical differentiation withrespect to the current bias IU —as a function of IU and IR ata back-gate voltage of VBG = 10 V. Two main features can beidentified in such a differential resistance map, similarly toprevious experiments [31,32,34,35,39,40,44,45]. First, in thecenter, around small current bias values an extended super-conducting region of zero resistance can be observed. Second,superconducting arms [labeled by 1, 2, and 3 in Fig. 1(b)]are spreading out from this central superconducting region inmultiple directions. Comparing differential resistance mapsobtained from the measurements of V1 [Fig. 1(b)], V2 (seeAppendixes), and V3 [Fig. 1(c)], it is easy to realize thatthe central superconducting region is present in all cases,indicating that the whole sample is superconducting and su-percurrent can flow between all of the terminals. On the otherhand, each of the superconducting (SC) arms correspond tosupercurrent flowing between only two terminals, resulting inzero resistance in only one of the differential resistance maps,while a finite voltage develops between the remaining pairs ofterminals [e.g., the SC arm labeled by 1 shows zero resistancein Fig. 1(b) and a finite voltage develops in Fig. 1(c)]. Thisindicates that both normal and supercurrents can flow in thesample simultaneously.B. RSJ simulationPrevious works [32,39,40,45] showed that MTJJs can bedescribed to a large extent by a network of RCSJ elements.Here, we neglect capacitive effects and model our three-terminal JJ with three resistively shunted junctions (RSJs),one between any pair of contacts, as shown in Fig. 1(d). First,we present the results of this model and highlight its limita-tions in comparison with our measurements. Later, we showthat the agreement between measurement and simulation can033143-2SELF-HEATING EFFECTS AND SWITCHING DYNAMICS … PHYSICAL REVIEW RESEARCH 6, 033143 (2024)be improved by including self-heating effects in the model.As detailed in the Appendixes, the differential equations ofthis network model can be constructed from the Josephsonequations and Kirchhoff’s laws. The necessary input param-eters of the model are the resistances (Ri with i ∈ {1, 2, 3})and the critical currents (Ic,i) of the individual junctions. Rican be obtained from the measured differential resistances inthe normal state, at large bias currents. For these, we obtainR1 = 420 �, R2 = 1355 �, and R3 = 815 �. Furthermore, as-suming that our junctions are in the short junction limit andusing Ic,iRi ∝ �, it is possible to extract Ic,i from the mea-sured differential resistance maps as well. For these, we getIc,1 = 545 nA, Ic,2 = 170 nA, and Ic,3 = 280 nA, respectively.(See Appendix B for details on the extraction of parameters.)By numerically solving the set of differential equations forthe network of Josephson junctions and resistors, we obtaindifferential resistance maps as shown in Figs. 1(e) and 1(f).The model is capable of capturing the most prominent featuresof the measured differential resistance map: (i) the centralsuperconducting region and (ii) the superconducting arms,corresponding to the coexistence of normal current and su-percurrent. In the context of this model, the SC arms canbe further discussed. The total current between any pair ofterminals (I1, I2, and I3) is determined by the Kirchhoff andJosephson equations for a given IU and IR. It can be shown thatfor arbitrary IR, a single value of IU exists for each junctionfor which the total junction current Ii = 0 (see Appendix B).The ratio of IU /IR for which Ii = 0 is determined solely by thenormal resistances and is independent of IU and IR. Therefore,we expect to observe superconductivity in the vicinity of lineswith slopes defined by the normal resistances. We also notethat in this particular geometry, due to Kirchhoff’s law, whichstates that the sum of voltages in a closed loop has to be zero,a single junction cannot switch to the normal state alone; avoltage drop has to appear on either two or all three junctionssimultaneously. Therefore, outside the central SC region, theSC arms correspond to a configuration where only a singlejunction is superconducting and the other two resides in thenormal state. On the other hand, several missing features canalso be identified in the simulated resistance maps. The mostprominent example is the decay of superconductivity that canbe observed in the measurements along the superconductingarms. While the width of these arms in the simulated maps isconstant towards higher current bias values, in the measure-ments a clear narrowing of the zero-resistance regions can beobserved. Furthermore, in the measured resistance maps anextended region exists where all three junctions switch to thenormal state simultaneously [e.g., marked by star symbol inFig. 1(b)], whereas in the simulated maps, this simultaneousswitching of all three junctions can only be observed for asingle bias value I∗R [marked also by vertical dashed line inFig. 1(e)]. Finally, multiple resonant features [e.g., marked bywhite arrows in Fig. 1(c)] are visible in the measurementsparallel to the superconducting arms that are attributed tomultiple Andreev reflections (MAR) [44] that cannot be ac-counted for by our simple model. We also note that MARfeatures are more pronounced for junction 1 than the othertwo junctions. This implies a larger contact transparency andmay explain the larger critical current and lower resistanceextracted for this junction.(a) (b)(c) (d)0.2 0.4P (nW)J150175200225I(nA)c1T (mK)107080060039024090401.5 1.75 2.0T (K)e150175200225I(nA)c1T (mK)10708006003902409040−2 −1 0 1 2I ( A)R μ−2−1012I(A)Uμ1 2 3T (K)e−2 −1 0 1 2I ( A)R μ−2−1012I(A)UμI*R0 0.5 1dV /dI (k )1 U ΩFIG. 2. Critical current of junction 1 Ic,1 along the correspondingsuperconducting arm as a function of (a) heating power PJ and(b) electronic temperature Te calculated assuming only phonon cool-ing. (c) Simulated differential resistance map taking into account theelevated electronic temperature due to normal current flowing in thedevice. (d) Simulated map of Te as a function of current biases.III. SELF-HEATING EFFECTSThe narrowing of the superconducting arms is attributedto Joule heating from the dissipative normal currents in thescattering region [39]. Due to the large superconducting gapof the MoRe that prevents hot electron diffusion towards theleads, the electron system can only dissipate heat via electron-phonon coupling. In this case the dissipated power towardsthe substrate is given by Pe−ph = �(T δe − T δ ) [46], where� is the electron-phonon coupling constant, and Te and Tare the electron and phonon bath temperatures, respectively.Following along the lines of Ref. [39], we determine � fromthe temperature dependence of Ic,1 along the correspondingSC arm. For this we measure the switching current Is,1 forjunction 1 by sweeping IU at different values of IR and bathtemperatures. Is,1 is then defined as the value of IU whereV1 crosses a certain threshold voltage (20 µV) correspondingto the switching from the SC to the normal state. As men-tioned earlier, in this current-biasing scheme, IU and IR donot directly correspond to the junction currents I1, I2, or I3.However, since along the SC arm supercurrent only flows injunction 1, it is possible to calculate the junction’s criticalcurrent Ic,1 from Is,1 (see Appendix B). Moreover, as it isdetailed later, the switching current of a Josephson junction isprone to fluctuations due to thermal effects. To eliminate thesefluctuations, we take the average of 10 000 measurementsto determine the average switching current Is,1. Next, wecalculate the power PJ dissipated in the normal regions fromJoule heating as PJ = IUV2. Figure 2(a) shows the measuredcritical current Ic,1 as a function of PJ for different T bath033143-3MÁTÉ KEDVES et al. PHYSICAL REVIEW RESEARCH 6, 033143 (2024)temperatures. As it can be seen from the figure, the increasedheating power leads to the decrease of the switching current,similarly to the increased bath temperature. Our assumptionis that Te is homogeneous in the device and the critical cur-rent value is defined by Te independently from whether itoriginates from bath heating or current dissipation. We thendetermine the value of � for which Ic,1 as a function of thecalculated equilibrium electron temperature Te scales onto asingle curve. As it is discussed later, by assuming δ = 4, weobtain � = 25 pW/K4. This is shown in Fig. 2(b), where allthe curves fall on top of each other. Although it is challengingto determine the exact active area of our device, we estimatethat � scaled by the graphene’s area yields ∼100 W/m2 K4.This is an order of magnitude larger than the value obtainedin Ref. [39] (∼10 W/m2 K3) and is significantly larger thanthe value obtained for large-area, nonencapsulated graphenedevices [46] (<50 mW/m2 K4). The authors of Ref. [39] alsospeculate that electron-phonon coupling can be enhanced bythe presence of the hBN substrate and by scattering at theedges of the graphene layer. Since our device area is aboutan order of magnitude smaller than the device studied inRef. [39], scattering at the edges could be even more signifi-cant and could explain the larger value obtained for � scaledby the graphene’s area.To take the effects of self-heating into account in our sim-ulations, we perform a fixed-point iteration based on the RSJmodel introduced previously. First, we perform the previoussimulation with the experimentally determined Ri and Ic,i pa-rameters for all IU and IR. We then calculate the Joule heatdissipated in the whole network as PJ = ∑i V 2i /Ri. From PJwe can obtain the equilibrium electron temperature Te usingthe electron-phonon coupling model for all IU and IR biascurrents. Finally, we take into account the elevated temper-ature using an Ic(Te) function, which we reconstruct fromthe temperature-dependent measurements shown previouslyin Figs. 2(a) and 2(b) and from the temperature dependenceof the central superconducting region (see Appendix B formore details). We then iterate this process to achieve a self-consistent solution using the modified Ic,i values in our RSJmodel, which now also depend on the applied IU and IRcurrent biases.Figure 2(c) shows the simulated dV1/dIU map obtained inour model with self-heating. Compared to Fig. 1(d), severalimprovements can be observed. First of all, the narrowingof the SC arms is qualitatively reproduced. The remainingquantitative difference could be explained by the incorrectreconstruction of the Ic(Te) function. Secondly, the improvedsimulation method is capable of producing an extended edgeon the contour of the central SC region where all threejunctions switch to the normal state simultaneously. It isalso worth noting that the simulated resistance map is in-version symmetric, in contrast to the measurements wherethe sweep direction of the bias currents results in a slightlyasymmetric central SC region. Finally, Fig. 2(d) shows themap of Te, illustrating that the heating outside the centralSC region is significant, increasing the equilibrium temper-ature to a few Kelvins, an order of magnitude above thebath temperature, in agreement with our measurements shownin Fig. 2(b).(a) (b)(d)(c)0.4 0.5 0.6 0.7I ( A)U μ0.000.050.100.15V 1(mV)I = 0.1 AR μ-0.2 Aμ0 0.5 1T (K)46σ(nA)T (mK)107060040−0.5 0.0 0.5I ( A)U μ-0.500.5I(A)iμIc, 1Ic, 2Ic, 3junction 1junction 2junction 30.20 0.35 0.50I ( A)s, 1 μ036counts(a.u.)J10.00.10.20.30.40.5I ( A)R μ0.38 0.44 0.50I ( A)s, 1 μ0.00.51.0counts(a.u.)J1T (mK)107080060039024040FIG. 3. (a) Current distribution between the three junctions in thecentral superconducting region obtained from the RSJ model (dashedlines) and from numerically minimizing the Josephson energy of thewhole network (symbols) for IR = 0.1 µA. Dotted lines show thecritical current of each junction. (b) SCD for junction 1 measuredat different IR in the central superconducting region obtained from10 000 measurements. (c) Temperature dependence of the SCD atIR = 0.1 µA. The narrowing of the SCD with increasing temperatureis consistent with phase diffusion. (d) Averaged I−V curves obtainedfrom 10 000 individual measurements in the central SC region (IR =0.1 µA) and in the SC arm of junction 1 (IR = −0.2 µA) for differenttemperatures. Inset: Standard deviation σ of the SCD as a functionof T for IR = 0.1 µA.IV. SWITCHING DYNAMICSIn the following we further investigate the interplay be-tween the three junctions in the regions where all junctionsswitch to the normal state simultaneously around I∗R [alsoshown in Fig. 2(c)]. Figure 3(a) shows the current in eachjunction as a function of IU for IR = 0.1 µA � I∗R [orangearrow in Fig. 2(c)], obtained from our simulation with self-heating (dashed lines). Although Ii cannot be obtained fromour measurements, we can also calculate the current in eachjunction as long as all junctions are superconducting by nu-merically minimizing the Josephson energy (symbols). As isvisible in Fig. 3(a), this method is consistent with our sim-ulation. The dotted horizontal lines show the critical currentof the respective junctions. As is visible also in Fig. 1(e),junction 1 is far below its critical current when junction 2 and3 reach their respective critical currents. Therefore, withouttaking self-heating into account, we only expect junction 2and 3 to switch together. However, when heating is included[Fig. 2(c) and measurements on Figs. 1(b) and 1(c)], all threejunctions switch at the same IU . From these we can infer thatjunctions 2 and 3 switch together and junction 1 switches im-mediately afterwards due to heating from the other junctions.On the other hand, for IR ∼ I∗R , all three junctions reach their033143-4SELF-HEATING EFFECTS AND SWITCHING DYNAMICS … PHYSICAL REVIEW RESEARCH 6, 033143 (2024)critical currents simultaneously and heating should play norole in the switching process, while for IR � I∗R [gray arrowin Fig. 2(c)], junctions 1 and 2 switch together, and junction3 switches due to heating. Based on the previous arguments,we emphasize that the observation of this correlated switchingof all three junctions in an extended region along the borderof the central SC region is strong evidence for self-heatingeffects.To gain insight into the dynamics of these correlatedswitchings, it is essential to investigate not only the averageswitching current but also its distribution. Figure 3(b) showsthe switching current distribution (SCD) obtained from themeasurement of V1 for different IR at 40-mK base temperature.The investigated values of IR are also indicated at the top ofFig. 1(b) by colored arrows. Each distribution is obtained bysweeping IU and detecting the switching current using thepreviously defined threshold voltage. This process is repeated10 000 times, and a distribution of switching current values isobtained. Interestingly, we observe that the width of the SCDis greatly tunable by IR [Fig. 3(b)]. We find that the standarddeviation σ of the SCD, which describes the width of thedistribution, increases by a factor of 2. This broadening of theSCD could be explained by the different junctions that switchsimultaneously at different IR. As junction 1 takes over the roleof junction 3 with increasing IR, the sum of the critical currentsof the two junctions that switch simultaneously increases,which could lead to a wider distribution. It is also important tonote that during the measurement of the SCD of junction 1, wesimultaneously recorded the SCD of junction 2, obtained fromthe appearance of a finite V2, and find that the distributionsare identical and the switching events of the two junctions areindistinguishable within the timescales of our measurement(see Appendix E). This suggests that the thermalization of thedevice is faster than our data acquisition.To further investigate the escape dynamics of our device,we measure the temperature dependence of the SCD alongthe contour of the central SC region. This is shown in Fig. 3(c)for IR = 0.1 µA, and a similar trend is observed for all inves-tigated values of IR inside the central SC region. It is clearlyvisible that the SCD gets narrower with increasing T , whichis in stark contrast to the thermally activated behavior, as theSCD is expected to broaden with temperature. This is furtherconfirmed by calculating σ as a function of bath temperaturefor IR = 0.1 µA, which is shown in the inset of Fig. 3(d).Here, an ∼40% decrease of σ is visible in the investigatedtemperature range. The observed narrowing of the SCD withincreasing T is a consequence of phase diffusion due tothermally activated escape and retrapping, and is consistentwith previous observations in moderately damped Josephsonjunctions [47] and planar Josephson junctions [48]. However,it is important to note that we observe the narrowing of theSCD in the whole available temperature range and do not findthe broadening of the SCD due to thermally activated escape,even for the lowest temperatures. This suggests that phasediffusion is significant even at base temperature.We performed similar measurements along the SC armof junction 1. Here, we find a different behavior and wecannot resolve a clear SCD. Figure 3(d) shows the aver-aged I−V curves of the 10 000 individual measurements forIR = −0.2 µA and IR = 0.1 µA [white arrows in Fig. 1(b)] fordifferent temperatures. In the central SC region, for IR =0.1 µA, a sharp transition between the SC and normal statescan be seen. In this case the curvature of the averaged I−Vcurves results from averaging curves with fluctuating switch-ing currents. On the other hand, for IR = −0.2 µA, along theSC arm of junction 1, a smooth transition is observed, indicat-ing that a finite voltage develops below the switching current.This is also consistent with the theoretical expectations formoderately damped Josephson junctions at higher tempera-tures [47]. As T increases, the thermally activated retrappingresults in a significant damping and the junctions becomeoverdamped. This is further confirmed by the T dependence ofthe curves. For IR = 0.1 µA, as T is increased, the switchingcurrent decreases [also visible in Fig. 3(c)]. However, forIR = −0.2 µA, along the SC arm of junction 1, the effect ofincreasing T is negligible; the increase of T rather makesthe transition between the SC and normal states smoother, asexpected for overdamped junctions. It is also consistent withthe self-heating picture, since increasing the bath temperaturehas less effect on the electronic temperature when a largeheating power is already present due to the normal currents inthe device. Therefore, we conclude that the switching of ourmultiterminal device is determined by phase diffusion at lowertemperatures along the contour of the central SC region andshow overdamped characteristics along the SC arm of junction1 due to the increased temperature.V. CHARGE CARRIER DENSITY DEPENDENCEFinally, we investigate the dependence of the differentialresistance maps on the applied back-gate voltage VBG. Asmentioned earlier, the exponent of the electron-phonon cool-ing power formula δ can be 3 or 4, depending on electronicmean free path lmfp and temperature [46]. At the relatively lowtemperatures accessed in our measurements, δ = 4 describesphonon cooling in clean devices where lmfp is large, whileδ = 3 corresponds to phonon cooling modified by impurityscattering in devices with small lmfp. Furthermore, the expres-sion for � is also different in the two limits. In the clean limit� = π2D2|EF |k4B/15ρM h̄5v3F s3, where D is the deformationpotential of graphene, which describes the electron-phononcoupling strength, ρM is the mass density of graphene, vF =106 m/s is the Fermi velocity, EF = h̄vF√πn is the Fermi en-ergy, and s = 2 × 104 m/s is the speed of sound in graphene.It can easily be shown that in the clean limit � ∝ √n, whilein the dirty limit the expression is modified and � becomesindependent of n [46]. Figure 4 shows the measured and simu-lated resistance maps for different VBG and δ = 4. We scale �according to the√n dependence, and n is calculated accord-ing to a planar capacitor model based on the hBN and SiO2dielectrics (see Appendix C). From Fig. 4 it is visible thatthe qualitative trend is reproduced well. Here, we assumedδ = 4, but note that a reasonably good agreement can alsobe achieved by taking δ = 3 and a constant � = 30 pW/K3(see Appendix C). Although the overall qualitative agreementbetween measurement and simulation is good, some differ-ences can still be observed. Most notably, some of the SCarms persist up to larger current bias values in the measure-ments, especially noticeable for VBG = 2 V. This could beexplained by the appearance of additional cooling paths. As033143-5MÁTÉ KEDVES et al. PHYSICAL REVIEW RESEARCH 6, 033143 (2024)−2−1012I(A)UμV = 10 VBG Simulation−2−1012I(A)Uμ6 V−2 −1 0 1 2I ( A)R μ−2−1012I(A)Uμ2 V−2 −1 0 1 2I ( A)R μ0 1 2dV /dI (k )1 U ΩFIG. 4. Back-gate dependence of the measured (left) and sim-ulated (right) differential resistance maps. Simulations were per-formed with δ = 4 and � = 25 pW/K4 for VBG = 10 V, and � wasscaled according to the√n dependence expected for the clean limitof electron-phonon coupling.Te is increased up to a few Kelvins, kBTe becomes comparableto �, allowing quasiparticles to diffuse into the MoRe leads.Furthermore, we assumed that Te is homogeneous in the wholedevice, which does not necessarily hold for large heatingpowers. The inhomogeneity of Te could significantly modifythe ratio of normal and SC segments of the scattering region,and as a result, the estimated input parameters of our modelwould become increasingly inaccurate with increasing heatingpowers.VI. CONCLUSIONSIn conclusion, we have measured three-terminal grapheneJosephson junctions and investigated the heating effects andjunction dynamics in this multiterminal system. We haveshown that a significant improvement can be achieved overexisting RCSJ models for MTJJs by incorporating heatingeffects into the simulation method. By considering only Jouleheating from the normal currents in the device and electron-phonon coupling as a cooling mechanism, we were ableto obtain the narrowing of the SC arms that is commonlyobserved in experiments and the simultaneous switching ofall junctions. By measuring the charge-carrier-density depen-dence of the differential resistance maps, we could infer thelimitations of our model and suggest that, for significantlyincreased electronic temperatures, new cooling mechanismsmight become available. We propose that by including ad-ditional cooling terms, such as the outflow of hot electronsvia the SC terminals, our model could be further improved.ABCDV3V1V2IRIU(b)(a)−10 0 10V (V)BG−101I(A)Uμ0 1 2dV /dI (k )1 U ΩFIG. 5. (a) Optical microscopic image of the device with theschematic illustration of the measurement geometry. Scale bar is2 µm. (b) Differential resistance of junction 1 dV1/dIU as a functionof current bias IU and back-gate voltage VBG for IR = 0.Furthermore, from the investigation of the SCD, we con-cluded that the switching from the central SC region to thenormal state is governed by phase diffusion, even at verylow temperatures. As the temperature is increased due toself-heating, this phase diffusion modifies the characteristicsof the device, resulting in smooth I−V curves resemblingoverdamped Josephson junctions. Building on these results,future experiments could focus on the phase-biasing of MTJJsand inductance measurements using rf techniques in the SCstate, where self-heating effects are absent.Raw measurement data and simulation results are availableat [49].ACKNOWLEDGMENTSThis work acknowledges support from the Topograph,MultiSpin, and 2DSOTECH FlagERA networks; OTKAGrants K138433 and K134437; VEKOP Grant 2.3.3-15-2017-00015; and an EIC Pathfinder Challenge grant QuKiT. Thisresearch was supported by the Ministry of Culture and Innova-tion and the National Research, Development and InnovationOffice within the Quantum Information National Laboratoryof Hungary (Grant No. 2022-2.1.1-NL-2022-00004), and bythe FET Open AndQC and SuperGate networks and by theEuropean Research Council ERC project Twistrain. We ac-knowledge COST Action CA 21144 superQUMAP. K.W. andT.T. acknowledge support from the JSPS KAKENHI (GrantsNo. 20H00354 and No. 23H02052) and the World PremierInternational Research Center Initiative (WPI), MEXT, Japan.M.K. and T.P. fabricated the device. Measurements wereperformed by M.K. and T.P. with the help of P.M. and Sz.Cs. M.K. and T.P. did the data analysis. Simulations wereperformed by G.F., M.K., and T.P., and K.W. and T.T. grew thehBN crystals. All authors discussed the results and worked onthe manuscript; M.K. and P.M. wrote the paper. The projectwas guided by P.M. and Sz.Cs.APPENDIX A: DEVICE GEOMETRYAND MEASUREMENT SETUPThe measured sample is shown in Fig. 5(a). The dry-transfer technique with polycarbonate/polydimethylsiloxanestamps was employed to stack hBN (20 nm, top)/single layer033143-6SELF-HEATING EFFECTS AND SWITCHING DYNAMICS … PHYSICAL REVIEW RESEARCH 6, 033143 (2024)−2−1012I(A)Uμ1 2 3−2 −1 0 1 2I ( A)R μ−2−1012I(A)Uμ1−2 −1 0 1 2I ( A)R μ2−2 −1 0 1 2I ( A)R μ3−1 0 1V (mV)i0 1 2dV /dI (k )i U Ω(a) (b) (c)(d) (e) (f )FIG. 6. (a)–(c) Raw measured voltages V1, V2, and V3 as a function of IU and IR. (d)–(f) Differential resistances calculated frompanels (a)–(c).graphene/hBN (35 nm, bottom). To fabricate electrical con-tacts, we used electron-beam lithography patterning followedby a reactive-ion-etching step using a CHF3/O2 mixture andfinally deposited MoRe (50 nm) by dc sputtering. As it is vis-ible on the optical microscopic image in Fig. 5(a), four MoRecontacts were fabricated; however, one of the contacts failedto contact the graphene layer, resulting in a three-terminaldevice as presented in the main text. The separation of neigh-boring contacts is around 150 nm. The heterostructure aroundthe cross-shaped region was etched away using reactive-ionetching with SF6/O2 mixture.Transport measurements were carried out in a Leiden dilu-tion refrigerator at a base temperature of 40 mK (unless other-wise stated). Measurements were performed using a NI USB6341 measurement card. In each measurement, contact A wasgrounded and the dc current biases IR and IU were applied via1-M� preresistors to contacts B and C, respectively. Duringthe SCD measurements for fixed values of IR, IU was rampedfrom 0 to 1 µA at a ramp rate of 100 µA/s while the voltagesbetween two different pairs of terminals were simultaneouslymeasured. Figure 6 shows the measured raw voltages as afunction of the current biases, corresponding to the differentialresistance maps shown in Fig. 1 of the main text.Figure 5(b) shows the differential resistance of junction 1as a function of VBG and IU , showing a highly tunable criticalcurrent with VBG, as is common for graphene devices. Thecritical current can be tuned to zero near the charge neu-trality point, and we observe a significantly smaller criticalcurrent for negative VBG that we attribute to doping from theMoRe contacts and formation of a p-n junction at the MoReinterface [50].APPENDIX B: RSJ SIMULATIONSAs discussed in the main text, we start our simulation bysolving an RSJ network model. For our three-terminal device,this consists of three blocks of resistively shunted Josephsonjunctions as shown in Fig. 1(d) of the main text. The ith blockis described by a resistor with resistance Ri and the phasedifference of the Josephson junction ϕi. The normal currentflowing in the resistor is given by IN,i = Vi/Ri, where Vi isthe voltage drop on the RSJ block. We employ a sinusoidalcurrent-phase relation, and the supercurrent flowing in theJosephson junction is given by Is,i = Ic,i sin ϕi. According tothe corresponding Josephson equation, the time derivative ofthe phase difference is given by ϕ̇i = 2eVi/h̄. With these, onecan obtain the differential equation of a single RSJ block:Ii = Ic,i sin ϕi + h̄2eRiVi,where Ii is the total current flowing in the ith block. In-troducing the external current biases IU and IR and thesuperconducting phases of the corresponding leads ϕU and ϕRaccording to Fig. 7(a), choosing the phase of the groundedterminal as zero, and applying Kirchhoff’s law, one can endup with a set of coupled differential equations for the completeRSJ network:da1dt= 2eh̄[IU − Ic,2 sin(−ϕU ) − Ic,3 sin(ϕR − ϕU )],da2dt= 2eh̄[IR − Ic,1 sin(−ϕR) + Ic,3 sin(ϕR − ϕU )], (A1)033143-7MÁTÉ KEDVES et al. PHYSICAL REVIEW RESEARCH 6, 033143 (2024)−2 −1 0 1 2I ( A)R μ−2−1012I(A)Uμ0.0 0.5 1.0dV /dI (k )1 U Ω(b)(a)−1 0 1I ( A)U μ0123dV/dI(k)iUΩI = 0 AR μjunction 1junction 2junction 3FIG. 7. (a) Differential resistance maps dV1/dIU . Dashed lines il-lustrate the obtained slopes of the SC arms. (b) Measured differentialresistances for IR = 0 (markers). Solid lines show the simulated dif-ferential resistances with our improved method, taking self-heatingeffects into account.wherea1 = ϕR − ϕUR3− ϕUR2,a2 = −ϕR − ϕUR3− ϕUR1,and we made use of the fact that ϕ3 = ϕR − ϕU . By numeri-cally solving equation system (A1), we obtain the stationaryϕi phase differences and Vi voltages from which both the nor-mal In,i and supercurrents Is,i in each block can be calculatedfor a given IU and IR.1. Determination of junction parametersAs mentioned in the main text, to quantitatively match thesimulations to our measurement, we determine Ri and Ic,i fromthe measured differential resistance maps. First of all, it is easyto show that the ratio IU /IR for which Is,i = 0, correspondingto the slope of the SC arms, is determined by the normalresistances asα = −R2 + R3R2,β = − R1R1 + R2,γ = R1R2, (A2)for junctions 1, 2, and 3, respectively. For these we obtainα = −1.6, β = −0.34, and γ = 0.31 from the measured dif-ferential resistance maps at VBG = 10 V. These are shownwith dashed lines in Fig. 7(b). Since these equations are notindependent, we also calculate the differential resistances inthe normal state where only normal currents are flowing asRI = dV1dIU= R1R2R1 + R2 + R3,RII = dV2dIU= R2(R1 + R3)R1 + R2 + R3,RIII = dV2dIU= R2R3R1 + R2 + R3. (A3)Combining equation systems (A2) and (A3), one can showthat R2 = RI (γ − α)/γ and the normal resistances can becalculated. For these, we obtain R1 = 420 �, R2 = 1355 �,and R3 = 815 �, respectively. Having obtained the normalresistances, it is also possible to calculate the junction crit-ical currents Ic,i. First, we calculate the superconductingcoherence length in graphene. Since the length of our junc-tions is smaller than 200 nm, well below the typical meanfree path for similar graphene devices, we assume ballis-tic conduction. Using � = 1.2 meV for the SC gap of theMoRe contacts [51,52], the coherence length is given by ξ =h̄vF /π� ≈ 200 nm. Therefore, we conclude that our junc-tions are in the short, ballistic limit, which implies that Ic,iRi ∝�. This allows us to calculate Ic,i from the measured differ-ential resistance maps using the previously calculated normalresistances. Using this, it can be shown that for IR = 0, thetotal critical current is given by Ic,tot = Ic,2 + Ic,3 = Ic,2(1 +R2/R3) and the individual junction critical currents Ic,i canbe calculated using the Ri normal resistances. We associateIc,i with the values obtained from the differential resistancemaps measured at base temperature. For these, we obtainIc,1 = 545 nA, Ic,2 = 170 nA, and Ic,3 = 280 nA.2. Determination of Ic,i(Te)As discussed in the main text, to include heating effectsin our simulations we perform a fixed-point iteration. Thepseudocode for this algorithm is shown in Algorithm 1. First,we solve the RSJ network model with the experimentallyobtained parameters and calculate the Joule heating power asPJ = ∑i V 2i /Ri and the equilibrium electron temperature asTe = 4√T 4 + PJ/�, using � = 25 pW/K4 as obtained fromALGORITHM 1. Iterative procedure for the self-consistent calculation of junction currents and electronic temperature.function CALCULATE_MTJJ niter = 10, δ, �, Ri, Ic,i(Te) for i ∈ {1, 2, 3}, TbathTe ← Tbath for all IU , IR � Initializationfor niter repetitions dofor all IU , IR in the rangeVi, Ii ← solve ODE set using Ic,i = Ic,i(Te) � Equation system (A1)P ← V 21 /R1 + V 22 /R2 + V 23 /R3T newe ← solve P = �(T δe − T δbath ) for Te � Assumes Pe−ph = PJend forend forreturn Ii(IU , IR),Vi(IU , IR ), Te(IU , IR )end function033143-8SELF-HEATING EFFECTS AND SWITCHING DYNAMICS … PHYSICAL REVIEW RESEARCH 6, 033143 (2024)0 1 2 3T (K)e0200400600I c,1(nA)T (mK)10708006003902409040Σ = 25 pW/K , = 44 δFIG. 8. The experimentally obtained Ic,1(Te) function.Triangles show Ic,1 obtained from the central SC region, andcircles correspond to the values extracted from the SC arm ofjunction 1.the temperature-dependent measurements (see Fig. 2. of themain text) and assuming homogeneous temperature distri-bution in the device. The next step is to take the effect ofthe elevated electron temperature into account via the Ic(Te)dependence. We construct this function from our temperature-dependent measurements. For this we have to consider twodifferent regimes. First, in the central SC region, as discussedpreviously, the individual junction critical currents can be cal-culated using the normal resistances. Assuming that the ratioof the resistances does not change with temperature, we canobtain Ic,1 by taking Ic,tot as the mean of the SCDs measured atIR = 0 for different T [Fig. 3(c) of the main text]. Moreover,since in this region all junctions are superconducting, we cantake Te = T , as there is no Joule heating.−2−1012I(A)Uμ#1 #3 #10−2 −1 0 1 2I ( A)R μ−2−1012I(A)Uμ−2 −1 0 1 2I ( A)R μ−2 −1 0 1 2I ( A)R μ0 1 2dV /dI (k )1 U Ω0 2 4ΔT (K)e(a) (b) (c)(d) (e) (f )FIG. 9. (a)–(c) Simulated differential resistance maps after 1, 3, and 10 iterations, respectively. (d)–(f) Change of electronic temperature�Te = Tn − Tn−1, where n is the iteration step.TABLE I. Charge-carrier densities n and � in the case of δ = 4,corresponding to the values of VBG for which the differential resis-tance maps were measured and simulated.VBG (V) n (1012 cm−2) �δ=4 (pW/K4)10 0.74 256 0.48 202 0.22 14Next, we consider the SC arm of junction 1. Utilizing theprevious definition of the slope α of the SC arm of junction1, for a given IR the supercurrent in junction 1 is zero forIU = αIR. Furthermore, since along the SC arm only junction1 is superconducting and the remaining two junctions arein the normal state, we can calculate the ratio of IU that isflowing towards junction 1. Combining these, the net currentof junction 1 is given by I1 = (IU − αIR)R2/(R2 + R3). In thiscase we define the average switching current of junction 1Is,1 as the value of IU for which V 1 exceeds the predefinedthreshold voltage (20 µV), where V 1 is the average voltageobtained from averaging 10 000 individual measurements.From this we calculate the critical current of junction 1 asIc,1 = (Is,1 − αIR)R2/(R2 + R3). The obtained values of Ic,1for different Te are shown in Fig. 8. To find the value of Ic,1for any Te, we linearly interpolate and extrapolate. Finally, toget Ic,2 and Ic,3, we simply scale the Ic,1(Te) function accord-ing to the ratio of normal resistances, based on our previousarguments.The simulated differential resistances for IR = 0 are shownwith solid lines in Fig. 7(c). As is visible, the simulated curves033143-9MÁTÉ KEDVES et al. PHYSICAL REVIEW RESEARCH 6, 033143 (2024)−2−1012I(A)UμV = 10 VBG 6 V 2 V−2 −1 0 1 2I ( A)R μ−2−1012I(A)UμSimulationδ=3−2 −1 0 1 2I ( A)R μ−2 −1 0 1 2I ( A)R μ0 1 2dV /dI (k )1 U ΩFIG. 10. Simulation with δ = 3 and constant � = 30 pW/K3.qualitatively match the measured points for IU > 0. For neg-ative IU , the retrapping to the SC state happens later in themeasurements than in the simulations. We attribute this also tothe elevated temperature due to self-heating, as the simulatedcurves do not take into account the sweep direction of thecurrent bias.3. Iteration processTo further illustrate the fixed-point iteration method, weshow the simulated differential resistance map dV1/dIU af-ter different numbers of iteration in Fig. 9. The first step[Fig. 9(a)] corresponds to the simulation without takingheating into account, also shown in Figs. 1(e) and 1(f) ofthe main text. After three iterations [Fig. 9(b)], the mainfeatures of the measured resistance maps are well repro-duced. Figure 9(c) shows the final result after ten rounds ofiteration, which only shows minor differences compared toFig. 9(b). Figures 9(d)–9(f) show the change of electronictemperature �Te = Tn − Tn−1, where n is the iteration stepand T0 = 40 mK is the base temperature. It can be seen thatwhile the electronic temperature is drastically modified forthe first step, later iterations only result in minor changes,indicating the convergence of our simulations.APPENDIX C: ADDITIONAL SIMULATIONSAs mentioned in the main text, we can also perform thescaling of Ic,1 along the SC arm of junction 1 using δ = 3,corresponding to the dirty limit of electron-phonon coupling.This scaling yields � = 30 pW/K3. We also construct theIc(Te) function using this modified � and simulate the dif-ferential resistance maps analogous to Fig. 4 of the main text.In this case the expression for � is modified; it is given by� = 2ζ (3)D2|EF |k3Bπ2ρM h̄4v3F s2lmfp. It can be shown that, in this case, � isindependent of n. The simulated resistance maps for δ = 3and � = 30 pW/K3 are shown in Fig. 10.(a) (b) (c)-0.5 -0.25 0 0.25 0.5V (mV)2-0.5-0.2500.250.5V(mV)1-0.5 -0.25 0 0.25 0.5V (mV)2-0.5 -0.25 0 0.25 0.5V (mV)30.0 0.3dV /dI (k )1 U Ω0.0 0.9dV /dI (k )2 U Ω0.0 0.7dV /dI (k )3 U ΩFIG. 11. Measured differential resistances dVi/dIU as a function of the measured voltages.033143-10SELF-HEATING EFFECTS AND SWITCHING DYNAMICS … PHYSICAL REVIEW RESEARCH 6, 033143 (2024)(a) (b)(c) (d)0.34 0.42 0.50I ( A)s, 1 μ0.00.51.0counts(a.u.)J1I = 0.2 AR μT (mK)1070800600390240400.2 0.3 0.4I ( A)s, 1 μ0.00.51.0counts(a.u.)J1I = 0.4 AR μT (mK)1070800600390240400.34 0.42 0.50I ( A)s, 2 μ0.00.51.0counts(a.u.)J2I = 0.2 AR μT (mK)1070800600390240400.2 0.3 0.4I ( A)s, 2 μ0.00.51.0counts(a.u.)J2I = 0.4 AR μT (mK)107080060039024040FIG. 12. (a), (b) Additional SCD data for junction 1 measuredat IR = 0.2 µA and IR = 0.4 µA, respectively, for different tempera-tures. (c), (d) SCDs for junction 2, simultaneously measured with theSCDs for junction 1.As detailed in the main text, for δ = 4, � is scaled ac-cording to an√n dependence. The � values for each VBGcan be found in Table I. We also present the charge-carrierdensities n for the different VBG values where the differentialresistance maps were measured and simulated in Table I.We determine the back-gate voltage of the charge neutralitypoint VCNP = −1.4 V from the gate-dependent measurementshown in Fig. 5(c). Using this, the carrier density is givenby n = αBG(VBG − VCNP). The lever arm of the back gate iscalculated according to a planar capacitor model as αBG =ε0/e(dSiO2/εSiO2 + dhBN/εhBN)−1, where ε0 is the vacuum per-mittivity, e is the elementary charge, εSiO2 = 4 (εhBN = 3.3),and dSiO2 = 300 nm (dhBN = 35 nm) are the dielectric con-stant and thickness of SiO2 (hBN), respectively.APPENDIX D: MULTIPLE ANDREEV REFLECTIONSFigure 11 shows the differential resistances dVi/dIU plot-ted as a function of the measured voltages Vi. We observeresonant features that are attributed to multiple Andreev−20 0 20B (mT)−101I(A)Rμ012R(k)ΩFIG. 13. Differential resistance R as a function of out-of-planemagnetic field B and IR for IU = 0. Orange and white dashedlines show the maximum of the switching and retrapping currents,respectively.reflections [44]. Each resistance map is plotted as a functionof the two voltages that were measured simultaneously.APPENDIX E: EXTENDED SCD DATAAs mentioned earlier, we performed the SCD mea-surements simultaneously for two different junctions.Figures 12(a) and 12(b) show additional SCDs for junction1, while Figs. 12(c) and 12(d) show the SCDs measured forjunction 2. As mentioned in the main text, we observe similartendencies for all investigated SCDs in the range of 0 µA