Phase-tunable thermoelectricity in a Josephson junction

Superconducting tunnel junctions constitute the units of superconducting quantum circuits and are massively used both for quantum sensing and quantum computation. In previous works, we predicted the existence of a nonlinear thermoelectric effect in a electron-hole symmetric system, namely, a thermally biased tunnel junction between two different superconductors, where the Josephson effect is suppressed. In this paper we investigate the impact of the phase-coherent contributions on the thermoelectric effect, by tuning the size of the Josephson coupling changing the flux of a direct-current Superconducting Quantum Interference Device (dc-SQUID). For a suppressed Josephson coupling, the system generates a finite average thermoelectric signal, combined to an oscillation due to the standard ac Josephson phenomenology. At large Josephson couplings, the thermoelectricity induces an oscillatory behaviour with zero average value of the current/voltage with an amplitude and a frequency associated to the Josephson coupling strength, and ultimately tuned by the dc-SQUID magnetic flux. In conclusion, we demonstrate to be able to control the dynamics of the spontaneous breaking of the electron-hole symmetry. Furthermore, we compute how the flux applied to the dc-SQUID and the lumped elements of the circuit determine the frequency of the thermoelectric signal across the structure, and we envision a frequency modulation application.


I. INTRODUCTION
The investigation of thermal transport in micro/nanoscale systems has attracted a growing interest in the last few decades [1][2][3][4][5][6][7][8], and is expected to have an impact on the performance of modern quantum technologies [9,10]. Heat dissipation is a key factor and limits the performance of classical computation platforms, but it is even more crucial in multi-qubit technology, where low operating temperatures further limit the heat exchange. Hybrid superconducting junctions [11,12] are ideal platforms for quantum devices [13][14][15][16], due to the well-established fabrication techniques and a precise modeling of the coherent electronic transport. In particular, they offer a tight control over thermal currents, with applications to electronic solid state cooling [2,3], phasecoherent modulation of thermal currents [17,18], and quantum sensing [19]. In the last few years, they have been also been extensively investigated for thermoelectricity, when superconductors are used in combination with ferromagnetic elements [20][21][22][23][24], and the interplay between thermoelectricity and the superconducting phase is being established [25][26][27][28][29]. In Refs. 30 and 31, we predicted an unexpected nonlinear thermoelectric effect occurring in a system with electronhole (EH) symmetry, paradigmatically a tunnel junction between two different Bardeen-Cooper-Schrieffer (BCS) superconductors. We observed that thermoelectricity arises due to a spontaneous breaking of EH symmetry determined by the electrode with the larger gap to have the higher temperature [30,31]. In the discussion, we focused only on the quasiparticle transport across the junction and we assumed to be able to suppress completely any phase-dependent contribution associated to the Josephson effect [12]. Purpose of this work is to investigate in details the impact of the phase-dependent * giampiero.marchegiani@nano.cnr.it † alessandro.braggio@nano.cnr.it ‡ francesco.giazotto@sns.it terms on the thermoelectric behavior. As we will show below, the generation of a finite thermoelectric voltage is still possible in the presence of Josephson contributions, but an additional oscillating behaviour is generated in accordance to the Josephson effect. We demonstrate that the spontaneous generation of a thermoelectric voltage/current can be controlled by tuning the strength of the phase dependent terms, which, for the setup we consider, can be modulated by changing the magnetic flux in a Superconducting Quantum Interference Device (SQUID) [12,32]. Moreover, we discuss the impact of the Josephson terms on the whole dynamics of the system. In particular, the frequency and the amplitude of the thermoelectricinduced oscillations are numerically computed, and approximate expressions are obtained in some limiting cases.

