Thermodynamics in topological Josephson junctions

We study the thermodynamic properties of topological Josephson junctions using a quantum spin Hall (QSH) insulator-based junction as an example. In particular, we propose that phase-dependent measurements of the heat capacity offer an alternative to Josephson-current measurements to demonstrate key topological features. Even in an equilibrium situation, where the fermion parity is not conserved, the heat capacity exhibits a pronounced double peak in its phase dependence as a signature of the protected zero-energy crossing in the Andreev spectrum. This double-peak feature is robust against changes of the tunneling barrier and thus allows one to distinguish between topological and trivial junctions. At short time scales fermion parity is conserved and the heat capacity is $4\pi$-periodic in the superconducting phase difference. We propose a dispersive setup coupling the Josephson junction to a tank LC circuit to measure the heat capacity of the QSH-based Josephson junction sufficiently fast to detect the $4\pi$-periodicity. Although explicitly calculated for a short QSH-based Josephson junction, our results are also applicable to long as well as nanowire-based topological Josephson junctions.

Intriguingly, thermal properties, such as the heat cur- rent driven across a Josephson junction, also depend on ϕ and can provide additional insight to distinguish between trivial and topological junctions [55,56]. Furthermore, it has been shown (employing Hill thermodynamics) that topological phases exhibit unexpected universalities in their thermodynamic signatures [57,58]. Motivated by this and by advances in the understanding of the thermodynamics of non-topological superconductors [59][60][61][62][63][64][65][66][67][68][69][70][71], we have recently applied these thermodynamic concepts to topological Josephson junctions to propose a topological Josephson heat engine [72]. In the same spirit, we discuss another way of utilizing coherent thermodynamics to demonstrate the peculiar nature of topological Joseph-son junctions: We show that measurements of the phasedependent heat capacity can provide distinct signatures originating from the topological ABS and thus represent an alternative property which can be investigated in topological Josephson junctions. Furthermore, we explicitly propose a dispersive measurement scheme, that is, a scheme coupling the Josephson junction to a tank LC circuit via whose response one can detect the characteristic features of topological Josephson junctions experimentally, both with and without parity constraints. As a proof-of-concept we consider a Josephson junction based on a quantum spin Hall (QSH) insulator in this paper, where the inclusion of a magnetic flux in the normal region allows for an additional, experimentally relevant, tuning parameter. The corresponding model and its ABS and continuum states are introduced in Sec. II. In Secs. III and IV, the thermodynamic properties of this system are discussed without and with constraints on the fermion parity, respectively. This is followed by a discussion on possible experimental realizations of measurements of the heat capacity in Sec. V. A short summary concludes the manuscript in Sec. VI.

II. MODEL, ANDREEV BOUND STATES AND CONTINUUM STATES
We consider a short, topological Josephson junction based on a quantum spin Hall (QSH) insulator, where the length L N of the normal (N) region is small compared to the Josephson penetration depth. In this setup, a QSH insulator is partially covered by s-wave superconductors, which proximity-induce pairing to the QSH edge states via tunneling and thus define the superconducting (S) regions [see Fig. 1]. The QSH system lies in the xy plane, with the direction of the superconducting phase bias denoted as the x direction. The Fermi level is situated within the bulk gap, where only edge states exist. Moreover, we assume that the two edges of the QSH system are separated by a distance W S large enough so that there is no overlap between states from opposite edges. In this case, the spin projection in z direction, s =↑ / ↓≡ ±1, and σ = t/b ≡ ±1, describing the top and bottom edges, are good quantum numbers. The corresponding Bogoliubov-de Gennes (BdG) Hamiltonian readŝ where τ j (with j = x, y, z) denote Pauli matrices of the particle-hole degrees of freedom (Appendix A).
Here we study a short junction with a N region of length L N , which we describe by a δ-like N region with h(x) = L N δ(x) [73]. In this model, the proximity-induced pairing amplitude ∆ [74] follows the constant profile ∆(x) = ∆ and the phase convention is describes the superconducting phase difference of the top and bottom edges, respectively. Furthermore,p x denotes the momentum operator, v F the Fermi velocity of the edge states, and V 0 is the potential difference between the N and S regions. The presence of an (orbital) out-of-plane magnetic field B = Be z with B ≥ 0 provides an additional tuning parameter and has two important consequences [53]: ϕ t and ϕ b differ, ϕ t/b = ϕ ± πBW S L N /Φ 0 , where ϕ is the superconducting phase difference at B = 0 and Φ 0 the magnetic flux quantum [75]. Likewise, B induces a Doppler shift described by the Cooper pair momentum p S = π BW S /Φ 0 . Employing a scattering approach, we determine the ABS and the continuum spectrum of Eq. (1) as detailed in Appendices B 1 and B 2.
This approach yields the ABS energies σ s (ϕ σ ) = s −σsgn sin which describes two ABS per edge. Without loss of generality, we will from now on focus on the top edge, where the phase dependence of all quantities is governed by Similar to Ref. [37], where a closely related system has been discussed, for finite p S the dispersion (2) exhibits discontinuities at ϕ t = 2πn with n being an integer (Fig. 2). At these discontinuities, the ABS merge into the continuum states of Eq. (1) described by the continuum density of states (DOS) of the top edge, While these continuum states reduce the effective superconducting gap, the ABS residing in this gap still exhibit a protected crossing [Figs. 2(a,b)] as long as |v F p S | < 2∆. These crossings also correspond to changes in the ground-state fermion parity [33,46]. If |v F p S | > 2∆, there are bound states, but these are all coexisting with continuum states as the superconducting gap is closed [ Fig. 2(c)].

