Thermodynamics of a Phase-Driven Proximity Josephson Junction

We study the thermodynamic properties of a superconductor/normal metal/superconductor Josephson junction in the short limit. Owing to the proximity effect, such a junction constitutes a thermodynamic system where phase difference, supercurrent, temperature and entropy are thermodynamical variables connected by equations of state. These allow conceiving quasi-static processes that we characterize in terms of heat and work exchanged. Finally, we combine such processes to construct a Josephson-based Otto and Stirling cycles. We study the related performance in both engine and refrigerator operating mode.


Introduction
Thermodynamic concepts have been recently considered at the nanoscale, conceiving and realizing systems where quantum coherent properties are mirrored in thermodynamic quantities at mesoscopic level [1][2][3][4][5][6][7][8][9][10][11][12][13][14] . Furthermore, one of the most impressive examples of quantum features reflected in macroscopic systems is represented by superconductivity, where quantum coherence is manifested at a mesoscopic scale. Therefore superconducting systems are interesting platform where investigating the interplay between thermodynamic concepts and quantum coherences.
In this paper, we review the equilibrium thermodynamic properties of a hybrid system based on a Superconductor/Normal metal/Superconductor (SNS) Josephson Junction in the diffusive limit.
The behavior of such a system is ruled by the proximity effect, which consists in a set of physical phenomena owing to the propagation of the superconducting electron correlations in the normal metal [68][69][70]. In particular, guided by a matter of thermodynamic consistency, we discuss a relation between the electronic and thermal properties of the proximized system. From this relation, we develop a basic investigation of the thermodynamic properties of such a system. These results are then exploited to investigate quasi-static processes and thermodynamic cycles. We focus within a semi-classical regime of Josephson coupling, i.e., we neglect non-commutativity between the phase and the number of pairs, as usually done in the thermodynamic limit.
We remark that, besides the system studied in this paper, many equilibrium thermodynamic properties have been investigated in different conditions, theoretically and experimentally: thermodynamics of rings interrupted by insulating Josephson Junction [71][72][73][74][75], heat capacity in SN systems [76][77][78], free energy in hybrid SN systems due to boundary effects with approaches different to the quasi-classical theory [79][80][81].
The paper is organized as follows. Section 2 describes the proximized system under study and introduces its thermodynamics, giving also an insight into the underlying microscopical mechanism. Section 3 studies the thermodynamic processes. Hence, these are combined in Section 4 to investigate two different thermodynamic cycles. Finally, Section 6 summarizes and discusses the main findings. For completeness, Appendix A discusses the thermodynamics of a Josephson junction close to the critical temperature.

Model
We consider a system as sketched in Figure 1, constituted by a superconducting ring interrupted by a Superconductor/Normal metal/Superconductor (SNS) proximity Josephson Junction. The superconductor gap depending on temperature T is ∆(T) and reaches ∆ 0 at T = 0. The critical temperature is T c . The phase difference ϕ of the superconducting order parameter across the junction is ruled by the magnetic flux threading the ring, owing to the fluxoid quantization relation ϕ = 2πΦ/Φ 0 , where Φ 0 = h/2e ≈ 2 × 10 −15 Wb is the flux quantum. (a) Sketch of the SNS proximized system. It consists of superconducting ring, L S long, pierced by a magnetic flux Φ. The ring is interrupted by a normal metal weak link. The electron system of the whole device is thermally and electrically isolated and at temperature T. The system is connected to a thermal reservoir at temperature T through a heat valve v. (b) Magnification of the SNS junction. The normal metal weak, L N long, is in clean electric contact with the superconducting leads. A j , N j are respectively the cross-section and the DoS at Fermi energy of the j = N or S metal. The phase drop ϕ of the superconducting order parameter takes place across the junction.
We assume that the system is thermally and electrically isolated and at a homogeneous temperature T, neglecting thermal gradients. We consider only heat exchange with a reservoir at temperature T through the respective heat channel connected by a heat valve v [33][34][35][36][37][38][39][40][41][42], as drawn in Figure 1.
The junction, magnified in Figure 1b, consists in the two S leads in electric contact with an N weak link. The superconductor has a critical temperature T c and BCS gap at zero temperature ∆ 0 . The N weak link and the S leads have respectively cross-sections A N and A S , conductivities σ N and σ S , Density of States (DoS) per spin at the Fermi level N N and N S . The length of the weak link is L N , resulting in a resistance R N = L N /A N σ N . The length of the superconducting ring is L S . The whole proximized system (ring+junction) volume is V.
We make the following assumptions about the junction, in order to make simple analytical predictions within the Kulik-Omel'yanchuk (KO) treatment [82][83][84]. We consider diffusive charge transport with diffusivity D for both the S and N parts. This requires that the weak link is longer than the mean free path L mfp : L N L mfp . The KO treatment holds when the whole junction can be treated in a quasi-1 dimensional approximation, i.e., when A S , A N ξ 2 . The diffusivity defines the coherence length ξ = √h D/∆ 0 [85]. Moreover, we consider a short constriction weak link respect to the superconducting leads. Quantitatively, using the parameters l = L N /ξ and a = σ S A S /σ N A N , we consider a short junction with l 1 and a constriction with al 1. The typical values for this kind of system are the following. The ring can be made of aluminium, with ∆ 0 ≈ 180 eV, corresponding to T c ≈ 1.2K [86][87][88]. The coherence length for hybrid Al-based devices is about ξ ≈ 150 nm [89][90][91]. In the following, we set the Boltzmann constant to k B = 1, implying that the temperatures have a physical dimension of energy while entropy and specific heat are dimensionless.

Hybrid Junction as Thermodynamic System
Before investigating the thermodynamic behavior of our system in detail, we discuss about the thermodynamic consistency under a general point of view that is valid for any Josephson Junction (JJ). In particular, we focus on the relation between the current transport and the junction entropy.
In a JJ, the Current Phase Relation (CPR) describes the dissipationless supercurrent I(ϕ, T) flowing across it as function of the phase difference ϕ and temperature T [70,83]. The precise form of the CPR depends on the geometry and on the materials of the junction, and can be calculated from the free energy ash where F(ϕ, T) has to be calculated within quantum statistical methods as a function of the state variables (ϕ, T). The CPR constitutes an equation of state connecting I, ϕ and T. Another equation of state is given by the entropy S(ϕ, T) as a function of phase difference and temperature The entropy and the CPR are necessarily linked by thermodynamic consistency. Indeed the two cross derivatives of F are identical, i.e., ∂ ϕ ∂ T F = ∂ T ∂ ϕ F, owing to the Schwarz theorem. Hence, the following Maxwell relation is universally valid Using this equation, the entropy of the JJ can be expressed as where S 0 (T) is the entropy at ϕ = 0 and δS(ϕ, T) is the phase-dependent entropy variation where E (ϕ, T) is the Josephson energy stored in the junction at a given temperature T, R 0 = h/2e 2 ≈ 12.9 kΩ is the inverse of the quantum of conductance. We note that the prefactor in Equation (6) is usually expressed as Φ 0 /2π. We chose the form eR 0 /2π to allow an easier comparison with the junction resistance R N .
The entropy S 0 (T) at ϕ = 0 cannot be determined from the knowledge of the CPR. Indeed, any function S 0 (T) of the temperature is dropped by the phase derivative in Equation (3), hence satisfying the Maxwell equation. The physical solution of S 0 (T) can be found within a microscopic model that we show in the next subsection.