II. MODEL AND RESULTS
A. Charge transport and thermoelectricity in a superconducting junction The charge current in a tunnel junction between two superconductors depends both on the phase bias (ϕ) and the voltage (V ) applied to the junction, as first predicted by Josephson [12,33,34] I(V, ϕ) = I qp (V ) + I j (V ) sin ϕ + I int (V ) cos ϕ, where I qp is the quasiparticle contribution, I j is associated with Cooper pairs tunneling, and I int gives the interference contribution associated with breaking and recombination process of Cooper pairs on different electrodes of the junction [35]. The explicit expressions of I qp , I j , I int are given in Appendix A, being well-known in the literature [12,34,36]. The current obeys the symmetry I(V, ϕ) = −I(−V, −ϕ). Thus, I qp , I int are odd in V and represent the dissipative (or active in the presence of thermoelectricity) components of the current, whereas the function I j (V ) is even and corresponds to a purely reactive contribution [37,38]. Indeed, in the presence of a phase bias, the junction can support an equilibrium (nondissipative) current even for V = 0, I = I c sin ϕ (dc Josephson effect), where I c = I j (V = 0) is called critical current. At a finite voltage V = 0, the phase across the junction oscillates in time according to the Josephson equation (AC Josephson effect) whereh is the reduced Planck constant, and −e is the electron charge. Namely, for a constant bias V (t) = V 0 the phase dependent terms oscillates in time with Josephson frequency f j = |V 0 |/Φ 0 and zero average value, where Φ 0 = πh/e ∼ 2 fWb is the flux quantum. In this case, the dc response of the junction is given by the quasiparticle contribution only. For our purposes, we consider a junction between two different superconductors (S, S ), and introduce an asymmetry parame- is the zero-temperature order parameter, and T c,i is the critical temperature of the i electrode. Figure 1a displays the voltage dependence of the three contributions when the electrodes have equal temperatures T S = T S = T and r = 0.75.
In the low temperature limit T T c,S (dashed lines), both I qp and I int are strongly suppressed for |V | < (∆ 0,S + ∆ 0,S )/e and finite at higher voltage. Note that I qp is positive and monotonically increasing for V > (∆ 0,S + ∆ 0,S )/e, where it asymptotically reads I qp = G T V (G T is the normal state conductance). On the other hand, I int is negative and monotonically decreasing for V > (∆ 0,S + ∆ 0,S )/e. The Cooper pairs term I j has a more complex evolution, and it is monotonically increasing for 0 < V < (∆ 0,S + ∆ 0,S )/e, where it diverges (the divergence is smoothed with the introduction of a finite phenomenological parameter Γ, see Appendix A), and it is monotonically decreasing at higher voltages. In the same plot, we display the evolution also for a finite value of the temperature, i.e., T = 0.6T c,S . While, the overall behaviour of the curves is similar, I qp , I int are now finite and displays a positive conductance G i = I i /V > 0, with i = {qp, int} at subgap voltages |V | < [∆ S (T )+∆ S (T )]/e, showing a nonlinear evolution characterized by a peak for V = ±V p = ±|∆ S (T ) − ∆ S (T )|/e, due to the matching of the BCS singularities in the density of states of the two superconductors [12]. With increased temperature, the Cooper pairs term is reduced due to the monotonically decreasing evolution of the superconducting gaps [∆ i (T ) ≤ ∆ 0,i ]. For the same reason, the voltage value where the Cooper pairs term has the peak and the other contributions have a sharp jump, i.e., V = [∆ S (T ) + ∆ S (T )]/e is reduced with respect to the low temperature limit.
Since we are interested in the description of thermoelectric phenomena, we consider a situation where a temperature difference is established between the electrodes, namely T S = T S . Note that in the absence of a temperature bias, i.e., T S = T S , the behaviour of the junction is purely dissipative, since (I qp + I int cos ϕ)V > 0 for every ϕ [38], as required by the second law of thermodynamics [30]. Conversely, with a thermal gradient it is possible to have a thermoelectric power generation with positive entropy production. In particular, a thermoelectric behaviour is characterized by a positive thermoelectric powerẆ = −IV > 0 produced by the junction. As discussed above, this definition mainly applies to the even-ϕ component of the current, i.e., I qp + I int cos ϕ, since the Cooper pairs term I j sin ϕ essentially describes a reactive component. In Ref. 30, we predicted the existence of a thermoelectric effect at subgap voltages in the quasiparticle component, namely we demonstrated that we can have I qp V < 0 for small values of V , provided that the superconductor with the larger gap is heated up [in our notation T S > T S and ∆ S (T S ) > ∆ S (T S ), since we assumed ∆ 0,S > ∆ 0,S ]. This is shown in Fig. 1b, where the subgap evolution of I qp , I int , I j is displayed for T S = 0.7T c,S and T S = 0.01T c,S . In particular, at a low voltage the quasiparticle curve displays a peculiar negative conductance G qp (V ) = I qp (V )/V < 0 and hence finite thermoelectric powerẆ = −I qp V > 0. In the absence of phase-dependent terms, the negative value of G qp for V → 0 implies a spontaneous breaking of electron-hole symmetry. This leads to the generation of a thermoelectric voltage due to the existence of finite voltage values ±V S where the current is zero [I qp (V S ) = 0], as discussed in Refs. 30 and 31.
Consider now the phase dependent terms. Interestingly, the interference term (green) behaves similarly to the quasiparticle term, also showing a negative conductance G int (V ) = I int (V )/V < 0 around the origin [39]. In particular, the zero-bias value of the differential conductance reads G 0,int = G 0,qp ∆ S (T S )/∆ 0,S in the limit T S → 0 (see Appendix A). The Cooper pairs term (yellow) is roughly constant for V < V p = [∆ S − ∆ S ]/e, where it sharply decreases, and rise monotonically at higher voltages. Similar jumps are observed also in the temperature evolution of the critical current [40]. Note that the size of the Cooper pairs term, which is finite at V = 0, is quite larger than the quasiparticle contribution. As a consequence, it is reasonable to expect the Josephson current to have a potentially crucial impact on the features of the thermoelectric effect.