III. THERMODYNAMICS OF A SINGLE EDGE WITHOUT PARITY CONSTRAINTS
From the ABS and continuum spectra of Eq. (1) one can calculate the free energy of the top edge, which then allows one to determine other thermodynamic quantities. First, we consider a situation where there are no parity constraints as is, for example, the case in Josephson junctions strongly coupled to external reservoirs. Then, the free energy is-up to additive ϕ t -independent contributions-given by [33,76] where k B is the Boltzmann constant and T the temperature. Equation (4) is composed of a first part generated by the discrete ABS and a second part due to the continuum DOS ρ t c ( , ϕ t ) = ρ t 0 ( , ϕ t ) + ρ S ( ), consisting of the ϕ t -dependent contribution from the Josephson junction given by Eq. (3) and a ϕ t -independent term originating from the superconducting electrodes, Here the energy scale E S = v F /L S is related to the total length L S of the superconducting QSH edge [33]. Equation (4) allows us to calculate the Josephson current via [76] where e is the elementary charge, and the entropy via From S t one can subsequently obtain the heat capacity which is directly accessible experimentally by measuring the temperature response to a heat pulse. Intriguingly, Eqs. (6) and (7) imply that the phase dependence of S t can be determined from I t via a Maxwell relation [69][70][71] (Appendix C). Without parity constraints, the only phase-dependent contributions to F t 0 arise from the ABS and ρ t 0 . Hence, we can ignore ρ S here. Inserting Eqs. (2) and (3) into F t 0 , it is convenient to split the Josephson current, I t = I t ABS + I t c , calculated from Eq. (6) into an ABS contribution, and a contribution from the continuum, For B = 0, p S = 0 and Eq. (10) vanishes, consistent with the expectation that in short junctions the current is driven by the ABS [76]. Then, Eq. (9) yields the wellknown expression for I t [37]. Only if p S = 0, there is a contribution from Eq. (10) as shown in Fig. 3(a), where I t is presented as a function of ϕ t for different Doppler shifts p S . For 0 < |v F p S | < 2∆, the current-phase relationship is non-sinusoidal and qualitatively similar to the case of p S = 0, apart from a phase shift proportional to p S . As p S is increased further, |v F p S | > 2∆, I t becomes more sinusoidal. In this limit, in fact, the ABS contribution I t ABS is unidirectional [37], while the addition of the continuum contribution due to ρ t 0 results in the sinusoidal behavior of I t seen in Fig. 3(a). For more details on the different contributions to I t , we refer the reader to Appendix D.
In addition to I t , we can also compute C t via Eqs. (7) and (8). We approximate ∆(T ) ≈ ∆(T = 0) and hence ∂∆/∂T ≈ 0 when computing C t , which is reasonable if T T c , where T c is the critical temperature of the parent superconductor. Differing from I t , the exact values of S t and C t at a given ϕ t also require knowledge of the ϕ tindependent contributions to the free energy [72]. To overcome this difficulty, we calculate the difference of C t  with respect to its value at ϕ t = 0, thereby canceling the offset due to the ϕ t -independent contributions. As before, it is convenient to split C t = C t ABS + C t c into contributions from the ABS and the continuum, respectively.
We present the ϕ t dependence of C t for different values of p S in Fig. 3(b). If p S = 0, a double-peak feature described by Eq. (11) develops around the zero-energy crossing of the ABS at ϕ t = π. Consequently, this feature remains also for 0 ≤ |v F p S | < 2∆ as a hallmark of the protected ABS crossing. When |v F p S | > 2∆, the ABS no longer cross at zero energy as shown in Fig. 2(c) and thus C t no longer exhibits the double peak. The sharp double-peak feature of C t is a direct consequence of the form of Eq. (11), which can be shown to always possess a pronounced minimum at t ↑ (ϕ t ) = 0, and thus provides a signature of any zero-energy crossing in the ABS. Since any short 1D topological Josephson junction exhibits an Andreev spectrum qualitatively similar to the one in Fig. 2(a), that is, a protected zero-energy crossing without additional ABS, the phase dependence of the heat capacity for these junctions is also qualitatively similar to the one shown in Fig. 3(b) [77].
The double-peak feature of C t thus allows one to distinguish between topological and trivial junctions: For a phase difference of π, the ABS of a trivial junction exhibit a gap around zero energy, the size of which is controlled by the potential difference V 0 between the S and N regions (which in turn controls the transparency of the junction). At low T , this gap, in turn, gives rise to a single peak of the heat capacity (Appendix E). On the other hand, V 0 cannot remove the zero-energy crossing of the ABS in a topological junction, as can be seen from Eq. (2), which like Eq. (3) does not depend on V 0 [78]. Hence, tuning V 0 and monitoring the phase dependence of the heat capacity always yields a double peak in topological junctions, whereas only a single peak arises in trivial junctions as V 0 increases.
While modulating V 0 thus provides a route to test the protected crossing via C t , the above quantities all exhibit a 2π-periodicity with ϕ t . In order to observe a 4π-periodicity in the thermodynamic quantities, external constraints on the fermion parity have to be imposed.

