Phase-coherent heat circulators with normal- or superconducting contacts

We investigate heat circulators where a phase coherent region is contacted by three leads that are either normal- or superconducting. A magnetic field, and potentially the superconducting phases, allow to control the preferential direction of the heat flow between the three-different temperature-biased contacts. The main goal of this study is to analyze the requirements for heat circulation in non-ideal devices, in particular focusing on sample-to-sample variations. Quite generally, we find that the circulation performance of the devices is good as long as only a few transport channels are involved. We compare the performance of circulators with normalconducting contacts to those with superconducting contacts and find that the circulation coefficient are essentially unchanged.


I. INTRODUCTION
With the miniaturization of electronic circuits and the possibility of fabricating devices operating at the nanoscale, control and management of heat flows [1] is becoming increasingly important. On the one hand, the performance of nanodevices can critically depend on ultracold temperatures [2]. Therefore, the ability to control heat flows is for instance very useful when developing microrefrigerators, making it possible to operate on-chip cooling [2,3]. On the other hand, from a different perspective, heat control in circuits is a fundamental requirement, when heat itself is used for operations [4,5], as in the field of coherent caloritronics [6][7][8] that has recently attracted a lot of attention.
In this context, it becomes important to design and analyze devices which provide versatile and tunable control over heat flows. Many such devices have been investigated so far, both theoretically  and experimentally [31][32][33][34][35][36], including thermal transistors [37][38][39][40], valves [41,42], interferometers [43,44] and a large variety of thermal rectifiers. Another key element for heat management, which has been less studied, is represented by circulators: these are multiterminal systems that are able to steer the heat conduction in a preferential direction (e.g., from a given terminal to the next one only in the anticlockwise direction). With this motivation, a three-terminal heat current circulator has recently been suggested [45]. Such a device works similarly as charge current circulators [46][47][48][49], which are crucial for electronics. Beyond that, three-terminal structures for coherent charge current control are also of interest for quantum networks [50,51].
The three-terminal device that has recently been suggested as a heat circulator consists of a ring-like structure with superconducting contacts, penetrated by a magnetic field, see Fig. 1 (b). The superconducting system is of special interest for several reasons. It has been shown in recent years that heat currents carried by quasiparticles in superconducting devices can be coherently controlled via superconducting phase differences [43,52,53].
This has started the field of coherent caloritronics, reviewed in [7], where a circulator element can become part of the basic toolbox. From a fundamental perspective, superconducting phases give an additional control knob and heat can be circulated independently of charge currents. However, at the same time, heat control in general and heat circulation more specifically is equally of interest in normalconducting devices. Still, heat circulation in three-terminal normalconducting systems has, to our knowledge, not been considered so far.
In this paper, we analyze three-terminal heat-current circulators from different perspectives. (1) We investigate the performance of heat circulation in a similar setup as in Ref. [45] with normalconducting as well as with superconducting contacts, see Fig. 1(a). This analysis actually shows that under certain conditions heat circulation can even be more effective in the absence of superconductivity. In contrast to the earlier proposal, where the gap was assumed to be constant, we fully include the self-consistent temperature-dependence of the superconducting gap. (2) The heat circulator of Ref. [45] consists of a simplified setup with three central sites tunnelcoupled to each other and to the superconducting contacts, penetrated by a homogeneous magnetic flux, see Fig. 1(c). We perform a detailed study on conductors deviating from the ideal ring-structure, see e.g. panel (d) and (e). A main consequence of this is that the enclosed flux varies with each trajectory that a particle can follow between the contacts. We show how this can lead to a deterioration and, in the most extreme cases, even to a full suppression of the circulation effect. We also discuss the conditions under which heat circulation is instead preserved. (3) We further investigate how realistic limitations can impact this ideal setting and analyze in detail the statistics of the circulating coefficient and its constituents. Previously, it has been shown how disorder can suppress phase-coherent control of heat flows [54,55]. Here, we analyze the role of sample-to-sample variation of onsite energies and coupling constants in the three-sites structure and in the modified ones in Fig. 1(d-e). In addition, we also consider as the central scattering region a chaotic cavity which we take as a model of an extended quantum system [56]. Here, by using random matrix theory methods [57], we study the sample-to-sample variations of heat conductances and rectifications, in the same spirit of universal conductance fluctuations studies.
This analysis is important for experimental realizations of such circulators. However, most importantly, it gives detailed insights into the working principles of the circulator device. To improve the understanding of the device, we do not only study the circulation coefficient, but also analyze the heat conductances between terminals as well as rectification coefficients between pairs of terminals, separately. This paper is organized as follows. In Sec. II we lay out the general model and define the quantities we use to characterize a three-terminal device with either normal or superconducting contacts. The simplest realization of such a device [ Fig. 1(c)] is investigated in Sec. IV, where normal and superconducting systems are compared. Then, more complicated setups, such as those in Fig. 1(d-e) and also a chaotic scattering region, are studied in Sec. V. Finally, in Sec. VI we draw conclusions. Two Appendices are dedicated to technical and/or complementary details.