B. Circuit dynamics modeling
In order to describe the impact of the phase-coherent contributions on the thermoelectricity of the junction, we consider a system where the size of I j , I int can be externally tuned. More precisely, we investigate the circuit displayed in Fig. 1c. The system features a superconducting ring made of two different superconductors, which are coupled by two tunnel junctions. This configuration is known in the literature as the direct-current Superconducting Quantum Interference Device (dc-SQUID) [12,32]. We assume that the superconducting ring is connected to an external circuit, which constitutes an idealized model for the electrical environment described in terms of lumped elements, i.e., an inductance L and a load R. Each of the two junctions displays the nonlinear currentvoltage characteristic I(V, ϕ) of Eq. (1), so that the total current in the dc-SQUID reads where ϕ 1 and ϕ 2 are the phase-differences across the two junctions, and α = G T 2 /G T 1 is the ratio between the conductances of the two junctions in the normal state. For simplicity, in the theoretical modeling we will consider a fully symmetric SQUID (α = 1), even if the results can be properly extended to asymmetric junctions [41]. For a proper description of the dynamics, we need to consider also the capacitance C of each of the two junctions. Due to the ring geometry, the superconducting phase differences across the two junctions are related by the fluxoid quantization, namely, where Φ is the total flux out of plane of the superconducting ring. The magnetic flux Φ coincides with the flux applied externally, since we assume the self-inductance of the ring negligible. By minimizing the free energy of the SQUID with respect to the superconducting phases ϕ 1 , ϕ 2 , and using the constraint of the flux quantization Eq. (4), we can rewrite the expression of Eq. (3) as whereφ = (ϕ 1 + ϕ 2 )/2. The circuit dynamics is finally expressed by [42]      which is an autonomous non-linear system of differential equations in the three variablesφ,V, I L . The first equation in Eq. (6) gives the current conservation in the circuit: the current I L which flows in the inductor L and in the load R is the sum of the current in the capacitances (first term in the right side) and in the two junctions (second term in the right side). The second equation in Eq. (6) gives the Kirchhoff voltage rule in the circuit: the voltage V across the SQUID is equal to the sum of the voltage drops in the inductance and in the load. The last identity in Eq. (6) follows from the Josephson relation between the phase bias and the voltage bias in a Josephson junction of Eq. (2). As can be seen from Eq. (5), the absolute strength of the phase-dependent contributions can be fully tuned by varying the magnetic flux Φ, as shown in Fig. 1d for the Cooper pairs term I j at V = 0. In particular, the evolution of the phase-coherent contributions is periodic with period Φ 0 : the Josephson current is maximum for Φ = nΦ 0 (with n ∈ Z) and it is exactly zero for Φ = (1/2 + n)Φ 0 [43]. For this reason, we will consider the evolution only in a single period Φ ∈ [0, Φ 0 ].

C. Flux modulation of the dc-thermoelectricity
We wish to investigate if the thermovoltage can be generated in the presence of the phase-dependent terms. First, we consider situations where these terms are suppressed, which mainly happens for Φ ∼ (n+1/2)Φ 0 (with n ∈ Z). Indeed, for Φ = (n+1/2)Φ 0 the phase-coherent contributions are zero. In this case, the dynamics of the variables V, I L is independent oñ ϕ. This limit corresponds to the one previously discussed in Ref. 30. In particular, the stationary time-independent solutions are obtained by solving the implicit equation [44] For a dissipative junction, I qp V > 0, and the only solution of Eq. (7) is V = 0 (and thus I L = 0). However, in the presence of a thermoelectric effect, the behaviour of the system depends on the size of the load R [30,31]. That is, for R < V p /(2I p ), the system may display a oscillatory behaviour with zero average value of I L and V . Conversely, for R > V p /(2I p ), the system admits stationary and time-independent solutions (V, I L ) = (V , 2I qp (V )), where |V | > V p . Note that, due to EH symmetry, each positive solutionV > 0 has a corresponding solution −V < 0. As a consequence, the system approaches eitherV or −V in the steady-state, depending on the specific initial condition [30].
Here, and in the rest of this work, we a set of realistic parameters, for an aluminum-based SQUID with T c,S = 1.6 K (and thus ∆ 0,S = 1.764k B T c,S ∼ 240µeV) and G T = (1kΩ) −1 .
We consider the thermoelectric situation displayed in Fig. 1b for scaled quantities, where r = 0.75 and V p ∼ 0.08∆ 0,S 19µeV. Figure 2a displays the absolute value of the stationary value of the thermoelectric voltage V 0 = |V | (solid) as a function of the load R computed through numerical solution of the system of Eq. (6) (see the discussion on Sec. II C 2) and the thermoelectric voltage obtained by solving the implicit equation Eq. (7) (dashed). Note that the two quantities coincide, except for a very narrow range 150Ω ≤ R ≤ 200Ω where the solution of the implicit equation is different from zero, while the result of the numerical integration is zero. This small different behaviour is associated to the stability of the V = 0 solution of Eq. (7) and will be discussed in more details after. As discussed above, the thermoelectric voltage is zero for low values of R and it is finite (and larger than V p ) and monotonically increasing for R > V p /[2I(V p )] ∼ 150Ω. In the limit R → ∞, the thermoelectric voltage approach the Seebeck voltage V S , i.e., the zero-current solution I(V S ) = 0 with finite voltage bias V S = 0.