IV. THERMODYNAMICS OF A SINGLE EDGE WITH PARITY CONSTRAINTS
As long as |v F p S | < 2∆, there is a crossing and one can define the ground-state parity of the top edge by In defining Eq. (13), we have chosen the convention that P (ϕ t = 0) = +1. If the fermion parity is kept constant, the free energy acquires an additional contribution and reads [33,46] with p corresponding to the lower (p = +1) and upper (p = −1) branches at ϕ t = 0, respectively. Again, we omit additive ϕ t -independent contributions to F t p in Eq. (14). The contribution originates from the superconducting electrodes and can be approximated by We now have to evaluate Eq. (14) using the ABS and ρ t c . The resulting expression for F t p is in general quite cumbersome and we do not provide it here explicitly. For the case of p S = 0, however, the problem is simplified significantly because ρ t 0 → 0 and The different thermodynamic quantities under parity constraints can again be obtained from Eqs.
For the corresponding Josephson current, we obtain If J S → 0, which is, for example, the case for T → 0, we recover I t p = p(e∆/2 ) sin(ϕ t /2) [33,45]. Due to the T dependence of J S , the expressions for S t p and C t p turn out to be very cumbersome even for p S = 0.
As F t p exhibits a 4π-periodicity with ϕ t , I t p and C t p inherit this periodicity as illustrated by Figs. 3(c,d), where we have chosen the parity branch p = +1. This 4πperiodicity is in contrast to the behavior of the thermodynamic quantities without parity constraints (Sec. III). Therefore, measurements of, for example, the heat capacity offer additional possibilities to confirm the 4πperiodicity in ϕ t . As shown in Figs. 3(c,d) this 4πperiodicity can also be observed for finite |v F p S | < 2∆. For |v F p S | > 2∆, on the other hand, we can no longer impose any parity constraints on the ABS.
We conclude our discussion of the thermodynamic properties of a single edge by mentioning that we have provided expressions for the top edge in Secs. III and IV. The corresponding expressions for the bottom edge are given by Eqs.

