Thermodynamics of Non-Hermitian Josephson junctions with exceptional points

At weak coupling, an external environment only induces broadening and small shifts to the levels of a quantum system. In contrast, the strong coupling limit is highly nontrivial and gives rise to many interesting concepts in e.g. quantum dissipation Weiss (2021), quantum information science Harrington et al. (2022) or quantum thermodynamics Rivas (2020), just to name a few. An interesting example is the emergence of spectral degeneracies in the complex spectrum (resulting from integrating out the environment), also known as exceptional point (EP) bifurcations, where eigenvalues and eigenvectors coalesce Berry (2004). During the last few years, a great deal of research is being developed in so-called non-Hermitian (NH) systems with EPs, in various contexts including open photonic systems Doppler et al. (2016), Dirac Zhen et al. (2015), Weyl Cerjan et al. (2019) and topological matter in general Shen et al. (2018); Gong et al. (2018); Kawabata et al. (2019); Bergholtz et al. (2021).

The role of NH physics and EP bifurcations in systems with Bogoliubov-de Gennes (BdG) symmetry has hitherto remained unexplored, until recently. Specifically, there is an ongoing debate on how to correctly calculate the free energy in open BdG systems, a question relevant in e.g Josephson junctions coupled to external electron reservoirs. Depending on different approximations, such "NH junctions" have been predicted to exhibit exotic effects including imaginary persistent currents Zhang et al. (2022); Li et al. (2021) and supercurrents Cayao and Sato (2023); Li et al. (2024) or various transport anomalies at EPs Kornich (2023). If one instead uses the biorthogonal basis associated with the NH problem Shen et al. (2024), or an extension of scattering theory to include external electron reservoirs Beenakker (2024), the supercurrents are real and exhibit no anomalies. At the heart of this debate is whether a straightforward use of the complex spectrum plugged into textbook definitions of thermodynamic functions leads to meaningful results or whether, on the contrary, NH physics needs to be treated with care when calculating the free energy.

We here present a well-defined procedure, valid for arbitrary coupling and temperature, which allows us to calculate the free energy, Eq. (4), without any divergences at EPs. Derivatives of this free energy, allow us to calculate physical observables such as entropy (6) or supercurrents (7). Interestingly, entropy changes of $\log 2/2$ can be connected to emergent Majorana zero modes (MZMs) at EPs Pikulin and Nazarov (2013); San-Jose et al. (2016); Avila et al. (2019). While such fractional entropy steps were predicted before in seemingly different contexts Smirnov (2015); Sela et al. (2019); Becerra et al. (2023), our analysis in terms of EPs allows us to obtain precise analytical boundaries for the temperatures at which they appear. We propose a novel Maxwell relation connecting supercurrents and entropy, Eq. (19), which would allow the experimental detection of the effects predicted here.

The starting point of our analysis is the description of an open quantum system in terms of a Green’s function

where $H_{\mathrm{eff}}(\omega)=H_{Q}+\Sigma^{r}(\omega)$ is an effective NH Hamiltonian which takes into account how the quantum system $H_{Q}$ is coupled to an external environment through the retarded self-energy $\Sigma^{r}(\omega)$. In what follows, we consider the case where an electron reservoir induces a tunneling rate ¹¹1We consider the simplest case where this selfenergy originates from a tunneling coupling to an external electron reservoir in the so-called wideband approximation where both the tunneling amplitudes and the density of states in the reservoir are energy independent, which allows to write the tunneling couplings as constant Meir and Wingreen (1992)., such that the complex poles of Eq. (1) have a well-defined physical interpretation in terms of quasi-bound states. If $H_{\mathrm{eff}}$ is a BdG Hamiltonian, electron-hole symmetry can be satisfied in two non-equivalent ways: (i) one can have pairs of poles with opposite real parts and with the same imaginary part $\epsilon^{\pm}=\pm E-i\gamma=-(\epsilon^{\mp})^{*}$; or, alternatively, (ii) two independent and purely imaginary poles $\epsilon^{\pm}=-i\gamma^{\pm}=-(\epsilon^{\pm})^{*}$. A bifurcation of the former, corresponding to standard finite-energy BdG modes with an equal decay to the reservoir $\gamma$, into the latter, two MZMs with different decay rates Pikulin and Nazarov (2013); San-Jose et al. (2016); Avila et al. (2019), defines an EP (Fig. 1a, inset).