Proximity Induced Minigap
In this subsection, we give an insight into the microscopic mechanism which determines the entropy in a hybrid junction. In particular, we show that the entropy dependence on temperature and phase is related to the presence of an induced phase-dependent minigap in the quasi-particle Density of States (DoS). In a hybrid NS, correlated electrons propagate from the superconductor into the normal metal, strongly modifying the properties of the latter with a set of phenomena called generically under the name of proximity effect [68,69,[92][93][94][95][96][97]. Among all possible consequences dictated by the proximity effect, here we focus on the induced mini-gap in the quasi-particle Density of States, it being responsible for the phase and temperature dependence of the entropy S in an SNS junction.
Let us consider the N weak link in an SNS junction. When not proximized, the weak link DoS is homogeneous and approximately constant at its Fermi level value N N in the energy range of interest of few ∆ 0 around the Fermi energy. Instead, when proximized by the superconducting leads, the DoS is no more constant neither on energy nor on position, but is given by N N N(r, ε, ϕ), where N(r, ε, ϕ) is the normalized local DoS [85,94,98,99] that is a function of the position r, energy ε and the phase difference ϕ.
One way to calculate the normalized local DoS is provided by the quasi-classical theory of superconductivity [85,[94][95][96][98][99][100]. Qualitatively, a result of this theory is that the normalized local DoS is characterized by an induced gap in the N weak link, whose amplitude∆ is smaller than the S bulk gap ∆(T). For this reason,∆ is dubbed induced minigap. This induced minigap has the following properties [85,98]: its width∆ at ϕ = 0 depends on the weak link length L N and reaches∆ → ∆(T) when L N is well below the coherence length ξ. Moreover,∆ depends on the phase ϕ through a function that is even and 2π periodic. The minigap is fully open at ϕ = 0 and shrinks till closure at ϕ = π. An analytical solution of the local normalized DoS is available for diffusive short junctions with rigid boundary conditions [101,102], yielding that∆ = ∆(T)| cos(ϕ/2)|. The proximity induced gap and its interesting properties have been observed experimentally by tunneling experiments [89,90,103].
An important feature of this microscopic proximity DoS modification is that it does not take place just in the N weak link, but also affects the S leads as well. The anti-proximization operated by the N weak link on the S leads is called inverse proximity effect and plays the role of a crucial correction in short junctions, since it gives an important contribution to the total entropy dependence on the junction phase [84].
A numerical example of the phase-dependence of the local normalized DoS N in a junction is reported in Figure 2, within the quasi-classical methods of Reference [84], calculated for a junction with parameters l = 0.1 and a = 10, ∆(T → 0) = ∆ 0 . The color plots show the evolution of the normalized local DoS N versus energy ε and spatial position x for four values of ϕ from 0 to π. The blue area corresponds to the gapped part of the local DoS; the white area is the saturation color that is associated to the divergence of the DoS at the gap edges. The position is normalized to the coherence length ξ: as shown in the first panel, the central zone x ∈ [0, 1] coincides with the N weak link, while the lateral zones are the superconducting leads. At ϕ = 0, the DoS is homogeneous and is approximatively given by the BCS form The spatial homogeneity is due to the fact that the calculation involves a short junction, otherwise the induced minigap would have been smaller than ∆ 0 [84,101,102]. Increasing ϕ, the induced minigap shrinks till the complete closure at ϕ = π. It is possible to appreciate also the inverse proximity effect in the S leads, outside the stripe delimited by the red dashed lines in Figure 2. Color plots of the quasi-particle local normalized Density of States (DoS) N in a Superconductor/Normal metal/Superconductor (SNS) junction, versus energy ε and position x, for ϕ = 0, π/3, 2π/3, π. The dashed lines separate the S regions (on the sides) to the N region (in the center), as shown by the junction sketch. The phase dependence of the DoS is mirrored in a phase-dependence of the junction entropy S. The numerical calculation has been obtained within the quasi-classical methods of Reference [84] with a = 10, The phase dependence of the quasi-particle DoS implies a phase dependence of the junction entropy. The total entropy is [84,[104][105][106][107] S(ϕ, T) = Vol S(ϕ, T, r)dV (8) where N r is N S or N N whether r is in the leads or the weak link. At this point, we have two ways to calculate the entropy variation δS(ϕ, T). One consists in calculating δS(ϕ, T) from I(ϕ, T) exploiting the Maxwell relation through Equations (3)- (6). The other way is by means of Equations (8) and (9) given by the statistical argument above concerning the quasi-particle density of states. It is a non-trivial result that the two approaches give results in agreement [80,84,108]. This is an equilibrium thermodynamic feature due to the fact that the equilibrium supercurrent is carried by the Andreev Bound States (ABS), whose spectral density is non zero below the superconducting gap |ε| < ∆(T) [85,101]. The quasi-particle DoS and the ABS spectral density are both functions of ϕ, ensuring that the two approaches are equivalent.
We conclude this discussion by calculating S 0 . As discussed in Section 2.2, this quantity can not be obtained by the Maxwell relation (3), constituting hence an undetermined function of the temperature T in Equation (4). However, S 0 can be determined with a statistical mechanics approach. Given the assumptions of Section 2.1 of short junction l 1 and a 1, the local normalized DoS at ϕ = 0 is given by the BCS expression (7) [84]. Hence In obtaining this expression from (9), we neglected that N S = N N in general. However, since the junction volume is negligible respect to the total volume, we have approximated the prefactor with N S V leads + N N V weak L ≈ VN S . Below, we discuss this result within the full dependence of the total entropy S on ϕ and T.

Kulik-Omel'yanchuk Theory
The Kulik-Omel'yanchuk theory, whose assumptions have been introduced in Section 2.1, provides the CPR [82,83,101] This expression [101] is equivalent to the Matsubara summation form presented in the first paper about the KO CPR [82]. Here we adopt the integral form that allows to find simple closed expressions in the limit T ∆ 0 . In the zero-temperature limit T → 0, the KO CPR reduces to [82] I(ϕ, We use as scale for the supercurrent the critical current I c at T = 0, obtained by maximizing (12). Numerical maximization returns that I c is where κ ≈ 1.33. The maximum is placed at phase ϕ ≈ 1.97 ≈ 0.63π. The KO CPR is shown in Figure 3a, normalized to I c . The T = 0 curve in Equation (12) is plotted in black dotted. As one can see the supercurrent decreases versus temperature, passing from a skewed shape to a more sinusoidal shape [83]. According to the prescription given in (6), the associated Josephson energy to the KO CPR is (14) where R 0 = h/2e 2 ≈ 12.9 kΩ is the inverse of the conductance quantum. The characteristics of E (ϕ, T) are plotted in Figure 3b. Being the integral of the supercurrent, the Josephson energy increases versus temperature. At T = 0, E reduces to The maximum Josephson energy is E 0 = E (ϕ = π, T = 0), given by From E it is possible to calculate δS. Figure 3c reports the entropy variation δS(ϕ, T) calculated numerically with δS(ϕ, T) = −∂ T E (ϕ, T), for chosen temperatures in legend. It can be noticed that δS decreases with the temperature, consistently with the third law of thermodynamics.
At low temperatures, where ∂ T ∆(T) → 0, a closed form of δS can be obtained by the temperature derivative of Equation (14), yielding [26,84,101] The behavior of the entropy can be qualitatively grasped with the minigap mechanism. Let us consider a fixed temperature T. Hence, the distribution function f log f in Equation (9) has a certain bandwidth of the order T. At low temperature T ∆ 0 and ϕ = 0, the DoS gap has width ∆ 0 and the distribution bandwidth is smaller than the gap. Hence, the lack of available states exponentially suppresses the entropy. When ϕ moves from ϕ = 0 to ϕ = π, the minigap shrinks giving new available states for the distribution f log f , increasing the entropy. At T ∆ 0 and short junction, it is approximately∆ = ∆ 0 | cos(ϕ/2)| [102], the matching phase between the minigap and the distribution bandwidth is 2 arccos(T/∆ 0 ), at which the entropy increases. This is particularly evident in the curve T = 0.1T c in Figure 3c, where δS is negligible except close to ϕ → π.

Total Entropy
Given the microscopic and the KO CPR considerations of the last subsections, we can study the total entropy, that is S(ϕ, T) = S 0 (T) + δS(ϕ, T) where S 0 is given by the BCS entropy in Equation (10) and δS = −∂ T E where E is given by expression (14). We note that the first term scales as ∆ 0 N S V, while the second as eR 0 I c /∆ 0 . For this reason, it is convenient to introduce a parameter α of the system that sets the ratio between these two quantities: α characterizes the relative influence of the phase-dependent term δS over the remaining term S 0 . The quantity α can be experimentally determined by heat capacity measurements, as explained in Section 3.2. Moreover, α controls the temperature of a first-order transition to the normal state when ϕ = 0, discussed in detail in Appendix A. Figure 4 reports the total entropy for α = 0.6. Different values ϕ in the legend are plotted, showing the increase of S from ϕ = 0 to ϕ = π. The four curves correspond to the DoS states in the frames of Figure 2. As expected, the closure of the minigap from ϕ = 0 to ϕ = π implies an increase of entropy. The scale of this increase is set by α. In the following, the calculations are obtained with α = 0.6. This value evidences the entropy variation and the related results while keeping a proximized volume negligible respect to the total volume, as shown below in this subsection, and keeping the unwanted first-order transition above the temperature 0.7T c , as discussed in Appendix A. Considering that N S ≈ 7 × 10 46 m −3 J −1 [109], α = 0.6 corresponds to a ratio I/V ≈ 20 mAµm −3 .
The behavior of the entropy can be studied in more detail at low temperature T ∆ 0 , where closed expressions can be obtained. At ϕ = 0, the DoS has the BCS form in the whole volume of the device, returning the exponentially suppressed behavior of entropy described by the red curve in Figure 4. Hence, at low temperatures T ∆ 0 and ϕ = 0 the entropy can be approximated by the expression [110,111] At ϕ = π, the minigap is closed and a proximized spatial region around the weak link has a metallic-like DoS. The entropy density (9) is then exponentially suppressed in the leads and with a linear-in-temperature dependence in the proximized region. This is confirmed by an analytical expression for δS that can be obtained at low temperatures at ϕ = π. Substituting ϕ = π in (17) we obtain Developing the logarithm around ε = 0 as log(1 − 2ε/(∆ 0 + ε)) ≈ −2ε/∆ 0 and substituting ε/T = z, For T → 0, we obtain The linear behavior of δS(ϕ = π, T → 0) allows to neglect the exponentially suppressed S 0 contribution to the total entropy S, allowing the following approximation Figure 4b reports the low-temperature behavior of the entropy for 0 ≤ T ≤ 0.4T c . The ϕ = 0 and ϕ = π curves show the exponentially suppressed and linear behavior respectively. The purple dashed curve report the analytical expression (23), revealing a good agreement at T < 0.2T c .
A qualitative explanation of the drastic change of the entropy behavior versus phase difference can be done within the mechanism presented in Section 2.3. When ϕ = 0, the local normalized DoS N(ϕ, ε, r) is homogeneous over space and gapped according to the BCS expression. As a consequence, the entropy density S in (9) is exponentially suppressed and independent on the position r. When 0 < ϕ < π, the DoS is altered: this alteration can be roughly described as an effective proximized volumeṼ, where the DoS is phase-dependent with minigap∆(ϕ), while in the non-proximized rest of the system the DoS is unchanged with the BCS gap ∆ 0 . As a consequence, the entropy contribution from the proximized region dominates over the entropy contribution from the non-proximized region, since from the proximized region it is S ∝ (∆(ϕ)/T) 1/2 e −∆(ϕ)/T while from the non-proximized region it is S ∝ (∆ 0 /T) 1/2 e −1/∆ 0 . This behavior can be noticed in Figure 4b, comparing the curves at ϕ = π/3 and ϕ = 2π/3 with the one at ϕ = 0. The curves at ϕ = 0 show a suppressed region in a low temperature interval whose width depends on∆(ϕ). Finally, when ϕ = π, the induced minigap is closed and the behavior is radically changed from exponentially suppressed to linear.
We conclude this Section with some remarks about α. In our treatment, α is a free parameter to be set to get numerical results. However, there is a physical upper limit to its value. Maximizing α can be done experimentally by maximizing I c and minimizing V. However, this can not be done at will, since our approach is based on the KO theory, that requires as assumption that the two leads are good reservoirs of electron coherence, i.e., the inverse proximity effect by the weak link does not spoil the bulk superconducting properties of the leads. We show the existence of this upper limit with the following two arguments.
The first is given by expressing α in terms of the system geometrical properties. Taking into account expression (13) for I c , the coherence length ξ 2 =hD/∆ 0 , the S conductivity σ S = 2e 2 N S D and that the volume is V ≈ A S L S , we have The requirement al 1, given in Section 2.1, implies hence that the only free parameter for increasing α is to decrease the length of the ring as most as allowed by the practical geometrical realization.
The second argument about the physical upper limit of α is that the volume V of the system must be in any case bigger than the effective proximized volumeṼ, involving both the weak link and the in inverse proximized leads. This can be obtained in a qualitatively by considering that when ϕ = 0, the minigap is closed and the DoS is modified in a region surrounding the weak link in such a way to return expression (23), that is linear like a normal metal. The volumeṼ of the proximized region can be estimated by comparing expression (24) with the entropy of a normal metal S N = 2π 2 N SṼ T/3: This expression suggests that the inverse proximized region is present in the leads for a length ξ/al. The volume of the proximized region does not coincide with the weak link region. In particular, they scale differently on the junction length L N , since the volume of the weak link is ∝ L N A N , while the proximized volume isṼ ∝ I c ∝ A N /L N . This point shows that the proximized region is not confined in the weak link but is extended also in the leads owing to the inverse proximity effect [84].
By imposing that the proximized volume is smaller than the system volume,Ṽ V, Equation (26) yields In our calculations, α = 0.6, corresponding to a ratioṼ/V ≈ 0.15. Finally, another argument that estimates the physical upper limit of α concerns the fact that a high ratio of I c /V decreases the critical temperature of the system. This point is discussed in details in Appendix A.

Thermodynamic Processes
In this section we discuss thermodynamic processes focussing on quasi-static situation, meaning that the device passes through a succession of equilibrium states. This condition can be met if one considers a sufficiently slow speed of the process under inspection. This speed is set by the leading (fastest) thermalization mechanism, that is the electron-electron (e-e) interaction.
Here, we study three different thermodynamic processes. The first is an isothermal one, where the phase is changed while the temperature is kept constant. The second is the isophasic where the temperature is changed while the phase is kept constant. Finally, we consider the isentropic process where the phase is changed while the system exchanges no heat with the universe, thus retaining entropy constant during the process. We will give particular attention to processes with phase variation only between ϕ = 0 and ϕ = π, for two reasons. First of all, for these two values of the phase difference, the circulating supercurrent is null and we can neglect any inductive contribution from the ring to the total energy when investigating thermodynamic cycles. Secondly, these particular values of ϕ admit simple and closed expressions, allowing for a simple and analytical discussion within KO theory.
For a quasi-static thermodynamic process, the heat flow between the initial and the final state can be written as is Q = P TdS (28) and the work released is where the integrals are meant to be line integrals over the path P in the space of the thermodynamic variables. For quasi-static processes, P lies in the surface of the equilibrium states. We will consider the three different paths corresponding to the isothermal, isophasic and isentropic situation. Hereafter, heat and work integrals are defined according to the following sign convention: the work W is positive when the system releases work to the universe, while the heat Q is positive when the system absorbs heat from the universe. According to this convention, the energy conservation over a closed loop path reads Q − W = 0.

Isothermal Process
Let us consider a isothermal process from an initial state i at (ϕ i , T) to a final state f at (ϕ f , T). This can be realized by keeping open the heat valve toward the reservoir sketched in Figure 1a. For notation simplicity, here we indicate both the system temperature and the reservoir temperature as T, implying that at thermal equilibrium T = T, where T is the reservoir temperature.
In this case the work released by the system is where E represents the Josephson energy in the junction defined in Equation (5). For a process where the universe has to perform a work on the system, the sign of W i f is negative, consistently with the convention of W defined above. The heat absorbed during this process is Heat is absorbed when ϕ goes from 0 to π, owing to the closure of the minigap. It is worth to note that in the isothermal process we do not explicitly rely on the BCS contribution S 0 (T). Hence, the thermodynamic consistency requires that an isothermal process must exchange heat. Interestingly, the supercurrent is not directly involved in this heat exchange, since the Cooper pair system carries no entropy and the supercurrent flow is dissipationless. Instead, the heat is absorbed by the quasi-particle (excited states of system) from an external system, i.e., in our scheme from the external reservoir at fixed temperature T. If heat is not supplied, the system undergoes an adiabatic transformation, treated in the next subsection. Below in Section 5 we discuss some strategies to measure this heat exchange. At low temperature T ∆ 0 , the heat absorbed and the work released in an isothermal process from ϕ = 0 to ϕ = π can be calculated exploiting the expression (23).
where the second equivalence is due to the fact that the temperature is constant during an isothermal process.
The released work at low temperature is obtained by calculating the expression of E at low temperature. From Equation (23) where we have used the expression (16) for E 0 . The work at low temperature for a ϕ = 0 → π isothermal is As expected, the work released scales as the critical supercurrent I c and increases in module by lowering the temperature.

Isophasic Process and Heat Capacity
In an isophasic process, the phase difference ϕ is kept constant while the temperature is changed. Considering Figure 1, this can be done by opening the thermal valve toward the reservoir while the threading flux Φ is fixed. The system passes from its initial temperature T i to the final temperature The work exchanged is then null, since dW = −eR 0 Idϕ/2π. The system exchanges energy only in the form of heat. In a process from (ϕ, T i ) to (ϕ, T f ), the exchanged heat can be written as At low temperature T i , T f ∆ 0 , using Equations (20) and (23), it is possible to obtain the isophasic heat in a closed form for ϕ = 0 and ϕ = π. At ϕ = 0 we get while for ϕ = π we have The heat exchanged in an isophasic process brings naturally to the concept of heat capacity. Indeed the heat exchanged can be expressed as a function of the initial and final temperatures as where C(ϕ, T) is the isophasic heat capacity: The importance of the heat capacity relies also on the fact that is an experimental observable quantity. Indeed, by definition, can be measured as the temperature response of the system to a heat pulse. From Equations (36) and (37) it is evident that the amount of heat exchanged for an isophasic process from T i to T f depends on the phase ϕ and, hence, the heat capacity is dependent on ϕ. From these expressions, it is possible to obtain the isophasic heat capacity at low temperature T ∆ 0 in a closed form. At ϕ = 0 Similarly to the entropy, the heat capacity assumes two different behaviour,passing from a suppressed superconducting-like to a linear metallic-like behavior whether the phase is ϕ = 0 or ϕ = π, respectively. Moreover, the ratio of (41) over (40) scales like α. Hence, an experimental measurement of C at ϕ = 0 and ϕ = π can be used to get an estimated value for the dimensionless parameter α discussed before.
Here, a sort of parallelism between the variables (I, ϕ) of our system and (p, V) of an ideal gas can be noticed. In the same analogy, the work differential eR 0 Idϕ/2π plays the role of the pdV differential for a classic gas. However, the phase difference ϕ variable is 2π periodic, differently from V.
In a generic situation, the heat capacity can be evaluated numerically from Equation (39). Figure 5 reports C(ϕ, T) for α = 0.6. Figure 5a is a color plot of C versus ϕ and T; panels (b) and (c) are cuts of the color plot versus T and ϕ. Looking at panels (a) and (b) one can note that C(ϕ, T) goes from a gapped-like behavior at ϕ = 0 to a linear behavior at ϕ = π. This is confirmed also by the dashed green line plotting the analyitical expression (41). Since C is the temperature derivative of S, its behavior can be explained qualitatively within the phase-dependent minigap mechanism, in the same fashion provided for the entropy in Section 2.

Isentropic Process
In a isentropic process the entropy of the system is conserved. Since we are considering quasi-static processes, an isentropic is also adiabatic in the thermodynamic meaning that no heat is exchanged with the universe. Indeed, for a quasi-static process, it holds dQ = TdS. The isentropic process can be physically realized when the system is thermally isolated, i.e., when the heat valve in Figure 1a is closed.
In a isentropic process, the constraint of constant S implies an implicit relation between the phase and the temperature. Let us consider an isentropic process that starts at the initial state (ϕ i , T i ). During the whole process, ϕ and T are related by the implicit equation For an isolated system, the phase and entropy (ϕ, S) are the independent variables that will specify the system state. The temperature is then a function T(ϕ, S). In detail, T(ϕ, S) decreases by increasing ϕ in the interval 0 < ϕ < π for fixed S. Indeed, since S(ϕ, T) is an increasing function in 0 < ϕ < π for fixed T and, as a consequence, the temperature of the system must decrease in order to keep S constant.
In particular, we focus on the isentropic temperature decrease for processes that start at phase ϕ = 0 and T i , where the initial state sets the entropy S(ϕ = 0, T i ) of the process. For these processes, we define the temperature decrease T f implicitly defined in Equation (42) as In Figure 6a,b we report the quantity T f (ϕ, T i )/T i , i.e., the relative temperature decrease, for a system with α = 0.6. T f /T i is enhanced toward low T i , since for T ∆ 0 the behaviors of S(ϕ = 0, T) and S(ϕ = π, T) are strongly different: the former is exponentially suppressed while the latter is linear (see Figure 4a,b).
A closed expression for T f /T i can be obtained for T ∆ 0 by exploiting equations (20), (23): The isentropic cooling is reminiscent of the adiabatic cooling process typical for the expansion of an ideal gas. The analogy goes forward when discussing in terms of available states. Indeed, when the gap reduces to closure in the process ϕ = 0 → ϕ = π the number of available states increases so that the temperature decreases to keep the entropy constant. The same thing happens in the case of an adiabatic expansion of a gas, where the position states are increased by the volume increase. Figure 6c plots T f /T i versus T i for different values of α. The relative cooling is more effective for higher α, since higher values of α correspond to a stronger weight of the proximized region, where the gap can be tuned, over the phase independent superconducting leads. In particular, T f /T i ∝ α −1 as shown in Equation (44). However, the passage from a gapped to a gapless state yields a strong temperature cooling even for moderate values of α, provided that T i is low enough. It is interesting to investigate how the CPR of a junction is modified by the assumption that during the change of the phase difference ϕ the entropy remains constant. For sake of simplicity, here we focus on isentropic CPR where the entropy is set by the initial state by S(ϕ i = 0, T i ). In such case the Josephson current can be calculated substituting the temperature T f (ϕ, T i ) in the isothermal CPR (11): In Figure 6d we show the comparison between the isothermal and isentropic CPR for α = 0.6 for the initial temperature T i = 0.6T c . It is worthy to notice that the isothermal CPR depends only on the nature of the junction, while for the isentropic case there is also a dependence on α, which includes the total volume of the system. The dashed lines are isothermal curves at the initial temperature T i = 0.6T c and the final temperature T f (T i , ϕ = π) = 0.51T c . We can notice that the two isothermal respectively overlap the isentropic at ϕ → 0 and ϕ → π. Moreover, the fact that the isentropic curve lies between the two isothermal curves indicates that the temperature is between the initial and final temperatures, since T evolves from T i to T f during the isentropic process. The work for an isentropic ϕ = 0 → π with initial temperature T i is given by Since the isentropic CPR is constrained between the isothermal CPRs at T i and T f , i.e., I(ϕ, T i ) < I S (ϕ, T i ) < I(ϕ, T f ) as shown in Figure 6d, the isentropic work is equally constrained between the isothermal works at T i and T f .

Thermodynamic Cycles
The combination of different thermodynamic processes, studied in the previous section, allows constructing thermodynamic cycles. In this section, we present two possible examples of thermodynamic cycles that can be built based on the various processes discussed above. In particular, we focus on two cycles that we call Josephson-Otto cycle and Josephson-Stirling cycle, thanks to their analogy with classic thermodynamic counterpart. We first explain their implementation and then we discuss their performances.
To this aim, we consider the hybrid system attached to two different reservoirs, identified as Left Reservoir (L) and right Reservoir (R), in the sketch of Figure 7. The two reservoirs are at fixed temperature T j and can release heat Q j to the system through a heat channel controlled by a heat valve v j , where the subscript j can be L or R, respectively. We consider Q j positive when the heat flows from the reservoir to the system, in agreement with the sign convention defined in Section 3. Thereafter, we will study cycle characteristics as a function of the temperatures (T L , T R ). In particular, we will show that there are regions of (T L , T R ) where the cycles can operate as engine or refrigerator. The reservoir roles depend on the operating mode: when a cycle operates as engine, the two reservoirs play the role of the Hot Reservoir (HR) and Cold Reservoirs (CR), at temperatures T hr > T cr respectively. In a cycle, the system absorbs an amount Q hr from the HR and releases |Q cr | < Q hr heat to the CR. In practical systems, the cold reservoir can be constituted by the ambient, i.e., the large substrate thermalized to the cryostat, while the hot reservoir can be a heated subsystem, like a large metallic pad heated by Joule effect. Conversely, when the cycle is considered as a refrigerator, the two reservoirs play the role of the Cooled Subsystem (CS) and Heat Sink (HS), at temperatures T cs < T hs respectively. In a cycle, the system absorbs an amount of Q cs from the CS and releases |Q hs | > Q cs to the HS. In practice, the CS is an isolated subsystem from which the heat is extracted, where the heat capacity is assumed to be large enough to consider the CS as a reservoir within one cycle. The CS can be constituted of a metallic pad that can be used as cooled substrate for nanodevices. In practical systems, the heat sink is typically constituted by the ambient, i.e., the substrate thermalized to the cryostat in our device.
The cycle performances are characterized by inspecting several figures of merit. In the case of the engine we investigate the work released per cycle W and its efficiency, defined as This quantity is physically limited by the Carnot efficiency In the following subsections, we show W and η versus both the temperatures T cr , T hr . We discuss in detail the dependence of W and η as a function of T hr for fixed T cr , since in real systems it is most likely possible to tune the HR temperature while the CR temperature T cr is fixed by the ambient.
In the refrigerator mode, the figures of merit we consider are the extracted heat Q cs from the CS per cycle, and the Coefficient of Performance (COP), defined as Like the efficiency, the COP is limited physically by the Carnot COP limit In the following subsections, we show Q cs and the COP versus both temperatures T cs , T hs . We discuss in detail the dependence of Q cs and the COP as a function of T cs for fixed T hs , since in real systems the HS temperature T hs is given by the ambient and can not be tuned, while T cs decreases from the ambient temperature in the refrigeration process. Notice that the work W and the heat extracted Q cs are quantities defined per cycle. Hence, at cycling frequency ν, the engine returns a PowerẆ = Wν and the refrigerator returns a Cooling Power CP = Q cs ν.

Josephson-Otto Cycle
Here we study the Josephson-Otto cycle, by starting with the engine mode for sake of simplicity. The Josephson-Otto engine is described by the scheme in Figure 8, where the panels a and b show respectively the processes in the (T, S) and (ϕ, I) planes. The cycle is constituted by two isentropic processes, i.e., 1 → 2 and 3 → 4, and two isophasic processes, i.e., 2 → 3 and 4 → 1, see Figure 8. We choose by convention that the state 1 and 3 are thermalized to the R and L reservoir, respectively. In this way, the R and L reservoirs play respectively the role of the HR and CR.
The cycle is given by the succession of the following processes: • Isentropic 1 → 2. All thermal valves are closed to make the system thermally isolated. The system is driven from the state (ϕ 1 = 0, T 1 = T R ) to (ϕ 2 = π, T 2 ), where T 2 = T f (ϕ = π, T 1 ). In this process the universe spends a work |W 12 | (W 12 < 0 according to the convention defined in Section 3). |W 12 | is represented by the green area in Figure 8b. No heat is exchanged, Q 12 = 0. • Isophasic 2 → 3. By opening the thermal valve v L , the system goes from the state (ϕ = π, T 2 ) to (ϕ = π, T 3 = T L ). The system releases heat |Q 23 | to the left reservoir (magenta area in Figure 8a). No work is performed, W 23 = 0. • Isentropic 3 → 4. All thermal valves are again closed to make the system thermally isolated. The system is driven from the state (ϕ 3 = π, T 3 = T L ) to (ϕ 4 = 0, T 4 ). By construction, if T 2 > T L then it is T 4 < T R . In this process the system returns a work W 34 (W 34 > 0 according to our convention), represented by the sum of the green and blue areas in Figure 8b. No heat is exchanged, Q 34 = 0. • Isophasic 4 → 1. By opening the thermal valve v R , the system goes from the state (ϕ = 0, T 4 ) to (ϕ = 0, T 1 = T R ). The system absorbs heat Q 41 from the reservoir at T R (magenta+pink area in Figure 8a). No work is performed, W 41 = 0.
The total work released per cycle is The heat Q hr absorbed from the HR (correspondent to R) is From the two schemes presented in Figure 8 it can be noticed that the cycle operates as an engine if T 2 > T 3 . This condition requires that T L < T f (ϕ = π, T R ), i.e., a temperature gap between the two reservoirs is required. When T L approaches T f (ϕ = π, T R ) the cycle tends to the degenerate case reported in Figure 9a, where the two adiabatic curves tend to superimpose. Also in the (ϕ, I) plane the two adiabatic curves tend to superimpose, meaning that the net work is W = 0 at T L = T f (ϕ = π, T R ). On the contrary, if T L > T f (ϕ = π, T R ), the cycle is inverted as in Figure 9b. In this case, the cycle works as a refrigerator and the work is W < 0, i.e., made by the universe on the system. Hence, the curve in the plane (T L , T R ) where W = 0 can be defined as the characteristic curve of the Otto cycle. It separates the regions where the cycle is in the engine or refrigerator mode and it is given by the equation Close to the characteristic curve, in the case shown in Figure 9a, it is evident that Q L , Q R tend to zero but their ratio tends to Q R /Q L → T R /T L . This property is exploited below to calculate the limits of η and COP close to the characteristic curve. Let us consider the refrigerator mode in the case T L > T f (T R ), represented in Figure 9b. In this case, the cycle is clockwise and operates as a refrigerator. The two reservoirs play a different role: the R reservoir represents the Heat Sink, while the L one represents the Cooled Subsystem. The case T L > T f (T R ) coincides with the following cycle • Isentropic 1 → 2. All thermal valves are closed to make the system thermally isolated. The system is driven from the state at the ambient temperature (ϕ 1 = 0, In this process, the universe spends a work |W 12 | (W 12 < 0 for of Section 3). No heat is exchanged, Q 12 = 0. • Isophasic 2 → 3. By opening the thermal valve v L , the system goes from the state (ϕ = π, T 2 ) to (ϕ = π, T 3 = T L ), removing the heat Q 23 from the CS (magenta area in Figure 9b). No work is performed, W 23 = 0. • Isentropic 3 → 4. All thermal valves are closed. The system is driven from the state (ϕ 3 = π, T 3 = T L ) to (ϕ 4 = 0, T 4 ). Now, T 4 > T R . In this process, the system returns a work W 34 . No heat is exchanged, Q 34 = 0. • Isophasic 4 → 1. By opening the thermal valve v R , the system goes from the state (ϕ = 0, T 4 ) to (ϕ = 0, T 1 = T R ). The system releases heat Q 41 to the reservoir at T R , since T 4 > T R , which correspond to the magenta+pink area in Figure 9b. The temperature T 4 plays an analogous role of the hot heat exchanger that is present in the refrigerators. No work is performed, W 41 = 0.
In the refrigerator mode, the work released is still given by W = W 12 + W 34 . The heat Q cs absorbed by the CS is Figure 10 is a summary of the work released W and the heat absorbed Q hr and Q cs . Panels a,b are color plots of these quantities versus (T L , T R ). The dashed red curve represents the characteristic curve defined in Equation (53), corresponding to W = 0. Above it, for T L < T f (ϕ = π, T R ), the cycle operates as engine, while below it (T L > T f (ϕ = π, T R )) the cycle works as refrigerator. The orange dot-dashed curve reports the thermal equilibrium T L = T R . We can notice that below this curve, i.e., for T R < T L , there is a region where work is spent by the universe to pump heat from the L reservoir (the hotter one) to the R reservoir (the colder one). Hence, work is spent to perform a process that can be performed spontaneously. We define this region as a cold pump, following the definition given in References [112,113]. Figure 10c reports the released work versus the HR temperature T R for different values of the CR temperature T L as reported in the legend. The curves reach the value zero corresponding to the characteristic curve plotted in Figure 10a. We observe that the general trend of the work is to increase with the temperature difference T R − T L between the two reservoirs. The order of magnitude of the work per cycle is ∼ 0.1eR 0 I c . Figure 10d reports the absorbed heat Q cs = Q L versus the CS temperature T L for different values of the HS temperature T R . The black curve reports the case of the heat absorbed Q cs at T L = T R . The curves with fixed T R are limited on the right at T L = T R , to not include the Cold Pump case, see color plots in Figure 10. The curves with fixed T R goes to zero in correspondence of the characteristic curve, defining the minimum achievable temperature of the refrigerator. The refrigerator can not physically cool below the minimum achievable temperature, since the absorbed heat reaches Q cs = 0. The black curve reporting Q cs at T L = T R is important since it reports the heat absorbed per cycle when the refrigerator starts to operate at the thermal equilibrium. Hence, for a cycling frequency ν, the corresponding cooling powerQ cs = Q cs ν for T L = T R gives the maximum heating power leakage that the refrigerator can sustain. If the heat leakage is above the cooling power at the thermal equilibrium, no net refrigeration can be accomplished. The heat absorbed per cycle has the same order of the work per cycle, ∼0.1eR 0 I c .
It is possible to find an analytic expression for Q cs valid for T R , T L ∆ 0 . Considering the scheme in Figure 9b, it can be noticed that at low temperature the heat absorbed by the CS is ruled by the purple area defined by the linear expression of entropy in Equations (23) and (24). Approximating the T 2 temperature to 0, due to the strong isentropic cooling at low temperatures, the Q cs is given at the leading order by the CS temperature This expression is plotted in Figure 10d as a violet dash-dotted curve. The agreement with the numerical results is good at T L < 0.2T c , corresponding to the agreement range in Figure 4b. From the characteristics of W, Q hr , Q cs in Figure 10 it is possible to calculate numerically the engine efficiency and the refrigerator COP. Figure 11 reports the efficiency and the COP for the studied Otto engine. Figure 11a shows a color plot of η, COP versus T L , T R . The two quantities are confined respectively in the engine and refrigerator regions of (T L , T R ). The gray area corresponds to the Cold Pump case. Figure 11b reports cuts of the efficiency η versus the Hot Reservoir temperature T R for chosen ambient temperatures T L . The curves end on the left in correspondence of the Otto characteristic curve, where the efficiency saturates at the Carnot limit. Figure 11c reports cuts of the COP versus the CS temperature T L for chosen HS temperatures T R , showing the evolution of the COP when the CS is cooled down toward the minimum achievable temperature, that delimits the COP curves on the left. The COP curves are limited on the right by the thermal equilibrium state T L = T R , where the COP reaches the theoretical Carnot limit. with different color palettes. The gray region represents the state where the cooled subsystem temperature is above the heat sink temperature. (b) Cuts of Otto cycle efficiency η versus T R for chosen T L in legend. The dot-dashed line reports the Carnot limit to efficiency. The curves end at the Otto characteristic curve, Equation (53), where the efficiency reaches the Carnot limit. (c) Cuts of Otto cycle COP versus T L for chosen T R in legend. The dot-dashed line report the Carnot limit to COP. The curves are limited on the right by the thermal equilibrium state T L = T R ; on the right, the curves are limited by the Otto cycle characteristic curve. On this curve, the COP reaches the COP Carnot limit.
An interesting property of the Josephson-Otto cycle is that close to the characteristic curve, both η and the COP reach the Carnot limit, even though the work released or the heat absorbed goes to zero. This point can be explained by referring to the degenerate case of Figure 9a. Close to the characteristic curve, the quantities Q L , Q R tend to zero but their ratio tends to |Q L /Q R | → T L /T R . Exploiting the energy conservation Q L + Q R − W = 0, it is that is the Carnot limit. With similar considerations, we obtain the analogous limit for the COP:

Josephson-Stirling Cycle
In this section, we analyze another possible thermodynamic cycle, i.e., a Josephson-Stirling cycle, that has in practice several practical applications, in particular as refrigerator [114]. The Stirling cycle is a different combination of the studied processes, being built with two isochorics and two isophasics. In an ideal gas system, it consists of two isochoric heat addition/rejection processes and two isothermal (compression + expansion). Real Stirling engines are eventually equipped by regenerators that increase the efficiency [115,116]; here we study the simple case without the regenerators.
First of all, let us consider the engine case and then move to the refrigerator one. The Josephson-Stirling engine is described by the scheme in Figure 12, where panels a and b show respectively the processes in the ST diagram and I ϕ diagram. The cycle is constituted by two isothermal processes (1 → 2 and 3 → 4 in Figure 12) and two isophasic processes (2 → 3 and 4 → 1 in Figure 12). The states 1, 2 and 3, 4 are respectively thermalized to the right and left reservoirs. When operating as Stirling engine, the left and right reservoirs play the roles of ambient and heat source respectively. In summary, the Josephson-Stirling engine is given by the succession of the following processes: • Isothermal 1 → 2. The thermal valves v R is open and v L is closed, so that the system is in thermal contact with the right reservoir. The system is driven from the state (ϕ 1 = 0, T 1 = T R ) to (ϕ 2 = π, T 2 = T R ). Here a work is spent |W 12 | represented by the green area in Figure 12b. The heat Q 12 is absorbed from the reservoir, represented by the green + dark purple area in Figure 12a. • Isophasic 2 → 3. By closing v R and opening v L , the system goes from the state (ϕ = π, T 2 ) to (ϕ = π, T 3 = T L ). The system releases heat Q 23 to the left reservoir, represented by the light purple + dark purple area. No work is performed, W 23 = 0.
• Isothermal 3 → 4. The valves are kept in the same state: v R open and v L closed. The system is driven from the state (ϕ 3 = π, T 3 = T L ) to (ϕ 4 = 0, T 4 = T L ). In this process the system returns a work W 34 represented by the sum of the green and blue areas in Figure 12b. The heat |Q 34 | is released to the left reservoir, represented by the blue area in Figure 12a. • Isophasic 4 → 1. By closing v L and opening v R , the system goes from the state (ϕ = 0, T 4 ) to (ϕ = 0, T 1 = T R ). The system absorbs the heat Q 41 from the reservoir at T R , given by the sum of the areas in blue, red and light purple in Figure 12a. No work is performed, W 41 = 0.
The total work per cycle is given by W = W 12 + W 34 . The heat absorbed from the Hot R (represented by the R reservoir) is In order to work as an engine, it must be T R > T L , as shown in Figure 12a. If T L > T R , the cycle is reversed as displayed in Figure 13. Panels (a) and (b) show the case of T R = 0.6T c and T L = 0.35T c and T L = 0.25T c respectively. In this case, the machine can work as a refrigerator with the CS represented by the R reservoir and HS represented by the L one (differently from the case of the Josephson-Otto cycle). Figure 13. Particular examples of the Josephson-Stirling cycle for T R < T L . (a) Stirling inverse cycle working as refrigerator. The heat absorbed from the R reservoir in the process 1 → 2, represented by the area defined by the related green arrow, is bigger than the heat released to R reservoir in the process 4 → 1, represented by the area defined by the related red arrow. (b) Stirling inverse cycle working as Joule pump, exploiting work to release heat to both reservoirs.
Even though the cycles in both panels are clockwise, only the cycle in panel (a) works as refrigerator. Indeed, there are further conditions that define the (T L , T R ) region where the cycle can work as a refrigerator. Let us consider the heat exchanged with the cold right reservoir, given by processes 4 → 1 and 1 → 2. From Figure 13 it can be noticed that in 4 → 1 the heat is released from the system to the R reservoir, while in 1 → 2 the heat is absorbed by the system. Cooling then can take place if Q cs = Q R > 0, i.e., if |Q 12 | > |Q 41 |. This is true when T L is closely below T R , T L T R ; then, when the CS is cooled down, |Q 12 | decreases, since δS(ϕ = π, T) in an isothermal heat exchange (31) decreases, while |Q 41 | increases with the increase of the temperature difference T hs 4 − T cs in the isophasic process. As a consequence, it exists a minimum achievable temperature T MAT that is characterized by a null cooling power Q R = 0, i.e., Note that T MAT is a function of the HS temperature T L . If T R < T MAT , the total heat Q R exchanged with the R reservoir is negative, and the CS is heated. This case corresponds to the (T, S) diagram in Figure 13b, where the red area representing the released heat to the right reservoir includes the green area of the absorbed heat from the right reservoir. As before, we call the curve (T L , T R = T MAT (T L )) the characteristic curve of the Josephson-Stirling cycle. We observe for completeness that when T R < T MAT and Q R < 0, also the left reservoir can absorb or release heat Q L = Q 23 + Q 34 , depending on (T L , T R ). If Q L > 0 the cycle absorbs work to transfer heat from the hot to the cold reservoir, constituting a Cold Pump similar to the situation described in the Josephson-Otto cycle. On the other hand, if Q L < 0, the machine releases heat to both the reservoirs, converting completely the work in heat. Following the definition of References [112,113], we call this operating mode as Joule pump.
In the refrigerator case, the total work is W = W 12 + W 34 and the heat extracted is Q cs = Q R = Q 12 + Q 41 , like the engine case.
The released work W and the heat absorbed Q cs , Q hr are summarized in Figure 14. In the color plots in panels a,b, the curves W = 0, Q R = 0, Q L = 0 separate the regions of the engine, the refrigerator, the Joule pump and the cold pump. The curve W = 0 corresponds to T L = T R . The refrigerator region is between the curve T L = T R and the characteristic curve T MAT (T L ).  Figure 14c reports cuts the released work per cycle versus the HR temperature T R for fixed CR temperatures T L . The curves reach zero at T L = T R . The general trend is that the work increases with increasing the temperature difference T R − T L between the two reservoirs.
An analytical expression for W can be calculated. Let us consider a Stirling cycle with T L T R ∆ 0 . The released work can be approximated by W = E (ϕ = π, T R ) − E (ϕ = π, T L ≈ 0). Using approximation (33) for E , we obtain This expression is plotted in Figure 14c and is in good agreement with the numerical results. Figure 14d reports the heat absorbed per cycle Q cs = Q R versus the CS temperature T R for fixed HS temperatures T L . The curves go to zero on their left in correspondence of the characteristic curve. The curves are limited on the right by the black curve of Q cs at T L = T R . The order of magnitude of the absorbed heat per cycle is ∼ 0.1eR 0 I c .
From the W, Q hr , Q cs characteristics it is possible to calculate the η and the COP, as reported in Figure 15. Figure 15a shows a color plot of η and COP versus (T L , T R ). The two quantities are plotted over the engine and refrigerator regions respectively. The gray area is where the cycle works as cold pump or Joule pump.  Figure 15b reports cuts of the efficiency η versus the hot reservoir temperature T R for chosen ambient cold reservoir temperatures T L . The curves end on the left at the state T R = T L . For T L → T R , both the cycle efficiency and Carnot limit tend to zero. Indeed, the work W = Q = Q R + Q L tends to zero, since Q R → −Q L as can be noticed from Figure 13a. Figure 15c reports cuts of the COP versus the CS temperature T R for chosen HS temperatures T L . The curves end on the left at T R = T MAT . For T L → T R , both the cycle COP and its Carnot limit tend to infinity. With the same geometrical argument used for the efficiency, it is W → 0 and Q R = T R δS(ϕ = π, T), implying that the COP = Q R /|W| → ∞.

Experimental Feasibility
Here, we briefly comment on some experimental aspects that have to be considered to implement and measure thermodynamic quantities discussed above, based on hybrid junctions. The two crucial assumptions of this paper are: (i) the processes are quasi-static and (ii) the system is thermally isolated.
As usually done in the thermodynamics, one is interested to investigate the performance of a thermodynamic cycles in the adiabatic limit (slow evolution) in order to avoid any additional irreversibility due to non-equilibrium processes. Anyway, in any practical realization, one need to develop a cycle in a finite time so it is fundamental to discuss which are the fundamental timescales for the validity of the quasi-static assumption. Hereafter we discuss the quasi-static assumption. It puts a limit on the speed of a process and hence to the cycling frequency ν. In our system, the equilibrium is determined by thermalization of the electron system, and hence the time of equilibration is set by the electron-electron thermal relaxation time τ e−e . Non-equilibrium experiments have been performed in superconducting systems for probing the time-scales of thermal relaxation [117][118][119]. For aluminium, τ e−e ≈ 1 − 10 ns close to T c and increases by decreasing the temperature [117]. Below T 0.1T c it has been measured that τ e−e saturates at ≈10 2 − 10 3 µs [118]. On the other side, a material with very low τ can be niobium nitride, where the electron-phonon relaxation time is 200 ps [119], suggesting a the same order of magnitude for the electron-electron relaxation time. The minimum time interval for a process to be quasi-static can range from the milli-second to hundreds of pico-second. Hence, we conclude that the rate at which a process or cycle can be performed depends largely on the material and temperature ranges and it is a fundamental issue related to the specific device realization. Even the presence of impurity scattering can alter the relaxation time [120]. Experiments suggest that this rate can range from the KHz to tenths of GHz [117][118][119].
The second assumption concerns the thermal insulation. In a superconductor, like in any metal, the electron system is in thermal contact with the phonon system. In the superconducting case, the heat flow is exponentially suppressed at low temperatures and scales like the volume V of the device [121][122][123]. Since the electron-phonon thermal conductance is an intrinsic property of the metal and can not be avoided, it is relevant to estimate a thresholdQ th to the heat leakage below which an isentropic process can be observed. In order to calculate this quantity, let us consider a process that drives the phase from ϕ = 0 to ϕ = π in a time τ. In the ideal case of an isolated system, the process is isentropic. In the real case, a certain amount of heat will be absorbed due to the closure of the minigap. If the electron system is well thermalized with the phonon system (corresponding to most of real cases), the process is isothermal and absorbs an average heat poweṙ where Q is the heat exchanged during the isothermal process, given by Equation (32). As a consequence, isentropic effects can be observed if the heat leakage of the electron system allows a power flow that is negligible respect to (61). In particular, for fast processes driven at the quasi-static limit, the thresholḋ Q th is maximized toQ th,m :Q If the heat leakage of the system is above this value, an isentropic process can not be observed. At T = 0.1∆ 0 and τ ≈ 10 ns it isQ th ≈ 3 × 10 −8 WA −1 I c . For I c ≈ 1 mA, the threshold isQ th = 30 pW. Based on these considerations, a first experimental setting for testing our findings consists of the measurement of heat capacity for different values of phase difference ϕ. In this case, no external reservoirs neither heat valves are required. The measure can be performed by heating up the device with a fixed amount of heat Q test (through a fixed power pulse) and subsequently measure the temperature increase ∆T. The quantity C = Q test /∆T, that approximates the heat capacity C(ϕ, T), is dependent on the phase ϕ according to the results of Section 3.2. In particular, the relative difference of heat capacity is C(ϕ = π, T)/C(ϕ = 0, T) ∝ α(T/∆ 0 ) 5/2 e T/∆ 0 , that is strongly enhanced at low temperatures and is proportional to the parameter α ∝ I c /V. Interestingly, for these experiments a perfect thermal isolation is not a crucial task, even though it would make the effect more evident. The presence of a certain thermal leakage would pull the device to the bath temperature after the power pulse, making the thermometry more difficult. In this case, a more complete thermal model is necessary to describe the device behavior. Another possibility is based on the measurement of temperature variation during an isoentropic process. In this case, no external reservoirs neither heat valves are required, but the thermal isolation is a crucial element. As discussed above, the heat leakage in a process ϕ = 0 → π in a time interval τ must be negligible respect toQ th in Equation (61). The heat leakage threshold can be increased by speeding up the process, i.e., reducing τ toward the limit τ e−e .
Finally, the most challenging but direct experiment is the cooling of a subsystem with a refrigeration cycle. In this case, external reservoirs, heat valves and the thermal isolation of the whole device (ring, valves, heat channels, cooled subsystem) are required. In detail, in order to observe the cooling effect, the spurious heat leakage must be lower respect to the cooling power at the thermal equilibrium. As an example, let us consider the Otto cycle removed heat per cycle in Equation (55). The cooling power is πν(T/∆ 0 ) 2 eR 0 I c /3κ, where ν is the cycling frequency. For ν = 100 MHz, T = 0.1T c , R 0 = 12.9 kΩ we have a cooling power of 0.2 pW.
Before concluding, a comment on the presence of defects at the SN interfaces is in order. Indeed, in SNS systems, the main source of defects is the quality of the SN interfaces. The opacity of the SN contacts can make the proximity effect weaker till disappearance in the tunnel limit [83,106]. However, several experiments on quasi-particle tunneling in SNS junctions show that the induced mini-gap has reached good quality features over time [89,103], suggesting the possibility of making good entropy variations in such devices. Finally, opaque interfaces imply that the CPR does not follow a KO form, so a different relationship between phase and entropy is expected. However, the Maxwell consistency relation is universal and must hold for every kind of CPR and the entropy variation can be obtained similarly as we proceeded for non-clean interfaces.

Conclusions
In this paper, we have reviewed and analyzed a proximized system with a phase-biased SNS junction under the thermodynamic point of view. By means of arguments of thermodynamic consistency, we have obtained the phase-dependent entropy of the system from its current-phase relation, that we assumed to a Kulik-Omel'yanchuk form. The entropy phase-dependence is related to the presence of an induced minigap in the density of states of the weak link; the minigap depends on the phase ϕ across the junction, yielding the phase dependence of the available states and hence of entropy. We obtained closed-form expressions of the entropy for low temperatures and ϕ = 0, ϕ = π. These expressions evidence a strong difference in the temperature dependence for the two phases, where the former is exponentially suppressed and the latter is linear. The entropy relative variation on the phase difference ϕ scales like the ratio of the critical current over the system volume. Hence, a stronger effect requires a higher critical current or smaller volume of the system.
We have discussed equilibrium thermodynamic quantities under quasi-static conditions, obtained by means of the Maxwell equations, investigating processes where phase-coherent properties are linked with thermal properties. In detail, this approach envisions two particular physically observable effects. First, the heat capacity of a proximized system is phase-dependent, passing from an exponentially suppressed behavior at ϕ = 0 to a linear behavior at ϕ = π. Second, the electronic temperature is subject to an isentropic cooling if the system is kept thermally isolated.
Finally, guided by the analogy of the Josephson-base thermodynamics with the classical thermodynamics, we discuss different kind of thermodynamic processes such as isothermal, isophasic and isentropic. After we combine these transformations to define the Josephson-Otto and the Josephson-Stirling cycles, which combine quantum coherence and Josephson effect. This requires the system to be connected to two different reservoirs through heat valves that can allow or stop the flow of heat from them. We characterized the cycle performances in terms of efficiency and COP. The Otto cycle, in particular, shows an interesting capability of having a cooling power till sub-milli-kelvin temperatures.
Further developments can be argued, including Processes that involve a current bias of the junction can be studied. In this case, the phase is a function of the temperature and iso-current processes can be conceived, that can be exploited in different cycles, like Brayton or Diesel cycles.

Conflicts of Interest:
The authors declare no conflict of interest. The funders had no role in the design of the study, in the writing of the manuscript, or in the decision to publish the results.

Appendix A. Thermodynamics Close to the Critical Temperature
In this appendix we briefly discuss a particular limit of our theory when the system temperature approaches the critical temperature. Apparently, our approach would bring to a thermodynamic inconsistency at temperature close to the critical temperature. Let us indeed consider the entropy scheme in Figure A1a. The red solid curve and the dashed green curve report respectively the quantities S(ϕ = 0, T) = S 0 (T) and S(ϕ = π, T) = S 0 (T) + δS(ϕ = π, T). For T > T c , the entropy is given by the normal metal form S N (T) = 2π 2 N S VT 3 (A1) plotted over the whole temperature interval as reference (blue dashed line). As can be noticed from Figure A1a, S(ϕ = π, T) has a decreasing jump ∆S = S(T → T + c ) − S(ϕ = π, T → T − c )at T = T c , corresponding to a first-order transition with negative latent heat T c ∆S. According this scheme, the system releases heat from the superconducting state to the normal state, that is unphysical. Figure A1. (a) Entropy scheme close to the critical temperature.S is equal to S where a phase transition is imposed at T c2 , calculated with the method in the text. S N (T) is the normal metal entropy. (b) Dependence of the critical temperature of the system T c2 at ϕ = π, normalized to the bulk critical temperature T c , versus α. This problem has been discussed and solved in References [73,74]. The unphysical states close to the transition are metastable and do not correspond to the minimum of the system free energy. The physical stable state is instead given by a null order parameter, that is the normal metal state. The physical solution, instead, consists in a first-order phase transition at a critical temperature T c2 , lower than the bulk critical temperature T c , that brings the system from a normal metal state to a superconducting state with absorption of latent heat. An experiment returned results in agreement with this theory [75].
Following the guidelines of [73,74], we discuss here a simplified macroscopical approach that yields physically acceptable result and allows to grasp the underlying physics.
Let us consider the free energy of both the normal state F N (T) and superconducting state F S (ϕ, T). At (ϕ = 0, T c ) the two free energies are equal, F S (ϕ = 0, T) = F N (T) = F 0 . In their neighborhood on the plane (ϕ, T), it is According the minimization of the free energy, the device is in a superconducting phase in the region (ϕ, T) where F S (ϕ, T) − F N (T) < 0. The boundary between the two region defines critical temperature T c2 (ϕ), that is a 2π periodic in ϕ with a minimum in ϕ = π. Over this boundary, the transition is of the first-order (but at ϕ = 0 where takes place the well known second order transition) with a jump in the entropy and associated latent heat. In Figure 4a, the dashed green line shows the entropy scheme for this first-order transition at ϕ = π for α = 0.6. In this case, the transition temperature is T c2 (ϕ = π) ≈ 0.7T c . The transition temperature depends on α, since it rules the ratio between S 0 and E in (A3). The dependence of T c2 (ϕ = π) versus α is reported in Figure 4b, normalized to the bulk critical temperature. As expected, T c2 (ϕ = π) → T c for α → 0. This calculation returns that the superconductivity is suppressed at α ≈ 1.5, setting practically an upper limit to the supercurrent density.