II. THREE-TERMINAL CONDUCTORS FOR HEAT CIRCULATION
We consider a general three-terminal device, composed of three contacts, connected to a central conducting region (see Fig. 1). Each terminal has a temperature T i (i = 1, 2, 3) and an electrochemical potential µ i . Our goal is to characterize the heat circulation properties of such a device, in the case where temperature biases (but no voltage ones) are present in the system. These transport properties depend on the details of the central region which is connected to the terminals and, in the framework of a scattering theory approach [58,59], are determined by a scattering matrix. In contrast to Refs. [60,61], we assume the central, ring-shaped device to be small with respect to the quasiparticle coherence length, guaranteeing coherent heat current control over the entire central structure. As Figs. 1(a-b) illustrate, we consider both the case where the terminals are normal metals and the one where they are superconductors. The main aspects distinguishing the two scenarios are that in the superconducting system heat and charge transport are very different and that the superconducting phases yield additional control parameters to tune the heat circulation. However, many of the properties we are interested in also occur in a normalconducting system and this is why we start with this case. In the remainder of this section we introduce the model describing the simplest realization of the scattering region [ Fig. 1(c)].

A. Normalconducting heat circulator
The most basic setup in which heat circulation is possible is a ring that can be modelled by a simple threesites conductor [45]. Any trajectory starting and ending at the same terminal encloses a magnetic flux nΦ which is an integer multiple (n ∈ Z) of the magnetic flux Φ = ring d 2 r B, given by the surface integral of the magnetic field B (perpendicular to the plane) inside the ring. For more complex structures, the enclosed flux depends on the details of the trajectories. This, including the extreme case, where the conductor is a chaotic cavity without any preferred trajectories, will be explored in Sec. V.
The setup represented in Fig. 1(a), together with (c), is a normalconducting version of what has been investigated in Refs. [45,62]. It consists of three sites (labeled by i = 1, 2, 3) with Hamiltonian W , that are connected by a hopping amplitude γ to their respective lead. The Hamiltonian W contains on-site energies i on the diagonal and hopping amplitudest ij (from j to i) elsewhere. In the presence of an external magnetic field the hopping amplitudes have the formt ij = t ij e iθij , where t ij can be taken real and symmetric, while the phases θ ij = − ie i j A · d (e > 0 is the elementary charge) encode the effect of the magnetic field. For any closed path, a particle thus gains a phase shift which is proportional to the flux enclosed by that path. In particular, one has the constraint θ 13 + θ 32 + θ 21 = 2πΦ/Φ 0 ≡ α, with Φ 0 = h/e the flux quantum. Quite generally, for such a system the scattering matrix S(E) at energy E is obtained in terms of W as (see Appendix A for details) (1) where Γ = πγ 2 ν and ν is the density of states in the normal leads. Considering the geometry in Fig. 1(c), we explicitly have Modifications of W , describing the setups in Fig. 1(d-e), are introduced in section V B. Note that the scattering matrix obtained according to Eq. (1) is 2π-periodic in the dimensionless flux α.

B. Superconducting heat circulator
In this subsection, we address the setup sketched in Fig. 1(b), where the three contacts are in a superconducting state. Each of them is therefore characterized by the (real-valued) magnitude of the gap ∆ i and a superconducting phase ϕ i . In order to have a meaningful comparison with a normalconducting system introduced before, we self-consistently take into account the temperature dependence ∆ i (T i ) of the gap. The scattering matrix S sc of this three-terminal superconducting device can be expressed in terms of the scattering matrix S of the central region as [59,63] and its expectation value Ĵ i = J i . Here, we measure all energies with respect to the electrochemical potentials µ i ≡ 0 which we assume to be equal for all i. The annihilation operatorsâ i,n for incoming fluxes are connected to operatorsb i,n for outgoing ones by the (elastic) scattering matrices given in the previous section. The subscript n indicates electron-and hole-like contributions to transport. To keep the notation simple, we write Eq. (7) for the single-channel case; for multi-channel cases, such as treated in Sec. V B 5, the sum over n needs to be extended to count channels in the contacts too. For a normalconducting system, the integral starts at E min = 0, while the lower bound is E min = ∆ in the superconducting case.
Here, ∆ = ∆(T i ) = max i (∆ i (T i )) is the magnitude of the largest gap ∆ i in the problem at given temperatures T i .

B. Heat conductance
To characterize the heat transport of the devices shown in Fig. 1, we consider the linear-response heat conductances in the presence of small temperature differences. A standard linear-response calculation for J i = j κ ij | Ti,j ≡T δT j yields the heat conductance matrix with elements κ ij ≡ κ ij | Ti,j ≡T . For the normalconducting case, one finds, taking together electron-and hole contributions, Here, f (E) = [1 + exp(E/k B T )] −1 is the Fermi function at temperature T . For the system with superconducting contacts, we obtain instead The trace appearing in this case is to be taken over the electron-hole degrees of freedom.