Calculating the free energy of an open quantum system is nontrivial since a direct substitution of a complex spectrum $\epsilon_{j}=E_{j}-i\gamma_{j}$ in the standard expression $F=-\frac{1}{\beta}\log Z=-\frac{1}{\beta}\sum_{j}\log\left(1+e^{-\beta(\epsilon_{j}-\mu)}\right)$, can lead to complex results and divergences after an EP Cayao and Sato (2023); Li et al. (2024). To avoid inconsistencies, one possibility is to use the occupation $\langle N\rangle=\int_{-\infty}^{\infty}d\omega G^{<}(\omega)=-\frac{\partial F}{\partial\mu}$, which is a well-defined quantity in open quantum systems, even in non-equilibrium situations where the so-called Keldysh lesser Green’s function $G^{<}(\omega)$ can be generalized beyond the fluctuation-dissipation theorem Meir and Wingreen (1992). In BdG language, $\langle N\rangle$ can be written as

where $f(\omega)$ is the Fermi-Dirac function and we have explicitly separated the total spectral function

in its particle ($\real(\epsilon_{j})>0$) and hole ($\real(\epsilon_{j})<0$) branches, $\rho^{p}(\omega)$ and $\rho^{h}(\omega)$, respectively ²²2A global factor of $1/2$ arises from the duplicity of dimensions in BdG formalism.. We have also added a Lorentzian cutoff $\Omega(D)=D^{2}/(D^{2}+\omega^{2})$ to avoid divergences in the thermodynamic quantities as $T\to 0$ (see Appendix B). The integral in Eq. (2) can be analytically solved by residues (Appendix B) and then be used to obtain $F$ as³³3Note that while the above discussion is limited to non-interacting BdG Hamiltonians, the formalism can be generalized to an arbitrary interacting system. Assuming the number of particles is conserved, each eigenstate $\Psi_{N,j}$ has a specific particle number $N$ and energy $E_{N,j}$. The single-particle Green’s function in Lehmann representation will then have a set of simple poles at $\pm(E_{N,j}-E_{(N-1),i})+i\eta$, that is, the difference between the energies of the states $\Psi_{N,j}$ and $\Psi_{(N+1),i}$. Adding to these poles a finite imaginary part coming from the coupling to a reservoir, as long as it only depends on the state, but not on the temperature $T$ or the chemical potential $\mu$, all the dependences with $T$ and $\mu$ will come from the Fermi function, just as for the non-interacting system.

with $\log\Gamma(z)$ being the log-gamma function and $h(T,E_{j},\gamma_{j})$ a generic function coming from the integration. We now perform the limits $\gamma_{j}=0$ and $T\to 0$ of the previous expression,

which, by comparison with well-known limits Zagoskin (2014); Hewson (1993), give $h(T,E_{j},\gamma_{j})=-\gamma_{j}/\pi-T\log(2\pi)$ ⁴⁴4Interestingly, this term is relevant for the entropy since it adds a physical contribution $\log(2\pi)$. In contrast, the term $-\gamma_{j}/\pi$ cancels out in the supercurrent (see Appendix B).. From now on, we fix $\mu=0$ but the complete derivation with full expressions, including $\mu\neq 0$, can be found in Appendix B. Derivatives of Eq. (4) allow us to obtain relevant thermodynamic quantities, as we discuss now.

Using the free energy in Eq. (4), the entropy, defined as $S=-\partial F/\partial T$, reads:

where $\psi(z)$ is the digamma function. From Eq. (6), we can define a crossover temperature $T_{j}=|\epsilon_{j}|/2$ as the inflection point when the eigenvalue $\epsilon_{j}$ begins to have a non-zero contribution ($S_{j}=\log(2)/4$) to the total entropy, Fig. 1a. Hence, two standard BdG poles $\epsilon^{\pm}=\pm E-i\gamma$ will contribute at the same temperature to the entropy of the system since their absolute values are equal, $|\epsilon^{\pm}|=\sqrt{E^{2}+\gamma^{2}}$, and thus a single plateau of $S=\log(2)$ can be measured. On the contrary, after an EP at $\epsilon^{+}=\epsilon^{-}=-i\gamma$, the poles bifurcate taking zero real parts and different imaginary parts, $\epsilon^{+}=-i\gamma^{+}$ and $\epsilon^{-}=-i\gamma^{-}$, and separating from each other a distance $\gamma^{-}-\gamma^{+}$. Then, each pole will contribute to the entropy at a different temperature $T^{\pm}=|\epsilon^{\pm}|/2=\gamma^{\pm}/2$, giving rise to a fractional plateau $S=\log(2)/2$ of width $T^{-}-T^{+}=(\gamma^{-}-\gamma^{+})/2$. Here it is very important to point out that this nontrivial behavior of the entropy seems absent in the free energy that we used for the calculation of $S=-\partial F/\partial T$. Indeed, $F$ is seemingly insensitive to any EP bifurcation even when varying the temperature over six decades, Fig. 1b.

Our method also allows to calculate the supercurrent of any generic NH Josephson junction (or similarly the persistent current through a normal ring Shen et al. (2024)) by just considering a phase-dependent spectrum $\epsilon_{j}(\phi)$ and taking phase derivatives of Eq. (4) as $I=\frac{\partial F}{\partial\Phi}=\frac{2e}{\hbar}\frac{\partial F}{\partial\phi}$, which gives

As $T\to 0$, the supercurrent carried by a pair of BdG poles simply becomes ⁵⁵5For a single pair of BdG poles, $\gamma$ is independent/dependent of the phase before/after the EP and the other way around for $E$ (see Appendix C).

which, for example, allows us to calculate the supercurrent carried by Andreev levels in a short junction coupled to an electron reservoir almost straightforwardly (see Appendix C). Note that, although $I(\phi)$ has a cutoff-dependent term, it cancels by the particle-hole symmetry of the problem (Appendix C). Eqs. (8) strongly differ from a calculation using directly the complex spectrum Cayao and Sato (2023); Li et al. (2024)

Eqs. (4), (6) and (7) are the main results of this paper and allow to calculate thermodynamics from generic open BdG models (arbitrary coupling and temperature) that can be written in terms of complex poles (Eq. (1)).

As an application we now consider a quantum dot (QD) array in a so-called minimal Kitaev model