Small Josephson contribution
In the presence of a small Josephson current, the picture described above is expected to be slightly modified. Indeed, in the presence of a finite voltageV , the phase evolves in timeφ(t) ∼ 2eV t/h due to the ac Josephson effect and so an oscillating term δV (t) with characteristic frequency fV =V /Φ 0 is superimposed to the dc thermoelectric voltage, i.e., V (t) ∼V + δV (t). In order to compute the perturbative contribution δV (t), we consider the first equation in Eq. (6), In the leading order of the perturbative expansion, I L ∼ −V /R and we can approximately set |I int (V )|, we can neglect the first term in the right side of the equation in first approximation, and obtain by integration is the amplitude of the Josephson current of the SQUID. The validity of the perturbative solution is good when the size of the correction is much smaller then the leading term, i.e.,hI j (Φ)/(4e|V |C) |V |. In terms of the Josephson current suppression, the previous relation requires The typical thermo-voltage is of orderV ∼ 0.1∆ 0,S /e, whereas the critical current is roughly I j (Φ = 0) ∼ G T π∆ 0,S /e , giving I j (Φ)/I j (Φ = 0) 0.01C∆ 0,S /(G Th ). Interestingly, this last inequality shows that the requirement on the Josephson coupling suppression to generate an average thermoelectric signal depends on the superconducting gap but not necessarily on the geometric area of the junction, since both C and G T are proportional to the area of the junction. For an aluminum based structure, characterized by ∆ 0,S ∼ 0.2 meV, specific capacitance of the barrier C/A = 50 fF/µm 2 and specific conductance G T /A = 1 mS/µm 2 , one obtain I j (Φ)/I j (Φ = 0) 0.15 in the worst case scenario. This means that sometimes a moderate suppression of the Josephson coupling is sufficient to generate a dc thermovoltage.

Flux evolution of the circuit dynamics
Now we wish to discuss the crossover between Φ = 0, where the phase-coherent contribution is maximum, to the case Φ = 0.5Φ 0 , where it is zero. We consider a few cases where R > V p /[2I(V p )], thus the load is large enough to induce a spontaneous symmetry breaking in the absence of Josephson contributions (see colored points in Fig. 2a). We solved numerically the system of Eqs.6: after a transient evolution, the steady-state solution is periodic, with a period T which depends both on the load and on the SQUID magnetic flux. Figure 2b displays the average value of the voltage signal in the steady-stateV = 1/T t 0 +T t 0 V (t )dt as a function of the magnetic flux for the three values of the load considered. For all the curves, the mean voltage is finite and very close to the value at Φ = 0.5Φ 0 for values of the flux where the critical current is strongly suppressed. These results show that even in the presence of a small phase-coherent (Josephson) contribution, the junction may still generate a breaking of EH symmetry and a net dc thermoelectric contribution. On the other hand, for large values of the critical current the mean voltage drops to zero. Note that the critical value of the flux where thē V switch from zero to a finite value depends on various parameters, and in particular on the load resistance. For a large load, the dc thermoelectric voltage is present even in the presence of a moderate Josephson current. In order to characterize the dynamics more completely, we computed the frequency f = 1/T and the amplitude, defined as A V = [max t V (t)−min t V (t)]/2, of the voltage oscillations in the steady-state. Figure 2c displays the flux evolution of the frequency. In particular, the frequency decreases by reducing the critical current in the region where the average value of the voltage is zero, whereas it is larger and it is exactly proportional to the mean valueV of the oscillations whenV = 0, due to the ac Josephson effect. The corresponding amplitude of the oscillations is shown in Fig. 2d. In the region where theV = 0, the amplitude slightly decreases upon decreasing the phase-coherent contributions, and then decrease sharply to a small value in the proximity of Φ = 0.5Φ 0 , where the amplitude is well described by the prefactor of the cosine term in Eq. (9) (see dashed lines in Fig. 2d).

D. Load dependence
Here we give a more general discussion of the impact of the load R on the dynamics of the junctions and hence on the thermoelectric features. We will consider two extreme values of the flux: Φ = 0, where the phase contributions is maximum and Φ = 0.5Φ 0 , where it is minimum. Figure 3a displays the frequency of the steady-state oscillations as a function of the load resistor. The corresponding amplitude of the voltage oscillations is shown in Fig. 3b.