C. Heat rectification and circulation
It has been shown in Ref. [45] that it is possible to achieve a situation where the heat preferentially flows in a given direction by tuning the magnetic flux and/or imposing a superconducting phase bias among the terminals. In linear response, this requires that the heat conductances present a rectification effect with κ ij = κ ji for i = j. To quantify this effect, we introduce the rectification coefficient The heat circulation is quantified via the coefficient It takes the value C = +1 for a perfect counterclockwise circulation and C = −1 for a perfect clockwise one, which arise for the limiting cases, where either κ ij = 0 or κ ji = 0. In the superconducting system, the definitions for R sc and C sc are identical to the previous ones, simply replacing each κ ij with κ sc ij . It is instructive to understand how the circulation and rectification coefficients C and R are related to each other. In the fully symmetric case, where all anticlockwise heat conductances are equal to κ cw and all the clockwise ones are equal to κ cw , the rectification and circulation are given by The relation between the two coefficients is more direct in the case of weak rectification, κ ij − κ ji κ ij + κ ji , where we find

IV. LINEAR-RESPONSE HEAT CIRCULATION IN AN IDEAL THREE-SITES SETUP
In this section, we show the results of our analysis for the ideal three-sites setup of Fig. 1(c), starting with the normalconducting system and then comparing its performance with the superconducting heat circulator.

A. Normalconducting contacts
We consider the specific setup combining Fig. 1(a) with (c), for which the model has been introduced in Sec. II A. Its properties are particularly simple in the symmetric situation with equal on-site energies i = and hopping amplitudes t ij = t; we will address this simple case in the present section. Then, all heat conductances in the same direction (clockwise/anticlockwise) are equal Figure 2. Circulation coefficient C for the basic three-sites ring. In panel (a) as a function of the dimensionless magnetic flux α for a hopping amplitude t = Γ and different temperatures (inset: as a function of temperature for different hopping amplitudes and α = 3π/2). (b) Low-temperature limit, where κij can be evaluated from Eq. (14). Here, we plot again C vs α for different on-site energies (inset: C vs t/Γ for α = 3π/2 and the same values of as in the main plot).
to each other. The circulation coefficient, obtained from Eqs. (8) and (11), is shown as a function of the magnetic flux in Fig. 2(a) and as a function of temperature (inset of (a)). Concerning the dependence on the magnetic flux, the ideal ring shape of the system constrains C to vanish at α = 0, π, 2π, meaning that for these fluxes it is equally probable for the heat current to circulate in both directions. Moreover, the circulation coefficient is antisymmetric about the point α = π and is maximal at α = 3π/2 for a large range of parameters. As far as temperature effects are concerned, Fig. 2(a) shows that by increasing T , the circulation coefficient is typically reduced. The typical scale on which the effect of a finite temperature becomes important is k B T ∼ Γ.
This means in particular, that in the regime k B T Γ the temperature dependence of the circulation coefficient (and of the underlying heat conductances) is negligible. The reason for this is that the scattering matrix S can in this limit be taken as energy-independent, setting E = 0 in Eq. (1). Then, Eq. (8) reduces to where κ 0 = π 2 k 2 B T /(3h) is the thermal conductance quantum. 1 The scattering matrix elements at E = 0, obtained from (1), have the simple form , where¯ = /Γ,t = t/Γ and y = 1 + i¯ . Note, that in this special case, heat and charge circulation are the same, as the thermal and electrical conductances are related by the Wiedemann-Franz law. We now show how C depends on and t, in this lowtemperature regime, k B T Γ. In Fig. 2(b), C is shown as a function of α for different on-site energies and the inset shows the dependence on t. We see that the optimal hopping amplitude is t = Γ but a close-to-maximal circulation coefficient is also found for a quite large range of values around t = Γ. In addition, Fig. 2(b) shows that the circulation coefficient is not very sensitive to : a deviation of the on-site energy from = 0 makes the curves more asymmetric, slightly shifting the value of α at which the maximum circulation coefficient is reached. Overall, as also shown in the inset, the effect of increasing | | is to reduce the circulation coefficient.