where $c_{i}^{\dagger}$ ($c_{i}$) denote creation (annihilation) operators on each QD with a chemical potential $\mu_{i}$. The QDs couple via a common superconductor that allows for crossed Andreev reflection and single-electron elastic co–tunneling, with coupling strengths $\Delta_{i,i+1}$ and $t_{i,i+1}$, respectively. Remarkably, only two QDs are enough to host two localized MZMs ⁶⁶6Here, it is important to clarify that MZMs in these minimal chains only appear in fine-tuned “sweet spots” in parameter space and without topological protection, which is why they are often called poor man’s Majorana modes. However, their non-local nature remains, where one MZM localizes in each QD of the minimal chain. Interestingly, longer chains, already at the three QD level, start to show a bulk gap and some degree of protection Bordin et al. (2024) $\eta_{1}$ and $\eta_{2}$ when a so-called sweet spot is reached with $\Delta_{1,2}=t_{1,2}$. This theoretical prediction Leijnse and Flensberg (2012) has recently been experimentally implemented Dvir et al. (2023); ten Haaf et al. (2024a). Let consider now a second double QD array (Majoranas $\eta_{3}$ and $\eta_{4}$) that forms a Josephson junction with the former array with a coupling $H_{\mathrm{JJ}}=-t_{2,3}e^{i\frac{\phi}{2}}c_{2}^{\dagger}c_{3}+\mbox{H.c.}$, with $\phi$ being the superconducting phase difference between both arrays and $t_{2,3}$ the tunneling coupling between inner QDs. If, additionally, the two inner QDs are coupled to normal reservoirs with rates $\gamma_{2}^{0}$ and $\gamma_{3}^{0}$ this system is a realization of a Non-Hermitian Josephson junction containing Majorana modes (see the sketch in Fig. 2a). In the low-energy regime, this model can be described in terms of four Majorana modes Pino et al. (2024). Assuming $\mu_{i}=0$ $\forall i$ and $\Delta_{12}=t_{12}$ and $\Delta_{34}=t_{34}$, only the inner Majoranas $\eta_{2}$ and $\eta_{3}$ are coupled and lead to BdG fermionic modes of energy $\epsilon_{\mathrm{inner}}^{\pm}=-\frac{i}{2}\gamma_{0}\pm\frac{1}{2}\Lambda(\phi)$, where $\gamma_{0}=\gamma_{2}^{0}+\gamma_{3}^{0}$, $\delta_{0}=\gamma_{2}^{0}-\gamma_{3}^{0}$ and $\Lambda(\phi)=\sqrt{2t_{23}^{2}(1+\cos\phi)-\delta_{0}^{2}}$, while the outer modes remain completely decoupled and $\epsilon_{\mathrm{outer}}^{\pm}=0$. For $\delta_{0}=0$, one recovers the standard Majorana Josephson term: $\epsilon_{\mathrm{inner}}^{\pm}=-\frac{i}{2}\gamma_{0}\pm t_{23}\cos(\phi/2)$. When $\delta_{0}\neq 0$, the spectrum develops EPs at phases $\phi_{\mathrm{EP}}=\arccos{[\frac{\delta_{0}^{2}}{2t_{23}^{2}}-1]}$. For an example of the resulting supercurrents and a comparison against Eq. (9), see Fig. 2. Similarly, when $\Delta_{12}=\Delta_{34}=\Delta$, $t_{12}=t_{34}=t$ but $t\neq\Delta$, an additional pair of EPs appears at ⁷⁷7Analytical expressions of eigenvalues, a similar calculation for a junction containing Andreev levels, as well as detailed discussions can be found in the Appendix C.

Using these analytics, the crossover temperatures $T_{j}=|\epsilon_{j}|/2$ read

To illustrate their physical meaning, we plot a full calculation of the entropy $S(T,\phi)$ using Eq. (6), Fig. 3a, together with the analytical expressions in Eq. (12) (solid lines). This plot demonstrates that changes in entropy can be understood from EPs, a claim that is even clearer by analyzing cuts at fixed phase (Fig. 3b). Interestingly, a universal entropy change of $1/2\log 2$ can be linked to the emergence of MZMs at phase $\phi=\pi$ as $T\to 0$. Alternatively, the entropy loss due to Majoranas can be seen in phase-dependent cuts taken at different temperatures, Fig. 3c, which show as an interesting behavior where a $S=2\log 2$ plateau at large temperatures becomes an emergent narrow resonance, centered at $\phi=\pi$ and of height $S=1.5\log 2$, as $T$ is lowered.

While in the minimal Kitaev chain model in Eq. (10) the pairing potential $\Delta$ is assumed to be a temperature-independent parameter, in a real experiment it comes from crossed Andreev reflections mediated by a middle segment separating both quantum dots. Such segment is typically a semiconducting region proximitized by a superconductor such that the parent superconducting gap $\Delta_{0}\gg\Delta$ (in the experiments of Ref.  Dvir et al. (2023); ten Haaf et al. (2024a), the gap of the parent aluminum superconductor is $\Delta_{0}\approx 270\mu$eV $\gg\Delta\approx 10\mu$eV). Thus, a regime where $T\gg\Delta$ without destroying superconductivity is possible. In this section, we include the BCS-like temperature dependence of $\Delta_{0}$ in a microscopic model of a minimal Kitaev chain to fully support this argument.