Small Josephson contribution.
Consider first the case Φ = 0.5Φ 0 , where the Josephson coupling is negligible. We can identify three main regions. For a small load, i.e., for R ≤ 10Ω, the voltage bias falls mainly across the inductor L, i.e., V (t) ∼ −Lİ L (t) and the system behaves as a LC oscillator of characteristic frequency GHz. The steady-state oscillations are characterized by a zero mean value of the voltage bias and a sizable amplitude which can be computed through the energy balance in the system. In particular, in the steadystate, the total energy dissipated in one cycle in the resistor (Joule heating) must be equal to the total energy provided by the superconducting junctions, In the steady-state, the previous equation is generally valid irrespectively on the strength of the Josephson coupling. In the case considered here (Φ = 0.5Φ 0 ), the SQUID current is given by the quasiparticle transport, I SQ (t ) = 2I qp (V (t )).
In general, the power generated by the thermoelectric effects in the junction is able to self-sustain an oscillatory behaviour in the steady-state. It is relevant to note that the thermoelectric power is given both by I qp and I int (if present). Assuming a quasi-sinusoidal oscillatory regime in the steady-state with zero average value, i.e., , the energy balance Eq. (10) yields an integralalgebraic equation in A V which can be numerically solved for each value of R. The result of this approximation is shown in Fig. 3b with a double-dotted dashed curve (lower) and describes very accurately the amplitude of the steady-state signal computed through the numerical solution of the differential equations in Eq. (6). At intermediate loads, i.e., 10Ω ≤ R ≤ 200Ω, the system shows a pure-dissipative behaviour and relaxes to a time independent zero-current state. We stress that in the range 150Ω ≤ R ≤ 200Ω, the junction also supports time-independent solutions with finite voltage [see dashed line in Fig. 2a, which represents the solution of the implicit Eq. (7)], and the timeevolution can be either dissipative and lead to the zero-current solution I L = V = 0 in the steady-state or may produce a finite dc-thermoelectric voltage, depending on the particular initial conditions.
For large loads, R > 200Ω, the system oscillates around the thermoelectric solutionV [see solid curve in Fig. 2a] with frequency f =V /Φ 0 and a small amplitude. The latter is better visualized in the inset of Fig. 3b, which shows a magnification of the amplitude, which for this example is of order 1µV. Note that the expression is well described by the coefficient of the cosine term in Eq. (9), displayed in the inset with a dashed curve. Thus, the frequency increases monotonically with the load [since the average thermoelectric voltageV is monotonically increasing, see Fig. 2a] and the amplitude is moderately decreasing.