B. Superconducting contacts
We now compare the results for the normalconducting device, analyzed above, with those of a superconducting circulator, introduced in Sec. II B. Before starting the discussion, it is worth noting that, in the absence of phase biases, also the superconducting device is completely symmetric and therefore all clockwise heat conductances are equal to each other (idem for the anticlockwise heat conductances). We will first consider precisely this case and show the effect of phase biases later on. Moreover, in order to keep the discussion as simple as possible, we assume from now on that the three superconductors have the same gap amplitude ∆ = ∆ 1 = ∆ 2 = ∆ 3 , the temperature dependence of which is calculated selfconsistently.
We show in Fig. 3(a) C sc as a function of α for = 0, different temperatures and t = Γ = 5k B T c , where T c is the critical temperature of the superconductor. Similarly to the normal system, in this low-temperature regime, it is possible to reach a close-to-maximal clockwise (C sc = −1) or counterclockwise (C sc = +1) circulation. Hence, as for the normal circulator, the most favorable regime for a good circulation coefficient is the low-temperature regime. Again this regime is fixed by k B T Γ, where now T should at the same time not exceed the critical temperature T c for the system to be in the superconducting state. Qualitatively, in this regime we find the same result as in Ref. [45], however at modified parameters and with small changes due to the selfconsistent evaluation of the temperature-dependence of the gap amplitude ∆.
When reducing Γ with respect to T and T c , C sc is typically reduced: an example of such a trend is shown in the inset, where Γ = k B T c . This considerable temperaturedependence is due to the fact that the energy dependence of the scattering matrix is enhanced compared to the one of the normalconducting system by the presence of the superconducting gap.
As further features, we observe that the value of α at which the maximal circulation coefficient |C sc | is reached in general depends on the interplay between the parameters t, T and Γ. In many cases, though, the maximum is found at α = 3π/2, as in Fig. 3(a). In addition, Fig. 3(b) shows that the maximal circulation coefficient  at α = 3π/2 is achieved for a hopping amplitude close to t = Γ (here, again, the exact value of the maximum as a function of t slightly depends on temperature). This plot is specific to the choice Γ/k B T c 1: at lower ratios, the behavior is quite different (not shown).
Comparing Fig. 2 and 3 only reveals small differences between the normal-and superconducting systems. It is natural to ask whether the superconducting device offers some advantages with respect to the normal one. Let us start with the case where no phase biases are imposed among the superconducting terminals and compare the performances of the normal-and superconducting devices, which we present in Fig. 4(a). Here, we show the ratio between the circulation coefficients C sc and C at α = 3π/2, where we have the maximal counterclockwise heat circulation for most values of t. We see that the normal system usually performs better than the superconducting one, at any temperature. 2 In Fig. 4 we chose Γ/k B T c = 5. We have verified that only for much larger ratios (over 30), it is possible that C sc > C for hoppings t ∼ 0.1Γ. In addition, one should be aware that the heat conductance in the superconducting device is typically smaller than that of the normal one, as shown in Fig.  4(b). This means that the amount of circulated heat is typically larger in the normalconducting compared to the superconducting setup. We can therefore conclude that, 2 When T = Tc, the two circulation coefficients are equal, since superconductivity is suppressed. in the absence of phase biases, there is no specific advantage of the superconducting with respect to the normalconducting device.
Nevertheless, the latter is more versatile as it offers an additional control parameter to tune the circulation, namely the possibility of imposing phase biases between different terminals. In Ref. [45] it was shown that when the heat circulation is controlled by just imposing phase biases (and no magnetic flux), an opposite behavior compared to Fig. 3(b) emerges. That is, a higher circulation coefficient is reached at lower hopping amplitudes t. Here, we show the combined effect of phase bias and magnetic field: in Fig. 5 we plot C sc as a function of α and ϕ 3 , fixing the other superconducting phases to ϕ 1 = ϕ 2 = 0. As we can see, the highest circulation coefficient is reached when t is close to Γ, namely when the dependence of C sc on the phase ϕ 3 is quite weak [ Fig. 5(b)]. In contrast, at low hoppings the circulation is more sensitive to variations of the superconducting phase, albeit the overall circulation coefficient is smaller [ Fig. 5(a)].

V. SAMPLE-TO-SAMPLE VARIATIONS
Having described the behavior of the simplest possible setup, in this section, we show how non-ideal operational conditions affect the performance of the heat circulator. Focusing on the low-temperature regime, where the energy dependence of the scattering matrix in κ can be neglected, see Eqs. (8) and (14), we consider two other mechanisms leading to deviations from the ideal condition. A first ingredient is represented by random variations of the parameters of the model (hopping amplitudes and on-site energies), in order to see whether sample-tosample variations limit the usefulness of the device. Secondly, we investigate what happens when the device does not have a ring-like structure, as shown in the modified setups in Fig. 1(d-e). Here, the main difference with the basic model in Fig. 1(c) is that a path starting and ending at the same terminal can include different fractions of the total magnetic flux penetrating the structure. To avoid unnecessary complications, we restrict the analysis to the normalconducting case. This choice is motivated by the fact that superconducting terminals exhibit qualitatively the same behavior, as shown in the previous section.

A. Variations of the hopping and onsite energies
As previously shown, the circulation coefficient of the ideal three-sites system is not particularly sensitive to the variation of the hopping amplitude and on-site energy (see Fig. 2). Therefore, C is also expected to be quite robust to sample-to-sample variations of the model parameters. In Fig. 6, we show that this is indeed the case, a result that was already anticipated in Ref. [45]. We do not restrict ourselves to the symmetric case i = and t ij = t, but we consider the general coupling matrix in Eq. (2). In particular, we have taken t ij , i to be uniformly distributed, i.e. |t ij − t ij | < δ(t ij )/2, with average t ij = Γ and full width δ(t ij ) = 2Γ and | i − i | < δ( i )/2, with i = 0 and δ( i ) = 2Γ. We find that the reduction of the circulation coefficient is roughly only about 15% (clearly by reducing the range of variations of the parameters, the performance of the device is even less affected). Furthermore, Fig. 6 shows that C varies only a little around its ensemble-average C . All these features are peculiar of the simple model of Fig. 1(c); in the following, we show that with increasing complexity of the device also the impact of the sampleto-sample variations grows.