The DQD model of Eq. (10) in the main text can be realized microscopically by means of a common superconductor-semiconductor hybrid region which couples both QDss. Specifically, they couple via a subgap Andreev bound state (ABS) living in the middle region. Such ABS coupling allows for crossed Andreev reflection (CAR) and single-electron elastic co-tunneling (ECT), with coupling strengths $\Delta$ and $t$, respectively. Physically, the spin-orbit coupling in semiconductor-superconductor provides a spin-mixing term for tunneling electrons. The spin-mixing term can lead to finite ECT and CAR amplitudes for spin-polarized QDs when the spin-orbit and magnetic fields are non-colinear. This ABS provides a low-excitation energy for ECT and CAR, which, starting from a microscopic model of two QDs connected by a semiconductor-superconductor middle region, can be obtained from a Schrieffer-Wolff transformation to obtain effective cotunneling-like terms Liu et al. (2022)

where $t_{L/R}$ are the local tunneling strengths between each QD and the central region and

is the energy of the Andreev state, with $z\equiv\mu_{\mathrm{ABS}}/\Delta_{0}$ being the chemical potential of the middle region, and BdG coefficients given by $u^{2}=1-v^{2}=1/2+\mu_{\mathrm{ABS}}/(2E_{\mathrm{ABS}})$. As discussed in Ref. Liu et al. (2022), this energy dependence of $\Delta$ and $t$ allows one to vary their relative amplitudes by changing the chemical potential $\mu_{\mathrm{ABS}}$ of the ABS and tune them to precise points, so-called "sweet spots", where $\Delta=t$, a precise tuning which has already been demonstrated in different experimental platforms and configurations Dvir et al. (2023); ten Haaf et al. (2024a, b).

In Eqs. (13), the gap $\Delta_{0}$ of the parent superconductor, is assumed to be temperature-independent Liu et al. (2022). We now extend this microscopic theory and explicitly consider the BCS-like temperature dependence of the Aluminium parent gap $\Delta_{0}\rightarrow\Delta_{T}$ which makes the crossed Andreev reflection and elastic cotunneling terms in Eq. (13) temperature-dependent too. Interestingly, their temperature dependence is highly nonmonotonic, with large regions where they remain nearly constant until they decrease near the Aluminium critical temperature $T_{c}\approx 141.651\;\mu\mathrm{eV}$ where the parent gap closes (see Appendix F).

Specifically, we study the temperature dependence of the difference $|t-\Delta|$ (Fig. 4a) for $\mu_{\mathrm{ABS}}=\Delta_{0}$, which corresponds to the sweet spot at zero temperature. If, additionally, each QD is coupled to normal reservoirs with rates $\gamma_{1}^{0}\neq\gamma_{2}^{0}$, this temperature evolution governs the emergence of exceptional point (EP) bifurcations, which are now temperature-dependent (Fig. 4b). As the temperature increases in Fig. 4a, both couplings remain equal $|t-\Delta|=0$ until reaching a temperature $T^{*}$ (lilac regions). Above $T^{*}$ both couplings start to differ $|t-\Delta|\neq 0$. Nevertheless, as long as the condition $|t-\Delta|<|\gamma_{1}^{0}-\gamma_{2}^{0}|/2$ is fulfilled, the system still hosts a pair of MZMs (green region). This can be explicitly seen in Fig. 4b where we plot the full temperature-dependent evolution of the complex spectrum. Specifically, the green region corresponds to a regime where two MZMs ($\real(\epsilon_{\mathrm{outer}}^{+})=\real(\epsilon_{\mathrm{outer}}^{-})=0$) are asymmetrically coupled to the reservoir ($\imaginary(\epsilon_{\mathrm{outer}}^{+})\neq\imaginary(\epsilon_{\mathrm{outer}}^{-})\neq 0$). This finite-temperature regime with MZMs persists until $|t-\Delta|$ is large enough to close the bifurcation. Specifically, the exact value of the temperature at which the EP forms, $T_{\mathrm{EP}}$, is determined by the point where the condition $|t-\Delta|=|\gamma_{1}^{0}-\gamma_{2}^{0}|/2$ is fulfilled