Large Josephson contribution.
For a large value of the Josephson current, the evolution is qualitatively different. In particular, the frequency is monotonically increasing for low values of the load R < 20Ω, where it reaches a maximum f max ∼ 17 GHz, and monotonically decreasing at larger values. The amplitude of the oscillations follows the inverse pattern, with a monotonically decreasing evolution for R < 100Ω, and a growth at larger values. It is possible to compute numerically with a good degree of approximation the load evolution of the frequency and the amplitude of the voltage oscillations (without integrating explicitly Eq.(6)) both in the low-load limit (roughly for R ≤ 10Ω, see the upper double-dotted dashed curves in Fig.3a and in Fig.3b) and in the large load limit (R > 200Ω, see the upper dotted dashed curves in Fig.3a and in Fig.3b). Here we give the fundamental elements of the theoretical approach, and we leave a more detailed discussion of the modeling to the Appendix C.  In the low-load limit, the frequency of the oscillations is increased with respect to the zero Josephson coupling case. In fact, the system still behaves approximately as a LC oscillator, but with a modified inductance L eff , which is the parallel of the external inductor L and the Josephson inductance [46] L j = Φ 0 /[2π2I j (0)] ∼ 0.55 nH, namely L eff = (L −1 + L −1 j ) −1 ∼ 0.35 nH and characteristic frequency f ∼ 19GHz. Note that the actual value of the frequency is slightly smaller and dependent on the load R. This behavior is related to the nonlinear phase dynamics of the junctions, which is associated with a frequency and amplitude dependence of the effective inductance of the circuit L eff ( f , A V ), as shown in Appendix C. Therefore, the approximate expressions for f , A V (upper double-dotted dashed curves in Fig. 3) are obtained by solving self-consistently for f = [2π 2L eff ( f , A V )C] −1 and the energy balance Eq. (10) in the circuit. In the latter, one can obtain an integral equation in terms of f , A V by assuming a quasi-sinusoidal regime, similarly to the zero Josephson coupling case. Moreover, the Cooper pairs term I j (V ) plays no role in the energy balance since it is purely reactive, and only affects the effective inductance of the circuit L eff ( f , A V ), as discussed above.
For a large load, the voltage drop occurs mostly in the resistor I L (t) ∼ −V (t)/R, and we can write the first of Eq. (6) as a second order nonlinear differential equation in the phase-bias where the effective external forces are cosφ. (12) For |V (t)| = |Φ 0φ (t)/2π| V p , one can approximate I j (V ) I j (V = 0) (see Fig.1b), and Eq. (11) yields a damped pendulum equation with zero-amplitude angular frequency ω 2 0 = 1/(2L j C) and additional nonlinear terms which involves both damping and power generation in the presence of the thermoelectric effect. Unfortunately, in the typical situation considered here, the amplitude of the voltage oscillations is larger than V p and the approximation I j (V ) I j (V = 0) is inaccurate. Hence, the system behaves as a nonlinear pendulum where the zero-amplitude frequency depends on V ∝φ [in the mechanical pendulum analogue, the length changes during the evolution similar to an elastic string]. The right-hand side of Eq. (11) contains both the damping and the driving force of the nonlinear pendulum. In particular, in Eq. (12) the first term in the right-hand side gives the damping associated to the Joule heating in the load. The second term gives the quasiparticle current, which is active when |V | < V S and dissipative otherwise, where V S = 0 is the Seebeck voltage, [we recall that at the Seebeck voltage the quasiparticle current is zero I qp (V S ) = 0]. A similar behavior, active at low voltage bias and dissipative at higher voltage bias, applies also to the interference term.
In order to evaluate f , A V as a function of the load R, one has to solve self-consistently the energy-balance in the steadystate Eq. (10) [which involves the terms in Eq. (12)] and the relation between the frequency and the amplitude in the nonlinear pendulum [47], Above, K[k] is the complete elliptic integral of the 1st kind, Aφ = [max tφ (t) − min tφ (t)]/2 is the amplitude of the phase oscillations, and we replaced The theoretical modeling exploits an highly accuracy approximate solution of the nonlinear pendulum equation [48] which includes the effect of higher harmonics (see Appendix C for a detailed discussion), and perfectly describes the motion of the system (see the upper dotted-dashed curves in Fig. 3). The load evolution of f , A V can be qualitatively understood as follows. By increasing R, the dissipation in the circuit for a given voltage bias V (t) is reduced, since RI L (t) 2 ∼ V 2 (t)/R, producing an increase in the amplitude of the oscillations. As a consequence, the frequency of the oscillations decreases, since in the nonlinear pendulum the frequency is monotonically decreasing with the amplitude of the oscillations [see Eq. (13)].
Finally, we note that the behaviour of the junction is chaotic at intermediate values of R (see filled regions). In particular, the system may relax to a zero-current time independent solution, depending on the initial conditions. This can be understood by inquiring the eigenvalues of the linearized equations which describe the dynamics of the system close to the stationary solutions (see Appendix B).