B. Trajectory-dependent enclosed flux
As a next step, we investigate the impact of a modified structure of the central scattering region on the behavior of the heat circulator device. In particular, we consider the possibility that different paths starting and ending at a given terminal enclose a different magnetic flux. First, we consider a modification of the simple toy model that we obtain by adding extra lattice sites to the minimal model, as shown in Figs. 1(d) and (e). Secondly, we consider a more realistic model of an extended central region that we describe by a chaotic cavity modelled by random scattering matrices. As we will see, while single realizations can still be tuned to act as heat circulators, these increasingly complicated systems have a considerable impact on the device performance.  Although the result is valid for any W, we restrict the discussion to nearest neighbor coupling only, as shown in the sketch of Fig. 1(d) and (e). The phases associated with the transition from any site ζ to any of its neighbors ξ are calculated according to θ ξζ = − ie ξ ζ A · d , using the symmetric gauge A = B(−y, x, 0)/2. As for the basic three-sites model, the resulting scattering matrix is periodic in the normalized flux α, although the periodicity is no longer 2π, due to the more complex possible paths that can be followed between any two terminals. It depends on the detailed geometry of the system, as we discuss in the following.

Hexagonal central structure
Let us start with the analysis of the "hexagonal" model in Fig. 1(e). We initially consider equal on-site energies ( i = µ = ), therefore having W 11 = 1 3 in Eq. (16).
Likewise, the other blocks are obtained by considering equal hopping strengths t µi = t (µ = a 1 ) between pairs of external sites and t a1i = t between the central site a 1 and its neighbors (see Fig. 1). Here, thanks to the symmetry of the setup, we have again κ 12 = κ 23 = κ 31 = κ cw and κ 13 = κ 32 = κ 21 = κ cw as for the basic three-sites ring. Moreover, κ cw is obtained from κ cw by simply reversing the magnetic flux; therefore we can just focus on one of them. In Fig. 7 we show the heat conductance κ cw as a function of the magnetic flux, for t = 0.5Γ and various parameters. For a disconnected center site, t = 0, the plot is 2π-periodic, as for the simplest three-sites model, because in this case any allowed path beginning and end-  Fig. 1(e).
ing at the same terminal encloses the whole flux too. However, as soon as t = 0, this is no longer true and different paths are available, thus changing the periodicity. For this geometry, it is 12π because the minimal flux enclosed by a path starting and ending at the same terminal is α/6 instead of α. The other important aspect emerging from the plots is that the behavior of the heat conductance is quite sensitive to the variations of both t and . This is reflected as well by the circulation coefficient, see Fig. 7(c), which depends in a highly nontrivial way on the model parameters. Note however, that despite this strong parameter dependence, in several instances a very high (and even maximal) circulation coefficient C can be obtained. At the same time, this behavior indicates less robustness in the circulation coefficient against sampleto-sample variations of the model parameters, compared to what we have previously illustrated in Fig. 6 for the simple three-sites model.
In Fig. 8, we show the averaged circulation coefficient over 2000 samples generated for a random choice of hoppings and on-site energies. As for the simple three-sites model, the parameters were allowed to vary independently from each other, meaning that the external and internal hopping amplitudes, as well as the on-site energies were not constrained to be equal among each other. In Fig. 8, we consider the average values t µi = Γ (µ = a 1 ), t a1i = 0.4Γ, i = 0.3Γ and full widths δ(t µi ) = 0.6Γ, δ(t a1i ) = 0.8Γ, δ( i ) = 2Γ. The dashed black line shows the circulation coefficient, calculated for a fixed value of the parameters. This indeed shows that variations in the parameters have a much more pronounced impact on the magnitude of the circulation coefficient (compared to Fig. 6). Also, parameter variations in the hexagonal model make the variance (purple line) to be of the same order of magnitude as the average value C . Note however, that the circulation effect is not fully suppressed but persists.

Square central structure
Let us now come to the analysis of the setup with a square structure, as sketched in Fig. 1(d). This configuration shows similar features as the hexagonal one, which we have just illustrated. Only details are different, due to the different symmetry of the system. For instance, the periodicity of the heat conductances of the square structure is 8π, because (as soon as the internal hopping t is non zero) the minimal flux that can be enclosed in a path starting and ending at the same terminal is α/4. Moreover, the position of the terminals is now such that not all heat conductances are equal to each other: κ 23 = κ 31 = κ 12 (and similarly for the counterclockwise direction). However, apart from these specific differences, the qualitative behavior is the same: the system is quite sensitive to a variation of the model parameters, which results in a considerable reduction of the circulation coefficient when sample-to-sample variations are introduced. This is shown in Fig. 9, where we observe a stronger drop in the circulation coefficient, even compared to the one of the hexagonal structure in Fig. 8. This can be attributed to the increased asymmetry of the device with respect to the three contacts.