Above $T_{\mathrm{EP}}$, namely $|t-\Delta|>|\gamma_{1}^{0}-\gamma_{2}^{0}|/2$ the system no longer contains MZMs (pink region). The lilac region, on the other hand, should be understood as the temperature $T^{*}$ below which the $T\to 0$ limit is valid

Even though the calculations in Fig. 4 have been performed for the particular (but reasonable) values of $t_{L/R}=2/3\Delta_{0}$ and $|\gamma_{1}^{0}-\gamma_{2}^{0}|=0.1\Delta_{0}$, we can extract some general conclusions that enable the possibility of an experimental demonstration of our claim: although the magnitude of $|t-\Delta|$ depends on the local tunneling strengths $t_{L/R}$, their functional form against the temperature remains the same, reaching a maximum value

at $\Delta_{T}=(\sqrt{2}-1)\mu_{\mathrm{ABS}}$, which is independent of $\gamma_{1/2}^{0}$ and $t_{L/R}$, and corresponds to a large value of $T_{\mathrm{max}}\approx 0.94T_{c}$ when $\mu_{\mathrm{ABS}}=\Delta_{0}$. Conversely, $T_{\mathrm{EP}}$ depends on the ratio between the reservoir coupling asymmetry $|\gamma_{1}^{0}-\gamma_{2}^{0}|$ and the tunneling strengths $t_{L/R}$, and it must be smaller than $T_{\mathrm{max}}$ such that the EP can develop. For the choice of parameters in Fig. 4, $T_{\mathrm{EP}}\approx 0.74T_{c}<T_{\mathrm{max}}$. Particularly, for values of $|\gamma_{1}^{0}-\gamma_{2}^{0}|/2$ greater than $|t-\Delta|_{\mathrm{max}}=\frac{t_{L}t_{R}}{4\Delta_{0}}$ ($\mu_{\mathrm{ABS}}=\Delta_{0}$), the system would never leave the non-trivial topological phase. By rewriting the couplings as $t_{L}=\eta_{L}\Delta_{0}$ and $t_{R}=\eta_{R}\Delta_{0}$, this condition reads $|\gamma_{1}^{0}-\gamma_{2}^{0}|>\eta_{L}\eta_{R}\Delta_{0}/2$. Since in realistic settings both $\eta_{L},\eta_{R}\ll 1$, this implies that there is no need of a huge coupling asymmetry for having MZMs.

We now calculate the entropy of a minimal Kitaev chain (Fig. 4d) including the temperature dependence of the effective parameters $\Delta$ and $t$, and hence of the eigenvalues that bifurcate. By tuning $\mu_{\mathrm{ABS}}=\Delta_{0}$, the system approaches the sweet-spot regime for temperatures below $T_{\mathrm{EP}}$, where the poles $\epsilon_{\mathrm{outer}}^{\pm}$ exhibit an EP (Fig. 4b). On the other hand, as explained in the main text, the entropy jumps associated with these poles occur at the crossover temperatures $T_{\mathrm{outer}}^{\pm}=|\epsilon_{\mathrm{outer}}^{\pm}|/2$. If the magnitudes of the poles differ, $|\epsilon_{\mathrm{outer}}^{-}|\neq|\epsilon_{\mathrm{outer}}^{+}|$, a fractional plateau $S=\log(2)/2$ appears, and its width is given by $T_{\mathrm{outer}}^{-}-T_{\mathrm{outer}}^{+}$. Thus, this fractional plateau is observable for temperatures below $T_{\mathrm{outer}}^{-}$. Importantly, this temperature is not directly given by $T_{\mathrm{EP}}$: although the system hosts MZMs up to $T_{\mathrm{EP}}$, the fractional plateau only emerges when the mode most coupled to the reservoir has not yet contributed to the entropy, which happens at $T_{\mathrm{outer}}^{-}$. Here, it is important to emphasize that since $\epsilon_{\mathrm{outer}}^{-}(T)$ depends on the temperature (Fig. 4b), the value of this crossover temperature becomes a self-consistent problem,