III. CONCLUSIONS AND DISCUSSION
In summary, we have discussed the dynamics of thermally biased Josephson junctions, in the presence of the nonlinear thermoelectric effect recently predicted in tunnel junctions between two different BCS superconductors. We investigated a system where the size of the Josephson coupling can be externally tuned, by modulating the flux inside a SQUID. The system displays a rich phenomenology, when inserted in a generic electric circuit, such as a RL circuit. Depending on the load, we focused on two relevant different regimes.
In the presence of a large load, the system generates a finite dc-thermoelectric voltage when the Josephson coupling is strongly suppressed but still finite, due to the spontaneous breaking of EH symmetry. In addiction, the system outputs an ac signal with frequency exactly proportional to the thermoelectric voltage, due to the ac Josephson effect. As a consequence, both the thermoelectric voltage and the ac signal can be ultimately controlled by changing the size of the load. When the Josephson coupling is stronger, the system generates a pure ac-thermoelectric signal. When the load connected to the system is small, the systems generates an ac signal, independently on the strength of the Josephson coupling. Interestingly, the modulation of the Josephson current induces a control of the effective inductance of the circuit, and hence of the frequency of the thermoelectric signal. The operating ranges depend on the inductance connected to the circuit and are in the GHz regime for a standard aluminum based structures. We may envision different applications for this system, taking advantage of the different regimes. Firstly, we note that when the system generates a dc thermoelectric signal, one has an autonomous system that convert a temperature gradient in a dc voltage signal which is perfectly tuned (by the Joseph-son relation) with the frequency of the ac component. This may find some value when one need to have a controlled generator that need to be galvanically disconnected from external circuits. Another application may involve the detection of radiation: by tuning the system very close to the transition point where the mean thermoelectric voltage switches from V = 0 toV V p , the system is highly sensitive to small parameters variation, such as the load or the temperature difference. Therefore, events such as photon absorption may trigger the spontaneous breaking of EH symmetry. Finally, one may envision an application as an high frequency oscillator controlled by the flux and feed with a thermal gradient only. We believe that the discussed system presents novel properties and functionalities that can be relevant in the field of superconducting quantum technologies.
where ℜ[. . . ] and ℑ[. . . ] denote the real and the imaginary parts, respectively. In the BCS model, the quasiparticle den- i are the anomalous Green's functions (here i = S, S and j is the imaginary unit). We assumed the electrodes in the quasi-equilibrium regime, hence the quasiparticle distributions are the Fermi functions f i (E) = [exp(E/k B T i ) + 1] −1 , where k B is the Boltzmann constant. The phenomenological parameters Γ i (typically called Dynes parameters) give a phenomenological representation of the finite quasiparticle lifetime [49,50] or the influence of the electromagnetic environment of a tunnel junction [51]. In all the calculations, we set Γ i = 10 −4 ∆ 0,i . Equation (1) (with the expressions Eqs. (A1),(A2), and (A3)) are derived in the tunneling limit for a constant voltage bias V [12,34]. In the presence of a time dependent voltage, the expression of Eq. (1) does not hold generally anymore, and must be generalized to include time-delayed effects also [12,37]. However, in this work we consider the adiabatic regime [12,46], where we can still use the expression of Eq. (1), replacing V → V (t). The adiabatic approximation holds when the voltage signal is small eV (t) ∆ 0,S , or the time variations of V (t) are small compared to the gap frequency ∼ (1 + r)∆ 0,S /h [12,46]. We consider realistic values of the circuit parameters where both these conditions are reasonably fulfilled. We expect that the main predictions are not crucially affected even beyond the adiabatic approximations. As discussed in Refs. 30 and 31, the quasiparticle current shows a thermoelectric behaviour, i.e., I qp (V )V < 0 for T S > T S , provided ∆ S (T S ) > ∆ S (T S ). In the limit T S → 0 and for small values of the bias V → 0, the currents is approximately linear I int ∼ G 0,qp V , where the zerobias differential conductance reads [30,31] A similar expression can be derived for the quasiparticle interference term in the same limit, where I int ∼ G 0,int V , with Note that the ratio of the two quantities is given by the expression quoted in the main text, where the inequality holds in the thermoelectric regime, where ∆ S (T S ) ≥ ∆ 0,S .

Appendix B: Linearization and stability analysis
In order to describe the different regimes of the dynamical system, it is convenient to work in scaled units, namely we consider i i = eI i /G T ∆ 0,S (with the subscript i = {L, SQ}), v = eV /∆ 0,S ,Φ = πΦ/Φ 0 and τ = t/ √ 2LC. The frequency of the oscillations are obtained by multiplying the scaled frequencyω by f LC = (2π √ 2LC) −1 ∼ 11.3 GHz for our parameters choice. The system of Eq. (6) in scaled units reads (B1) Note that the dynamics of the system depends on three dimensionless parameters: ε, κ, ξ . More precisely, κ = 2∆ 0,S √ 2LC/h is the ratio between the gap frequency 2∆ 0,S /h and the angular frequency of the LC oscillations 2π f LC . As discussed in Appendix A, the validity of Eq. (1) [with the expressions Eqs. (A1), (A2), and (A3)] is restricted to the adiabatic regime, where the time variations are much smaller than the gap frequency, i.e., κ 1 [in our calculation, we set κ ∼ 10.5]. The other parameters are: ξ = G T R = R/R T , which is the ratio between the load and the normal state resistance R T and ε = G T L/(2C), which is proportional to the strength of the thermoelectric effect and thus also characterizes the coupling to the nonlinear terms of the system of equations. It is convenient to have a small value of ε, to avoid strong non-linearities in the dynamics [in the calculations, we set ε ∼ 0.07]. It is worth noting that the values of κ and ε adopted are obtained by considering realistic values for typical Josephson junctions realized through standard nanofabrication techniques. The stationary and time independent solutions are obtained by settingv =φ =i L = 0 and read v = i L = 0,φ = nπ (with n ∈ Z). The stability analysis can be inquired with a standard linearization procedure, which leads to the matrix equation is the scaled zero-bias differential conductance and i j,0 = i j (v = 0). In particular, a necessary condition for the stability of the stationary and time independent solutions is that the real parts of all the eigenvalues of the matrix in Eq. (B2) must be negative [52]. The eigenvalues λ can be obtained by solving the characteristic equation, obtained by setting det (M − λ I) = 0, where I is the 3 × 3 identity matrix. The explicit expressions of the eigenvalues in terms of the various parameters of the system are obtained by using the cubic formula (not shown here). Figure 4 displays the load evolution of the real part of the eigenvalues λ 1,2,3 as a function of the load resis-tor for the set of parameters used in the main text, both for odd values [ Fig. 4a] and even values [ Fig. 4b] of n. Note that, for odd values of n (Fig. 4a) the real part of λ 2 (solid red) is positive irrespectively of R. As a consequence, the stationary time independent solutions v = i L = 0,φ = nπ (with odd n) are always unstable. The situation is different for even values of n [ Fig. 4b].
In particular, the plot shows that for 60Ω ≤ R ≤ 230Ω (filled region) all the real parts of the eigenvalues are negative [the real parts of λ 2 and λ 3 coincide since λ 2 = λ * 3 in this case]. As a consequence, in this region the stationary time dependent solution characterized by v = i L = 0,φ = nπ (with even n) is stable, and the time dependent evolution of the system depends on the initial conditions, as we verified numerically by solving Eq. (B1) for different values of i L (τ = 0), v(τ = 0),φ(τ = 0). More precisely, if the system is slightly perturbed around the stationary solution, the evolution relaxes to it. For larger perturbations, the systems approach a limit cycle characterized by a periodic oscillation whose amplitude and frequency are shown in Fig. 3 of the main text. In the presence of a strong Josephson current, the electronhole symmetry breaking is only obtained in the timedependent domain, and both the mean value of the currentī and the mean voltagev are equal to 0. We focus on two differ-ent regimes, related to the value of the load. For simplicity, we consider the case of zero-flux, but the results can be extended to Φ = 0.