Three-terminal chaotic cavity
Given the fact that, as we have seen, the performance of the system deteriorates considerably when increasing the complexity of the minimal model only slightly, it is an important question to understand how much circulation can be expected in more realistic models for extended scattering regions. To answer to this question, we study the case of an extended central system (cavity), e.g., a large quantum dot. The dynamics in such a system (with irregular boundaries) is chaotic, see Ref. [56] for a review. Of course, this example is quite far apart from the initial simple ring-shaped system considered in Sec. II A. Nonetheless, we will see that heat circulation is still possible under certain conditions.
We start by investigating the system in the case where only a single mode (M = 1) of each lead is coupled to the central cavity. Such a system can be realized by having a quantum point contact, tuned to the first conductance step, in between the cavity and each of the leads. A successful method to address chaotic and disordered systems is random matrix theory [57]. According to this approach, the scattering matrix S for a system with broken time-reversal symmetry is a random matrix distributed in the Circular Unitary Ensemble (CUE). However, random matrices from the CUE do not carry information on the magnetic field. An extension of the CUE has therefore been developed to include the intensity of the magnetic field B as a parameter, resulting in the formula [66] S(B) = U 11 + U 12 1 N − R(B)U 22 −1 R(B)U 21 . (17) Here, U ij are the four blocks of a (3 + N ) × (3 + N ) random matrix U distributed according to the Circular Orthogonal Ensemble (unitary and symmetric matrices). Furthermore, we have R(B) = exp(BQ), Q being an arbitrary real and antisymmetric matrix. As long as N 1, the detailed choice of the matrix Q has been shown to be irrelevant and the result only depends on the parameter Tr(Q 2 ), related to the Thouless energy E th and the mean level spacing δ in the cavity [66]. Therefore, it is convenient to use for R the parametrization R(x) = exp(xD) [62], D being an antisymmetric matrix with Tr(D 2 ) = −1, while x is a dimensionless quantity related to the magnetic flux piercing the cavity as x ∝ α E Th /δ. The exact proportionality coefficient is a numerical factor of order 1 that depends on the precise shape of the cavity [66]. With this parametrization, the distribution of the scattering matrix S interpolates between the Circular Orthogonal Ensemble (COE) at x = 0 (time-reversal symmetric system) and the CUE at x 1 (when time reversal symmetry is fully broken). In contrast to the previous models, the scattering matrix obtained from Eq. (17) is no longer periodic in α as the magnetic flux enclosed in the path between any two terminals of the system can assume arbitrary values.
The trend indicated in the two setups with extra lattice sites is continued in the chaotic system: in the same way as the variations of the model parameters produced quite different results for both the heat conductance and the circulation coefficient, here any different realization of the random scattering matrix (17) results in a completely different outcome for the same quantities. Some examples illustrating this behavior are shown in Fig. 10(a). As a result, after averaging over a large number of random samples, the circulation coefficient is suppressed to zero. This is shown in Fig. 10(b), where the ensemble-averaged circulation coefficient C is plotted as a function of x. The average is performed over 10000 samples, generated according to Eq. (17), with N = 40. It should be emphasized, however, that the sample-to-sample variations of C are quite large, as shown by the variance (purple curve); this confirms that it is likely that a given random realization produces a good circulation coefficient for some fine-tuned values of the magnetic flux, which depend on the specific device. In particular, Fig. 10(a) shows some realizations where a close-to-optimal circulation coefficient is reached for several values of x ∝ α E Th /δ. Let us now focus on the evolution of Var(C) as a function of x: we observe a transition from a vanishing variance at small x towards a more or less constant value when x increases. This evolution is the result of the progressive breaking of time reversal symmetry. Indeed, at zero magnetic flux, we have C = 0 for every realization of the scattering matrix, since Onsager's reciprocity guarantees that κ ij = κ ji . As a result, also the variance Var(C) vanishes in this case. On the other hand when time-reversal symmetry is broken, it is possible to have Figure 11. (a) Averaged heat conductance κ13 (orange) and its variance (purple). The dashed lines highlight the analytical predictions from Eqs. (18) and (19). (b) Averaged rectification coefficient R13 (orange) and its variance (purple). The dashed line shows the analytical result from Eq. (20). In both panels the averages are taken over 10000 random samples, generated according to Eq. (17), with N = 40. Finally, the shaded bands correspond to 95% (darker) and 99% (lighter) confidence intervals. κ ij = κ ji and therefore the circulation coefficient varies. In order to have an independent check of the magnitude of Var(C) at large x, it is useful to find analytical expressions in limiting cases. We discuss in the following how to achieve this, considering that at large x the scattering matrix is distributed in the CUE.
With this goal in mind, we first investigate the sampleto-sample variations of the heat conductance and the rectification coefficient. They are reported in Fig. 11, showing κ 13 /κ 0 and R 13 , together with their variances. The values of the ensemble-averaged heat conductance can be calculated analytically for x = 0 and x 1. Indeed the ensemble averaging amounts to an integration in the unitary group and, in the case of the heat conductance, Eq. (8) shows that we have to integrate a polynomial function of the scattering matrix S. It is known how to perform such integrations (see for instance [67]) and we get (for i = j) where we used that the size of the scattering matrix S is in our case d = 3. These values are in good agreement with the numerical average in Fig. 11(a) for x large enough. In a similar way, it is possible to obtain analytical results when x = 0 and the scattering matrix is distributed in the COE. In this case we find again in agreement with Fig. 11(a).
Let us now consider the rectification coefficient and the circulation coefficient. At zero magnetic flux, we have already observed that R ij = C = 0 for every realization. We then consider the case of large x and do the averaging over the unitary group. It is easy to conclude on a general basis that R ij = C = 0 thanks to the possibility of relabelling the indices in the definitions (10) and (11) when computing the integration over the unitary group (see App. B 1). Concerning the variance, we take a specific parametrization of U (3) in terms of trigonometric functions [68] in order to compute ensemble averages. For the rectification we obtain (see App. B 1) for any i = j. Notice that the numerical value found in Fig. 11(b) at large x perfectly matches with the above analytical result. Finally, we have checked that the outcome of the integration over U (3) is consistent with the value of Var(C) found in Fig. 10(b), although we are not able to provide an analytic result for this quantity.