Thus, $T_{\mathrm{outer}}^{-}$ will be given by the crossing point between the curve $|\epsilon_{\mathrm{outer}}^{-}(T)|/2$ and the line $y(T)=T$ (and the same for the rest of poles), see Fig. 4c. Importantly, in this figure all the crossover temperatures lie in the lilac region, where the limit $T\to 0$ is valid. Under this limit, we have in general $T_{\mathrm{outer}}^{-}=\max(\gamma_{1}^{0},\gamma_{2}^{0})/2$, being completely valid as long as $\gamma_{1}^{0},\gamma_{2}^{0}<2T^{*}$.

Specifically, for the particular choice of parameters in Fig. 4 (cyan line of panel d for $\mu_{\mathrm{ABS}}=\Delta_{0}$), this temperature is one order of magnitude smaller than $T_{c}$ ($T_{\mathrm{outer}}^{-}=0.088T_{c}<T^{*}$), so that the $T\to 0$ approximation is valid for this and the rest of the poles, Fig. 4c. Using the Aluminium $T_{c}$ this gives a temperature $T_{\mathrm{outer}}^{-}\approx 15\mu\mathrm{eV}$ (which in real temperature units corresponds to a temperature of $T_{\mathrm{outer}}^{-}\approx 170mK$ well within the experimental range of observability) where the fractional plateau (cyan line for $\mu_{\mathrm{ABS}}=\Delta_{0}$) develops. Importantly, in the experimental demonstrations of minimal Kitaev chains based on quantum dots  Dvir et al. (2023); ten Haaf et al. (2024a, b) $\Delta\approx 10-20\mu$eV, in the same energy range of our estimation of the temperature range below which the fractional plateau in the entropy emerges. This supports our claim that unreasonably low temperatures $T\ll\Delta$ are not needed to observe our predicted fractional plateaus.

Recently it has been demonstrated that one can measure entropies of mesoscopic systems, either via Maxwell relations Hartman et al. (2018); Child et al. (2022) or via thermopower Kleeorin et al. (2019); Pyurbeeva et al. (2021). The Maxwell relation method relies on continuously changing a parameter $x$ (e.g. chemical potential or magnetic field) while measuring its conjugate variable $y$ (e.g. electron number or magnetization, respectively) such that $y=\partial F/\partial x$. Then, the Maxwell relation yields $\partial S/\partial x=-\partial y/\partial T$. Here we propose a novel application of this procedure, employing the Josephson current $I(\phi)$, which gives

Since the phase difference on the Josephson junction can be controlled (Fig. 2a) by e.g. embedding it in a SQUID loop Bargerbos et al. (2023); Pita-Vidal et al. (2023), one can integrate $dI/dT$ between $\phi_{1}=0$ and $\phi_{2}=\pi$. From Fig. 3b we expect $\Delta S$ to change from zero at high-T to $\log 2/2$ at low-T, an unequivocal signature of MZMs in the junction.

Our estimates using a microscopic model of a minimal Kitaev chain based on two QDs coupled through an Andreev bound state localized in a hybrid semiconductor-superconductor segment Liu et al. (2022), and including the temperature dependence of the effective parameters for cross Andreev reflection ($\Delta$) and single-electron elastic cotunneling ($t$), give an estimated temperature to observe the fractional plateau of $T\approx 170mK$, well within the experimental range of observability (see previous Section).

Note added While finishing this manuscript, two recent preprints in Arxiv Shen et al. (2024) and Beenakker (2024) also pointed out the subtleties of calculating the free energy in a NH Josephson junction. Eq. (9) in Ref. Shen et al. (2024) for the supercurrent agrees with our Eq. (7) in the limit $D\to\infty$. Moreover, Eq. (16) in Ref. Beenakker (2024) agrees with our Eq. (7) in the regime without EPs. This latter case, in particular, results in a reduction factor $\frac{2}{\pi}\arctan\left(\frac{E}{\gamma}\right)$ in the $T\to 0$ supercurrent, see Eq. (8), owing to the coupling with the reservoir.