V. EXPERIMENTAL REALIZATIONS
Finally, we discuss a potential experimental realization to measure the heat capacity C t (ϕ t , T ) of one (top) QSH edge of the topological junction. Measuring small heat capacities requires a detection scheme similar to lowenergy calorimetry [79][80][81][82][83][84]. As the measurement should be conducted quite fast (see below), we also take inspiration from fast thermometry techniques [85,86]. Our basic idea here is to subject the topological Josephson junction to radiative heating [87][88][89]: An incident electromagnetic radiation with photon energies ω is absorbed by the junction. When ω < ∆ S , where ∆ S is the gap of the parent superconductor S, the radiation is mainly absorbed at the proximitized single-edge Josephson junction as shown in Fig. 4. For sufficiently short Josephson junctions, such as considered here, the time of flight τ f = L/v F could be very short, of the order of 10 −11 s, if the total junction length L is a few micrometers. Hence, for time scales t τ f , one can define a temperature T QSH for the junction.
We assume a single-edge topological Josephson junction irradiated with a low-intensity flux of photons such that the time between single-photon absorption events is longer than the dead time, that is, the time required for the junction to be restored to its initial state. During the absorption of the single-photon energy ω one observes an increase δT in the temperature such that the final temperature becomes T f QSH = T + δT , with the bath temperature T assumed to be the equilibrium temperature before absorption. The increase of the temperature can be estimated as [84] T +δT For sufficiently small photon energies ω, it can be approximated as δT ≈ ω/C t (ϕ t , T ), with the increase of temperature inversely proportional to the thermal capacity. In Fig. 4(b), we represent the behavior of T QSH after absorption of the photon energy ω at t = t 0 . The junction relaxes back to the bath temperature T with a characteristic time scale τ T (which here plays the role of the dead time) [see Fig. 4 . This time scale is mainly determined by the strengths of the different mechanisms which restore thermal equilibrium such as electron-phonon coupling, radiative losses and heat leakage through the lateral S banks. In the relaxation dynamics, the photon energy ω, which is initially converted fast to thermal energy, is progressively lost to the environment due to thermal losses. Similarly to what happens during the discharging of a capacitor due to a load in an electrical analogy, the thermal response time can also be connected to the thermal capacity [84] via τ T ≈ C t (ϕ t , T )/G th , where G th is the thermal conductance describing all the different thermal losses in the junction. Although this appears to show another way to measure the heat capacitance, it is hard to think of measuring the thermal capacity from τ T since this requires detailed knowledge of the thermal losses G th of the topological Josephson junction, something that-to the best of our knowledge-still needs to be investigated theoretically and experimentally [90]. The relaxation mechanism could become quite slow however. Indeed, in a Josephson junction the thermal losses are reduced due to the thermal opacity of the superconductors, that is, G th is expected to be small. This time scale could become even longer than the time scale for quasiparticle poisoning (see below). Hence, we propose to measure the thermal capacitance only using the first method we have discussed.
In order to make the measurement of δT , one has to perform a fast and precise measurement of the temperature T QSH on time scales t < τ T . Taking inspiration from fast thermometry techniques [85,86], we propose the setup shown in Fig. 4(a). The topological Josephson junction is inserted in a superconducting ring where the superconducting phase difference is mainly driven by the flux Φ [91]. This resembles the configuration of an rf-SQUID for the topological junction. The superconducting ring with inductance L S is coupled with a mutual inductance M to a tank LC circuit. The resonant frequency of the LC circuit, f = (2πL T C T ) −1 , is affected by the coupling to the superconducting ring. In particular, the effective tank circuit inductanceL T becomes [92] renormalizing the intrinsic tank inductance L T because of the superconducting ring, where the inductance is dominated by the kinetic inductance L K (Φ, T ) of the topological Josephson junction [see Fig. 4(a)] since L K L S . The kinetic inductance of the Josephson junction is where I t (ϕ t , T ) is the Josephson current of the junction computed from Eq. (6) as in the analysis in Secs. III and IV. In this setup, the temperature measurement is done in a dispersive way, that is, the LC tank circuit resonance can be investigated with a small rf-signal which probes the resonance. Indeed, any change of the junction temperature T QSH will affect the Josephson current I t (ϕ t , T QSH ), the kinetic inductance L K and finally the resonance frequency f T of the tank circuit. Following a calibration stage where the characteristics of the resonance spectrum are investigated in terms of the temperature of the junction, one can define the optimal strategy and setup for fast measurements of the junction temperature [93]. Now, one needs to relate the measurable temperature T QSH and its radiation-induced change δT to the heat capacity of the QSH edges. It is useful to provide a rough estimate of the expected temperature increase in the proposed setup. We first assume a realistic value for the proximitized gap of the topological junction ∆ ≈ 40 µeV, corresponding to 500 mK [94,95], and a bath temperature of T = 0.1∆, corresponding to around 50 mK, which represents exactly the case considered in Secs. III and IV. For a (proximitized or intrinsic) superconductor of gap ∆ the thermal capacitance at low temperature can be written as C S (T ) ≈ k B N 0 √ 2π(∆/k B T ) 3/2 e −∆/k B T ∆ with N 0 the superconductor density of states [70,96]. This analytical form shows that the thermal capacity is exponentially suppressed with decreasing temperature as ∝ e −∆/k B T . Thus, even including the contribution from the parent superconductors taken as Al films with a superconducting gap of 200 µeV, the phase-independent term of the thermal capacity hardly surpasses a few k B at such low temperatures. Looking at Fig. 3, one can discern that the phase-dependent term of the thermal capacitance, δC t (ϕ t , T ) = C t (ϕ t , T ) − C t (0, T ), contributes with a comparable amount ∼ k B . Assuming photon absorption in the range of ω ∼ 2 GHz, one finds a temperature increase of δT ∼ 10 − 100 mK, having assumed a thermal capacitance of just C t (ϕ t , T ) ∼ 1 − 10k B . These temperature differences are appreciable with respect to the fixed thermal bath temperature T . It is crucial in this setup to observe that the flux Φ fixes the superconducting phase difference between the Josephson junction, fixing ϕ t , and that the temperature calibration has to be done for every value of Φ. This technique can thus be utilized to study the full phase dependence of C t (ϕ t , T ) with its characteristic features discussed in Secs. III and IV, for standard and constrained thermodynamics, respectively.
Finally to measure the 4π-periodic heat capacity of a system with fermion parity constraints, the measurement of C t (ϕ t , T ) needs to be done faster than the quasiparticle poisoning time. With typical quasiparticle poisoning times of the order of 1 µs [97,98] and a resonant frequency of the tank circuit of a few GHz, we expect a temperature detection on the scale of ns or even less. Such time scales can in principle be reached as demonstrated in qubit and SQUID technology [99]. Hence, even measuring a 4π-periodic heat capacity appears feasible [72]. However, even without conserving the groundstate fermion parity, the phase dependence of the heat capacity exhibits pronounced signatures originating from the topological nature of the Josephson junction studied here (see Sec. III).

VI. CONCLUSIONS
In this work, we have analyzed key properties of the thermodynamics in topological Josephson junctions based on quantum spin Hall insulators.
The 4πperiodicity of the free energy in the phase bias ϕ is replicated in other thermodynamic observables when fermion parity is conserved. As a consequence, measuring the ϕ dependence of the heat capacity offers an intriguing alternative to confirm the topological nature of the Andreev bound states. If, on the other hand, there is no fermion parity conservation, the thermodynamic observables will exhibit only a 2π-periodicity in ϕ. Even in this case, however, the protected zero-energy crossing of the Andreev spectrum manifests itself in a sharp doublepeak feature in the ϕ dependence of the heat capacity. Since this double-peak feature is robust against changes in the strength of the tunneling barrier only in a topological junction, one can distinguish between topological and non-topological Josephson junctions by modulating the tunneling-barrier strength, for example, via a gate on top of the normal region.
The inclusion of an out-of-plane magnetic flux inside the normal region induces a Doppler shift and allows for an additional tuning parameter to control the protected zero-energy crossing in the Andreev spectrum: While this Doppler shift reduces the effective gap of the Andreev spectrum, the crossing remains as long as this gap is not completely closed. The gap closing and disappearance of the zero-energy crossing in the Andreev spectrum with increasing magnetic flux is also reflected in the heat capacity. Although our results have been obtained for a quantum-spin-Hall-based junction, they are also applicable to other topological Josephson junctions, such as nanowire-based junctions.
Finally, we have also provided a dispersive setup to measure the small thermal capacity at fixed phase bias which is sufficiently fast to detect the 4π-periodicity in the absence of quasiparticle poisoning. Comparison with the 2π-periodic equilibrium heat capacity also provides a method to indirectly investigate thermal losses and realize novel calorimetric applications.

Appendix A: Simplified Hamiltonian
In this section, we provide additional details on the derivation of the simplified BdG Hamiltonian (1), which serves as the starting point of our calculations. With the basis order ψ t,↑ ,ψ t,

the BdG Hamiltonian describing the Josephson junction introduced in Sec. II is given bŷ
Here s j , σ j and τ j (with j = x, y, z) denote Pauli matrices describing spin, top/bottom edge and particle-hole degrees of freedom, respectively. Note that unit matrices are not written explicitly in Eq. (A1). The N region is described by the profile h(x) and the potential V 0 , which can also be viewed as describing the difference between the chemical potentials in the S and N regions, µ S and µ N = µ S − V 0 . Here we use a δ-barrier model with h(x) = L N δ(x) and ∆(x) = ∆. A more general approach to the junction would be to use a finite N region with h(x) = Θ(L N − x)Θ(x) and ∆(x) = ∆ [Θ(x − L N ) + Θ(−x)]. Such a model is capable of describing both a short as well as a long junction. As explained in Appendix B below, the δ-barrier model is still suitable to capture the essential physics of a short junction based on QSH edge states. Since Ĥ BdG , s z = Ĥ BdG , σ z = 0, the wave functions can be described by the good quantum numbers s =↑ / ↓≡ ±1 for the spin projection in z direction and σ = t/b ≡ ±1 for the top and bottom edges. Hence, an eigenstate Ψ s,σ (x) of the BdG Hamiltonian (A1) can be written as Ψ s,σ (x) = ψ s,σ (x) ⊗ η σ ⊗ χ s , which satisfieŝ whereĤ s,σ is given by Eq. (1) in the main text. The spinors χ s and η σ are the eigenvectors of the Pauli matrices s z and σ z , respectively. They satisfy s z χ s = sχ s and σ z η σ = ση σ .

Appendix B: Computing the Andreev and continuum spectra
In order to compute the ABS as well as the ϕ-dependent part of the continuum DOS ρ c ( , ϕ), we employ a wave-function matching approach to solve Eq. (A2). For the following derivations, it is convenient to introduce the spin-dependent energy variable s = − sv F p S /2.

Andreev bound states
For ABS, that is, states with | s | < ∆ and thus states localized around the N region, we can make the piecewise ansatz for a junction with a finite N region. The ansatz for a δ-junction is similar to Eq. (B1), with the states ψ s,σ (x < 0) given by the first line of Eq. (B1) and ψ s,σ (x > 0) given by the third line of Eq. (B1). In Eq. (B1), C (s) , and L tot = L N + L S is the unit length of the entire S/N/S edge. The ansatz (B1) for ψ s,σ (x) has been chosen in such a way that Eq. (A2) is satisfied in each S and N region separately and that lim x→±∞ ψ s,σ (x) = 0. The coefficients a s,σ , b s,σ , e s,σ , and h s,σ have to be determined from the boundary conditions at the S/N interfaces.
For the δ-barrier, the boundary condition can be obtained by integrating Eq. (A2) from x = −η to x = η with η → 0 + . The corresponding procedure [55,[100][101][102] yields where Z 0 = V 0 L N /( v F ). For a finite barrier, on the other hand, the boundary conditions require ψ s,σ (x) to be continuous at the S/N interfaces, First, we consider the δ-barrier model, for which we invoke the boundary condition (B2) and require a nontrivial solution for the coefficients a s,σ and b s,σ . This procedure yields a transcendental equation for the energy variable s in the form of where n is an integer that has to be chosen in such a way that 0 ≤ arccos ( s /∆) ≤ π. Solving Eq. (B4) and using s = − sv F p S /2, we obtain the ABS given by Eq. (2) in the main text. For a finite N region, we use the boundary conditions (B3) and require a nontrivial solution for the coefficients a s,σ , b s,σ , e s,σ , and h s,σ . Now, we obtain the transcendental equation where n is again an integer chosen in such a way that 0 ≤ arccos ( s /∆) ≤ π. The second term on the left-hand side of Eq. (B5) takes into account the finite width of the N region. This term becomes only important in long junctions, ∆ v F /L N , where it causes multiple subbands to appear in the Andreev spectrum. These additional bound states in long junctions cannot be captured with a δ-model, which always yields two ABS per edge. Still for short junctions, ∆ v F /L N , Eq. (B5) tends to Eq. (B4).

Continuum states
In the previous section, we have determined the discrete Andreev spectrum. If we now turn to the continuous spectrum with s > ∆, we have to modify the ansatz (B1) to account for propagating wave functions. An incident electron-like quasiparticle propagating from x → −∞ to x → ∞ with s > ∆ can be described by the ansatz if a finite N region is considered. The ansatz for a δ-junction is again similar to Eq. (B6), but with ψ s,σ (x < 0) given by the first line of Eq. (B6) and ψ s,σ (x > 0) given by the third line of Eq. (B6). In order for the quasiparticle to be a right mover, sσ = 1 in Eq. (B6), that is, s =↑ at the top edge and s =↓ at the bottom edge. Corresponding equations can be set up for hole-like quasiparticles propagating to the right as well as for quasiparticles propagating to the left. In Eq. (B6), u 2 , and all other quantities are the same as defined in Eq. (B1) [103].
By invoking the boundary conditions (B4) or (B5), one can then obtain the reflection and transmission coefficients r s,σ eh and t s,σ ee in Eq. (B6). For hole-like incident quasiparticles, one can likewise obtain r s,σ he and t s,σ hh . These states then allow us to set up the S matrix as and The scattering matrix S σ SN S ( , ϕ) allows us then to compute the (ϕ-dependent part of the) continuum DOS of the Josephson junction as Since S SN S is block-diagonal in the top/bottom-edge degree of freedom, det (S SN S ) = det (S t SN S ) det S b SN S can be decomposed in a contribution from the top edge (σ = 1) and one from the bottom edge (σ = −1). Then, we obtain If we consider a finite N region, the above discussion still applies, but now for the top edge if ↑/↓ > ∆. The contribution from the bottom edge can be obtained by replacing ϕ t → ϕ b , ↑ → ↓ , and ↓ → ↑ .

Appendix C: Free energy and thermodynamic observables
From the ABS and continuum spectra of Eq. (A2) one can calculate the free energy F σ (ϕ, T ) at a single edge σ = t/b. These free energies can in turn be used to obtain the Josephson current and other thermodynamic quantities of a given edge.

No parity constraints
Without parity constraints, the free energy is-up to an additive ϕ-independent contribution-given by [33,76] where k B is the Boltzmann constant and T the temperature. The sum over n describes the contribution from the ABS with discrete energies σ n (ϕ), while the integral describes the contribution from the continuum states with DOS ρ σ c ( , ϕ). At this point, it is important to note that due to the lack of spin degeneracy the degeneracy factor is just g = 1 in Eq. (C1) [33]. Whereas we can directly insert the Andreev spectrum σ n (ϕ) in Eq. (C1), we have to compute the continuum DOS The ϕ-dependent continuous spectrum of the Josephson junction is described by ρ σ 0 ( , ϕ), which we compute from the scattering matrix via Eq. (B10). Following Ref. [33], we also include in Eq. (C2) a ϕ-independent term originating from the superconducting electrodes, see also Eq. (5) in the main text. Here the energy scale E S = v F /L S is determined from the length L S of the superconducting electrodes. While ρ S ( ) is independent of ϕ and does not affect the Josephson current in the absence of parity constraints [76], it can play an important role if the parity is kept constant.

Parity constraints
If the fermion parity p = ± is conserved at a given edge σ = t/b, the free energy acquires an additional contribution due to the parity constraint and reads [33,46] (C4) Here p = ±1 and the function P (ϕ) describes the ground-state fermion parity as ϕ is tuned, see also Sec. IV in the main text. From the form of Eq. (C4) it is clear that even ϕ-independent terms in ρ σ c are important in determining the ϕ dependence of F σ p . This is why, for example, the contribution from the superconducting electrodes can also affect other quantities such as the Josephson current. If we split the integral over ρ σ c in Eq. (C4) into an integral over ρ σ 0 and ρ S , define and insert the ABS for σ n , we arrive at Eq. (14) in the main text.

Josephson current and entropy at a single edge
Having determined F σ , given either by Eq. (C1) or by Eq. (C4), we can then calculate the Josephson current via where e is the elementary charge, as well as the entropy via

Due to the Maxwell relation
the phase dependence of the entropy can be determined from the Josephson current as Thus, one can use the temperature dependence of I σ to obtain the ϕ dependence of S σ . In addition, the heat capacity provides a quantity that is directly accessible experimentally.
Appendix D: Decomposing the Josephson current and heat capacity into contributions from the Andreev bound states and the continuum As mentioned in the main text, one can separate the contributions of the ABS and the continuum states to the free energy F t 0 from each other if there are no parity constraints. Consequently, one can also split the Josephson current, the entropy and the heat capacity in two such separate contributions [see Eqs. (9), (10), (11), and (12) in the main text].
For illustration, Fig. 5 shows the separate contributions I t ABS and I t c as well as the full Josephson current I t = I t ABS + I t c for different Doppler shifts v F p S . Whereas in Figs. 5(a,b) with p S = 0 the total current is given by I t = I t ABS , the situation changes for p S = 0: For finite p S , some weight of the DOS is transferred from the ABS to the continuum states [104,105], and a finite I t c develops. If 0 < |v F p S | < 2∆, the effect of the continuum states is to mainly shift the total current I t compared to I t ABS as shown in Figs. 5(c,d). As p S is increased further, such that |v F p S | > 2∆, I t ABS becomes unidirectional [see Fig. 5(e)], in the sense that I t ABS flows in the same direction for any phase difference ϕ t . This is similar to Ref. [37], where the unidirectional, chiral current carried by the ABS has been shown as a hallmark of a 4π-periodic current-phase relation. The total Josephson current, on the other hand, is a smooth function of ϕ t and no longer unidirectional due to the contribution from the continuum states. These total Josephson currents are the ones shown in Fig. 3(a) of the main text. Similarly, the total heat capacity of the top edge computed via Eqs. (7) and (8) can also be split into contributions C t ABS from the ABS and contributions C t c from the continuum [see Eqs. (11) and (12) in the main text]. Figure 6 shows C t ABS and C t c as well as the full heat capacity C t = C t ABS + C t c for different values of p S . Assuming ∆(T ) ≈ ∆(T = 0), we present C t ABS and C t c up to a ϕ t -independent constant, while the total heat capacity of the top edge C t is measured with respect to its value at ϕ t = 0. For v F p S = 0, the heat capacity in δ-junctions is again exclusively due to the ABS as shown in Figs. 6(a,b). Importantly, one can see that around the zero-energy crossing of the ABS, ϕ t = π, a sharp double-peak feature given by Eq. (11) develops.
For finite p S , the continuum contribution plays a crucial role: As one can discern from Figs. 6(c,e), C t ABS is discontinuous due to the discontinuities of the ABS spectrum (2) if ϕ t is an integer multiple of 2π. These apparent discontinuities are artificial, however, since the ABS merge into the continuum at these values of ϕ t . Indeed, including the continuum contribution (12) lifts these discontinuities in the total heat capacity C t as demonstrated by Figs. 6(d,f). If 0 < |v F p S | < 2∆, C t is qualitatively very similar to the case of p S = 0 with a sharp double-peak feature around the phase ϕ t , where the ABS cross at zero energy. If |v F p S | > 2∆, the ABS no longer cross at zero energy as shown in Fig. 2(c) and consequently C t no longer exhibits the sharp double-peak feature characteristic of the zero-energy crossing [see (Color online) Phase dependence of (a,c) the relative heat capacity δC t (ϕt) = C t (ϕt) − C t (0) and (b,d) C t (ϕt) = d 2 C t (ϕt)/dϕ 2 t for the top edge of a topological Josephson junction and for a trivial Josephson junction with different transmission probabilities τtriv. Here pS = 0 and no parity constraints have been imposed on the topological junction. The temperature has been chosen as kBT = 0.1∆ in panels (a,b) and as kBT = 0.02∆ in panels (c,d).

Appendix E: Comparison between the heat capacities of a topological and a trivial Josephson junction
In the absence of parity constraints, the heat capacity C t (ϕ t ) of a topological Josephson junction exhibits pronounced double-peak features as a function of the phase bias ϕ t for k B T ∆, providing a signature of the protected zero-energy crossing. Here we show that this is in contrast to the behavior of the heat capacity in a trivial Josephson junction. For simplicity, we compare the topological Josephson junction discussed in the main text without parity constraints and for p S = 0 to a one-dimensional (1D) trivial Josephson junction.
We consider a short 1D trivial Josephson junction based on a quadratic Hamiltonian. Such a junction exhibits a pair of ABS per spin with energies [76] ± triv (ϕ t ) = ±∆ 1 − τ triv sin 2 ϕ t 2 , where 0 ≤ τ triv ≤ 1 describes the junction transparency (modulated, for example, by a potential difference V 0 between the N and S regions). Since we are studying the limit of short junctions, we assume that the ϕ t dependence of the free energy originates solely from the ABS. Then, the heat capacity (per spin) of the trivial junction is also described by Eq. (11) in the main text, with t ↑ (ϕ t ) replaced by + triv (ϕ t ). For a perfectly transparent junction, τ triv = 1 and Eq. (E1) reduces to ± triv (ϕ t ) = ±∆ cos(ϕ t /2), similar to the case of the ABS in a short topological junction for p S = 0, see Eq. (2) in the main text. In contrast to a topological junction, the zero-energy crossing at ϕ = π is not protected in a trivial junction and is removed if τ triv = 1. This results in a gap of δ = 2∆ √ 1 − τ triv . In Fig. 7, we show C t (ϕ t ) (measured with respect to ϕ t = 0) and its second derivative, C t (ϕ t ) = d 2 C t (ϕ t )/dϕ 2 t , for a topological Josephson junction as well as for a trivial junction with different values of τ triv . As long as k B T δ , the splitting δ can be thermally resolved and a peak develops at ϕ t = π for a trivial Josephson junction as illustrated by Fig. 7(a) for τ triv = 0.8 and τ triv = 0.9. If k B T is too large to clearly resolve δ , as is the case in Figs. 7(a,b) for τ triv = 0.99, C t (ϕ t ) resembles that of a topological junction. If the temperature is decreased, however, the doublepeak feature around ϕ t = π merges into a single peak at ϕ t = π as can be seen by comparing Figs. 7(a) and (c). In a topological junction, on the other hand, the double-peak feature with a minimum at ϕ t = π remains as T is decreased. Hence, probing the T dependence of C t (ϕ t ) around ϕ t = π allows to distinguish between topological and trivial Josephson junctions: C t (ϕ t ) always exhibits a minimum at ϕ t = π in a topological junction, while a peak develops at ϕ t = π in a trivial junction as T is decreased.