Chaotic cavity with multi-mode leads
So far, we have seen that when a single conduction channel connects the reservoirs to the chaotic cavity, the circulation is still highly efficient for some sample-specific values of the magnetic flux. Moreover, we have provided analytic results for the average and variance of heat conductances and rectification coefficients. It is natural to ask whether the circulation effect survives even when there are M channels of each lead connected to the central cavity. As before, the average heat conductances and their variance can be obtained analytically. In the time-reversal-symmetric case (x = 0) we have whereas for broken time-reversal symmetry the result is Notice that for large M there is no difference in the average heat conductance (M/3 in both cases), while a tiny difference persists in the variance (4/81 in the CUE and 5/81 in the COE). Next, we look at what happens to the rectification and the circulation coefficient by increasing M . We directly work in the limit x 1, so that the scattering matrix is distributed in the CUE and perform a numerical simulation by generating random matrices of increasing sizes. The result is shown in Fig. 12, where we see that both Var(R ij ) and Var(C) decay as M −2 for large M . This behavior signals that the different conducting channels are not independent of each other 3 and the coupling among them results in a faster decay of the variances of R ij and C. In addition, the decrease of the variances implies that sample-to-sample variations around the average become smaller and smaller by increasing M , meaning that, even for a single realization, the circulation coefficient is strongly suppressed if many channels are present. More precisely, starting from a finetuned value of the magnetic flux at which a given device has a good performance with one conducting channel, the circulation coefficient will decrease as 1/M by increasing the number of modes in the leads. This indicates that in order to maintain a good performance, the reservoirs have to be connected to the central scattering region via quantum point contacts, in such a way that at most a few conduction channels are open.
Finally, we investigate whether the M −2 behavior observed in the numerical simulation, Fig. 12, can also be understood analytically. In the following we show that indeed the decay of Var(R ij ) and Var(C) can be found exactly in the large-M limit. Starting with the rectification, we can consider the approximation with X = κ ij −κ ji and Y = κ ij +κ ji . This approximation is expected to work well when X/Y does not depart too much from its average, which corresponds to the large-M limit in our case. We clearly have X = 0; moreover Cov(X, Y ) = 0 as can be verified via direct substitution. We have already calculated Y = 2 κ ij = 2M/3 and the only thing left is the evaluation of Var(X), for which we obtain .
The approximation (23) then yields This is exactly the behavior found with the best fit (blue line) in Fig. 12. Notice also that the above formula predicts Var(R ij ) = 3/16 = 0.1875 for M = 1, to be compared with the exact result 3 − 4 ln 2 = 0.2274 [see Eq. (20)]. Concerning the circulation, we can apply the same method, with X, Y = κ 13 κ 32 κ 21 ∓ κ 12 κ 23 κ 31 and evaluate Var(X) exactly for all M . However, the result is cumbersome (see Appendix B 2) and in the large-M limit, we can adopt an alternative strategy, by considering a function g(x 1 , . . . , where each x i is one heat conductance, according to the definition of C. Its variance is estimated as with all derivatives being evaluated at x i = x i . The result of this calculation gives the estimate which again matches with the best fit in Fig. 12 and captures the M −2 decay with the right prefactor.

VI. CONCLUSION
In conclusion, we have presented a detailed analysis of three-terminal conductors with normal-or superconducting contacts acting as heat current circulators. We have shown that the presence of superconducting terminals as considered in a previous proposal [45] is not an essential ingredient, even if they introduce further tunability on the device. Normalconducting systems have a similar (and often even improved) circulation coefficient. The essential requirement is the presence of a magnetic flux which breaks time-reversal symmetry.
Importantly, we have also investigated to what extent non ideal devices affect the circulation coefficient. In slightly modified setups compared to the proposal in [45], introducing the possibility of trajectories enclosing different amounts of magnetic flux, we found a much more important sensitivity on the system parameters. Therefore, the device is less robust with respect to sample-to-sample variations of these parameters, even though it is possible to fine-tune it to obtain high circulation coefficients. Finally, we addressed the extreme case of a chaotic scattering region and investigated the statistics of the heat conductances and the circulation performance. Here, while on average the circulation effect is completely suppressed, specific realizations still exhibit high performances, provided that the number of conducting channels is low.
An interesting issue to be still addressed is to understand the behavior of the heat current correlators (noise) in such multi-terminal devices and what information can be extracted from them. Despite being less studied with respect to its charge counterpart, heat current noise is an interesting topic for the community and is being more and more investigated [69][70][71][72][73]. We leave this issue for future works.

the hopping Hamiltonian
The scattering region is formed by N sites (0 1 , . . . 0 N ), 0 i being connected to the i-th chain, and by M additional sites that are not directly connected to the chains (a 1 , . . . , a M ). The central region is then made of N + M sites, coupled to each other via an (N + M) × (N + M) matrix W. We write the corresponding Hamiltonian as where the block W 11 (W 22 ) of the matrix W describes the subspace of the sites 0 1 , . . . , 0 N (a 1 , . . . , a M ) only. Diagonal elements of these two blocks contain the onsite energies of the N (M) sites. The blocks W 12 and W 21 take cross couplings into account. Finally, the coupling between the scattering region and the leads is The free spectrum of the tight-binding chains is E q = −2v cos q. Now, in order to find the scattering matrix of the system, we consider an incoming wave from the chain j and write the scattering state in the chain i as (for n ≤ −1) ψ ni = δ ij e iqn + S ij e −iqn . (A4) Next, one has to solve the Schrödinger equation ξ ζ|H|ξ ψ ξ = E q ψ ζ , where ξ and ζ can take values n i and a 1 , . . . , a M . For ζ = −1 i we have For ζ = 0 i one finds Finally, for ζ = µ = a 1 , . . . , a M , By eliminating ψ µ and ψ 0i and recalling (A4) one eventually finds the matrix equation −v(e −2iq 1 N +Se 2iq )−γ 2 W 11 + W 12 (E q 1 M − W 22 ) −1 W 21 − E q 1 N −1 (e −iq 1 N +Se iq ) = E q (e −iq 1 N +Se iq ) , (A8) which is solved by (neglecting a global phase factor) where the notation A/B stands for B −1 A. Next, we linearize the spectrum E q ≈ 2v(q − π/2) (which amounts to approximate the density of states ν in the leads as a constant) and consider the wide-band limit, obtaining where W = W 11 + W 12 E1 M − W 22 −1 W 21 and Γ = γ 2 /v = πγ 2 ν, ν being the density of states in the leads. Finally, in the case where just the N sites connected to the chains are present (and no additional ones), only the block W 11 exists and the scattering matrix is found by simply letting W = W 11 . In this way, and also taking N = 3, one recovers Eq. (1). |U ji | 2 /(|U ij | 2 + |U ji | 2 ). It is clear that f 1 (U ) = f 2 (U 0 U ), U 0 being the unitary matrix which swaps rows and columns (i, j) in the matrix U . Therefore By the very same reasoning, one concludes that C = 0.
We now evaluate Var(R ij ) = R 2 ij and derive Eq. (20). First of all, following the same argument given above, one easily shows that Var(R ij ) is the same for every i = j. Next, to get Eq. (20), we consider a single channel for each lead and thus the integration is performed over the unitary group U (3). We use the parametrization of Ref. [68], according to which the group measure can be written as dφ i d(cos 4 θ 1 )d(cos 2 θ 2 )d(cos 2 θ 3 ) , (B2) with 0 ≤ θ 1 , θ 2 , θ 3 ≤ π/2 and 0 ≤ φ i ≤ 2π. The parametrization of the elements U ij is given in Eq. (2.10) of Ref. [68]. As we said, Var(R ij ) does not depend on i and j and then we can take R 13 which is the most convenient for the calculation. By applying to the matrix U the unitary transformation that exchanges the first two columns one gets Var(R 13 ) = dV |U 13 | 2 − |U 32 | 2 |U 13 | 2 + |U 32 | 2 2 . (B3) The parametrization for the elements entering the last expression is [68] U 13 = cos θ 1 sin θ 2 e iφ4 and U 32 = cos θ 1 sin θ 3 e iφ5 . By using these relations in the previous formula, together with Eq. (B2), one finds Var(R 13 ) =  for M ≥ 2. As we can expect, this is a poor estimate for small M , but it captures exactly the large-M behaviour 3M −2 /2 and already for M = 8 the error with respect to the true value is less than 5%.