Small Load
In the presence of a small load, we look for a perturbative solution for the current in the circuit where i We are interested in the steady-state oscillatory evolution of the system, characterized by an unknown angular frequencỹ ω = 2π/T (hereT is the scaled period). By assuming a quasi-sinusoidal oscillation in the voltage, neglecting higher harmonics of the oscillations, one gets v(t) =Ã V sin(ωt) where we definedÃφ = κÃ V /ω. In order to computeÃ V ,ω, we insert these expression in the current conservation equation Eq. (C1), and obtaiñ where we defined ψ =ωt. We obtain two coupled equations through multiplication by either cos ψ or sin ψ and integrating over a period. We get ξ πÃ V 2ε 2ω 2 + where we have divided the either active/dissipative components of the current (related to i qp and i int , in phase with the voltage bias) by the reactive component i j (shifted by π/2 with respect to the voltage bias), exploiting the different symmetries inφ, v of the three contributions. Equation (C6) is related to the energy balance in the system, since at the steady-state the energy dissipated in the load during a period must be equal to the total energy produced in the junction for each cycle. In fact, it can be rewritten in general as which is exactly Eq. (10) in scaled units. Equation (C7) gives the relation between the frequency and the amplitude of the oscillation. The Josephson current affects the effective inductance of the circuit, and produces an increased frequency of the oscillatory behaviour with respect to the case where i j ∼ 0. The second term in the square brackets in Eq. (C7) can be interpreted as the frequency dependent correction of the circuit inductance due to the Josephson term (in units of 1/L). In fact, for small values of the phase-oscillationsÃ V , κÃ V /ω 1 (which is never properly met in our case), the integral gives a frequency independent result: 2εω A V π 2π 0 i j (Ã V sin ψ) sin Ãφ cos ψ cos ψdψ 2κεi j,0 = L L j (C9) where L j = Φ 0 /[2π2I j (0)] is the Josephson inductance, and we usedÃφ = κÃ V /ω. Finally, the amplitude and the frequency of the oscillation are obtained by solving selfconsistently Eq. (C6) and Eq. (C7).

Large Load
In the presence of a large load, we can neglect the voltage drop across the inductor and write i L (t) ∼ −v(t)/ξ . Upon substitution in the current conservation equation, we can write down a pendulum-like equation with self-forcing and dissipation [Eq. (11) in scaled units] (C10) As discussed in the main text, the mechanical analogue of this equation is a pendulum where the pendulum length depends on the phase derivativeφ and it is subjected to driving and dissipative forces (right-hand side of the equation). Since the sine term changes with time during the evolution, we replace i j (v) with a value averaged over the dynamics, i j (v) →ī j = We approximate the frequency by using its relation to the amplitude of the phase oscillationsÃφ by the well known expression for the standard nonlinear pendulum [Eq. (13) in scaled units] whereω 0 = 2κεī j and K[k] is the complete elliptic integral of the 1st kind. We verified numerically that the frequency of the steady-state oscillations is well described by this expression also for our case, upon inserting the values ofÃφ obtained from the numerical computation. The amplitude is still related to the energy-balance in the circuit in each cycle (C8). In order to properly describe the oscillations even for large values ofÃφ , we exploit the high precision approximate solution of the pendulum equation with initial amplitudẽ Aφ which is derived by differentiating Eq. (C12) and using v = ϕ/κ. Note that, in the sinusoidal limit, the second term in the square brackets is much smaller than one and the relation between the phase and the voltage oscillations reduces tõ A V Ãφω /κ, as in the analysis of the previous section. In the main text, we discussed how these approximate expressions compare with the solutions obtained through the direct numerical integration of Eq. (6).