Classical and quantum harmonic mean-field models coupled intensively and extensively with external baths

We study the nonequilibrium steady-state of a fully-coupled network of $N$ quantum harmonic oscillators, interacting with two thermal reservoirs. Given the long-range nature of the couplings, we consider two setups: one in which the number of particles coupled to the baths is fixed (intensive coupling) and one in which it is proportional to the size $N$ (extensive coupling). In both cases, we compute analytically the heat fluxes and the kinetic temperature distributions using the nonequilibrium Green's function approach, both in the classical and quantum regimes. In the large $N$ limit, we derive the asymptotic expressions of both quantities as a function of $N$ and the temperature difference between the baths. We discuss a peculiar feature of the model, namely that the bulk temperature vanishes in the thermodynamic limit, due to a decoupling of the dynamics of the inner part of the system from the baths. At variance with usual cases, this implies that the steady state depends on the initial state of the particles in the bulk. We also show that quantum effects are relevant only below a characteristic temperature that vanishes as $1/N$. In the quantum low-temperature regime the energy flux is proportional to the universal quantum of thermal conductance.


Introduction
In recent years we witnessed a growing interest in the out-of-equilibrium properties of classical and quantum many-body systems, motivated by the remarkable progresses in the research on nanosized structures and cold atoms or ions physics. From the fundamental point of view, this research theme is part of the aims of nonequilibrium statistical mechanics to justify transport properties from microscopic interactions. A key quantity to understand the transport properties of a system is the thermal conductivity κ, which is defined by Fourier's law: where ∇T is the temperature gradient and J is the heat flux flowing through the system. For normal (diffusive) systems thermal transport follows the above law, that is, κ in an intensive quantity that does not depend on the size of the system. At the microscopic level, it is well known that the harmonic approximation of interparticle forces yields a violation of Eq. (1): transport is ballistic since quasi-particles travel undisturbed throughout the system bulk. This was demonstrated by Rieder, Lebowitz and Lieb [1] in their seminal study of the non-equilibrium steady state of the simplest example: the harmonic chain with nearest-neighbor coupling in contact with two external reservoirs. They found out that the thermal conductivity scales as κ ∼ N , where N is the number of particles in the chain, and the temperature profile in the bulk of the system is flat, whereas Fourier's law would yield a linear profile. Many later works considered variants of the harmonic crystal, since in this case one has to deal with linear problems, one can exploit a variety of established techniques like the approach based on transmission coefficients or non-equilibrium Green's functions, see the dedicated sections in the reviews [2,3] and also [4] for a more recent account. Quantum harmonic lattices outside equilibrium have also been studied in various contexts, starting from the simplest case of the one-dimensional chain [5], see e.g. [6] and [7][8][9].
Besides this, it is also established that Fourier's law may be violated in lower dimensions (d = 1, 2) also in presence of nonlinear interactions. This is the phenomenon of anomalous (superdiffusive) heat transport that has been thoroughly studied in the last two decades [2,3,10,11]. In one dimension, the thermal conductivity may scale as κ ∼ N α . There are strong evidences (both numerical and theoretical) that the exponent α can be used to categorize different universality classes [11]. In those cases, Eq. (1) should be replaced by its fractional version, yielding nonlinear temperature profiles across the system [12,13]. It is now understood that one and two-dimensional, non-linear, nonintegrable, short-range models conserving energy, momentum and stretch, superdiffusive transport occurs generically. Admittedly, there exist also instances where diffusive regimes are unexpectedly found numerically [14]. If one conservation law is broken, as in the case of coupled rotors [15,16], standard diffusion is restored. Finally, nonlinear integrable models like for instance the Toda chain [17] should generically display ballistic transport mediated by solitons (see [13,18] and references therein for some recent results). However, there may be exceptions to this intuition, as in the cases of hard-point particles [19] and the classically integrable Landau-Lifshitz spin chain [20].
One can then wonder how, and if, does the picture change for long-range interacting systems, that is, systems in which the inter-particle interaction is of the form V (r) ∼ r −d−σ , in dimension d. There is a vast literature regarding the physics of these systems, see for example the reviews [21] and [22], for classical and quantum systems, respectively. At equilibrium they are characterized by non-additivity and critical exponents depending on the value σ. In particular the critical exponents for σ < 0 are the mean-field ones obtained by putting σ = −d, where d is the dimensionality of the system. At nonequilibrium, the dynamics of long-range systems presents metastable states whose lifetime scales as N [23][24][25] and even lack of thermalization upon interaction with a single external bath [26].
Moreover, since perturbations may propagate with infinite velocities [27,28] anomalous transport properties are expected. Indeed, numerical studies of energy transport in classical long-range interacting systems support this expectation [29][30][31][32][33][34] The long-range version of the XY model [30,33] and the Fermi-Pasta-Ulam-Tsingou chain were analyzed in [31,33,35,36], and violations of Fourier's law were observed depending on the exponent σ of the interaction. On the quantum side, transport in a mean-field harmonic model was studied in [37]; we also refer once again to [22] and references therein for further results. From an analytical perspective, a classical model with stochastic momentum exchange was studied in [38], while a hydrodynamic approach for transport in long-range quantum magnets was developed in [39] and [40].
In the present work, we study one of simplest cases of long-range interacting nonequilibrium systems: a network of mean-field coupled harmonic oscillators in contact with two external heat reservoirs, modeled as large ensembles of bosonic oscillators at different temperatures. Models with flat interaction have a paradigmatic importance in the study of long-range systems due to the fact that they may allow for an analytic solution. Moreover, they may be seen as providing the limit σ → −d of systems with spatial couplings of the form ∼ 1/r d+σ , with r being the inter-constituents distance, and a generally interesting question is whether this limit is well defined or singular. Furthermore, spin models with flat interactions have been experimentally realized in trapped-ions systems [22] and for cold atoms in double wells described by the Lipkin-Meshkov-Glick model [41], with the stationary state of the Hamiltonian Mean Field model simulated with cold atoms in a cavity [42]. Recent progress in such systems allows for the study of dynamical and transport properties, despite the effective inherent absence of spatial distance [22].
The goal of this paper is two-fold. On one hand, we want to have an analytically treatable model for which one can study both classical and quantum regimes. On the other hand, in presence of long-range interactions a natural question is how to couple the system with the baths. In fact, while in short-range systems the coupling concerns only the boundary particles, in the case of long-range interactions, one could ask what is the effect of the interaction between the bulk of the system and the bulk of the baths. Since within our models we can study both intensive and extensive (clarified below) couplings of the system with the external baths, they provide a useful playground to investigate such question. To be specific, we will consider these baths to be Ohmic and coupled to the system in two different ways: 1) We attach each baths to a finite subset of particles, for instance just one as it is usually done for short-range systems. We will refer to this case as intensive coupling. This case also includes the situation in which a finite number of sites is coupled to the two baths.
2) Since the system is long-range, it also makes sense to ask what happens if we couple an extensive fraction of the sites of the system to the baths, since the interaction between the system and the environment may well be itself long-ranged. Therefore, we will connect both baths to a macroscopic number (i.e., scaling with N ) of sites.
We will refer to this case as extensive coupling.
In Sec. 2 we introduce the model, and the tecniques we are going to use to compute the flux and the temperature profile. In Sec.3, we study the classical case with intensive couplings, and we discuss the peculiarities of the temperature distribution, while in Sec. 4 we tackle the quantum case. In Sec. 5 and 6 we turn to the case of extensive couplings, in the classical and quantum regime, respectively. Finally we review our results in the different cases, discuss a physical interpretation of our results and draw our conclusions.

The model and the coupling to external baths
We consider a Hamiltonian system describing a network of N fully-connected harmonic oscillators with displacements x i . The Hamiltonian of the model is: where p i = mẋ i , k is the coupling constant and N is the Kac factor, introduced to make the energy extensive. The Hamiltonian can also be cast in the following form, which we will use throughout the paper: where the matrix Φ is given by: We connect the system to two external heat reservoirs at different temperatures, T L and T R . The subscripts here would refer to "left" and "right" bath, but of course, in the mean-field model there is no notion of spatial ordering. Nonetheless, in order to fix the notation we will keep this nomenclature throughout. In the more general situation in which a subset of N L sites is coupled to the left bath, at temperature T L , while another subset of N R sites is coupled to the right bath, at temperature T R . We will also assume that these two baths are two ensembles of bosonic harmonic oscillators, linearly coupled to the system, as explained in [6]. To study the transport properties of this system we will use the Green's function method, which we will briefly illustrate below (for a full account see [6]). By integrating out the oscillators, we get a set of quantum Langevin equation of motion. Since we are interested in the properties of the stationary state, we then switch to Fourier space. The final result is: where ξ(ω), the noises introduced by the baths, and Γ L,R (ω), the Fourier transforms of the memory kernels, are related by the fluctuation-dissipation relation: and an equivalent relation holds for the right bath. The specific form of the Γ matrices depends on the spectral density of the baths. From now on, we will consider the case of Ohmic baths, whose memory kernel is purely imaginary: Also, for simplicity of notation, we conventionally assume to couple to the left-hand bath the sites i = 1...N L and to the right-hand bath the sites i = N − N R + 1...N , Eq. (7) defines the intensive coupling case, whereas (8) refers to the extensive coupling case when both N L and N R are taken to be proportional to N . By plugging (7) in (6) we recover the familiar fluctuation-dissipation relation, for example the classical limit of (6) becomes: The solution of the equations of motion (5) can be written as: where the matrix G ij (ω) is the Green's function, given by: Once G ij is computed, we can use the noise correlation (6) to compute all the correlators of interest. In the following we will report the relevant formulas focusing on the computation of the Green's function (11). Before proceeding further, let us mention that another way to compute the correlators is to explicitly solve the Fokker-Planck equation in the stationary regime. Since the system is quadratic, this amounts to compute the covariance matrix of the canonical variables.
As explained in [1], the covariances satisfy a Lyapunov equation, which can be solved numerically [43]. Note that this approach can be used also for the case of a long-range interaction with σ > −1 (see below). We anticipate that in the mean-field case we tested the results of the calculations obtained in the large N limit against the numerical solutions finding an excellent agreement.

Heat flux
In the classical case the flux is given by: where ∆T = T L − T R . Plugging (7) in (12) we get: Now we need to compute the Green's function by inverting the following matrix: In order to do it we can employ the Sherman-Morrison formula [44], that reads as follows. Given a matrix M , with known inverse M −1 , and two vectors u and v, the inverse of A = M + uv T is: It is easy to show that the Z int matrix (14) can be written as −Z int = D int + uu T , with: Using (15) we can then compute the Green's function (11) exactly. In particular we can compute G 1N and plug the result in (12) to obtain the heat flux: where the function I 1 is given by the following integral (y = ω 2k/m): where we introduced the dimensionless coupling constant: Consider now the second factor of the denominator of the integrand in (19), as a polynomial in s = y 2 : As N → ∞ (21) has a vanishing root: where o(x) indicates a quantity that goes to zero faster than x. By decomposing the integrand in (19) in partial fractions, it is easy to see that the dominant contribution for large N is the one coming from (22): and so the heat flux at leading order in N −1 reads as: We can compare this analytical prediction with the numerical calculation of the integral I 1 as a function of k 1 for several values of N , which is plotted in fig. 1a. It is clear that the two prediction do not match since I 1 has a maximum for some optimal value of the coupling and then goes to zero for large k 1 ‡. To better understand the problem, we also plot N I 1 (k 1 , N ) as a function of k 1 in figures 1a and 1b for different values of N . As N grows, the region of agreement between the predicted scaling with N and dependence on k 1 (given by 23) grows as well. This fact can be understood in the following way: as N grows the maximum of I 1 moves to the right on the k 1 axis, as can be seen from figure 1 so when N → ∞ the maximum is virtually at k 1 = ∞ so that only the linear region of I 1 (k 1 , N ) is visible.

Temperature profile
We consider the kinetic definition of the temperature of the i th site: The classical velocity-velocity correlator can be expressed in terms of the Green's function as [6]: ‡ It is worth noting that this dependence on the coupling constant is qualitatively the same as the one in short-range case [1]. and by substituting in (7) we get: For the first site, i = 1, in the large N limit we get: where the term O(N −1 ) is proportional to the heat flux. The same formula also holds for i = N , with L → R. For all the other sites, i = 1, N , the velocity-velocity correlator is given by: with: This integral can be decomposed in partial fractions and dealt with in the same way as I 1 : it turns out that at leading order for large N we have I 2 = I 1 /k 1 . Therefore the temperature of the generic i-th site is given by: This result may appear at first glance unphysical. The temperature of bulk oscillators do not equilibrate to the average of the temperatures of the baths, but rather vanishes in the thermodynamic limit. This is a peculiarity of the model and depends on the choice of the initial conditions. In the next paragraph we give an explanation of this fact based on the analysis of the equation of motion.

Analysis of the equations of motion in the stationary regime
Let's start from the equations of motion in the time domain of the system coupled to the baths: (we set m = 1 for simplicity) where Φ is defined in (4), and ξ L,R are Gaussian noises with correlation given by (9). We now introduce the "total magnetization" M (t) = i x i /N , and S = x 1 + x N . The equations of motion (31) can then be cast in the following form: Switching to Fourier space, we find the following solution for M (ω) and the position x i (ω) of the uncoupled sites: Notice how the pole on the proper frequency of the system ω 2 = 2k does not give any contribution, as if the baths were unable to properly interact with the system. It is convenient to recast Eq. (36) as: that has no pole on the dispersion law ω 2 = 2k. The mean square velocity of the i th site then reads as: which reproduces exactly formula (28).
To get a better understanding of the physics of the model, let us introduce the relative coordinates z i = x i+1 − x i . Then, equation (34) entails that: so the relative coordinates of the uncoupled particles follow a harmonic motion without being influenced by the baths. This, in turn, means that the initial conditions of the system are essential to determine the properties of the stationary state at long times. Indeed, to solve the equations of motion we should use the Laplace-rather than the Fourier-transform. The use of the latter, made in the previous sections, implicitly assumes x i (0) =ẋ i (0) = 0 for i = 2, ..., N − 2. All our results are thus valid, provided we make this assumption on the initial conditions. To support the above considerations in fig. 2 we report the kinetic temperatures as measured in a Langevin simulation of the equation of motion for two different initial conditions. In the case x i (0) =ẋ i (0) = 0 the results coincide with the result (30). On the other hand, for random initial data the temperatures of the particles not connected to the baths remain at their starting value and do not thermalize at all. We remark that the crucial point is the cancellation of the pole in the dispersion relation of the system, which stems from two properties of the model. The first one is the linearity of the system, that allows the equations of motion to be solved exactly, in terms of the Green's function, which in this analysis is given by Q(ω)/ω. The second one is the conservation of the total magnetization M in absence of external baths, which stems from the mean-field nature of the system: from a mathematical point of view, this is related to the (N − 1)-fold degeneracy of the spectrum of the matrix Φ. In order to demonstrate that the breaking of the matrix Φ may suffice to restore thermalization, we considered two variants of the model. The first one is the quadratic chain with a power-law decaying interaction. In the second one we add a nearest-neighbors coupling term to (4). The Hamiltonians corresponding to these two choices are, respectively: Note that in the case of Hamiltonian (40)  does not appear to be straightforward. For these reasons, we decided, in both cases, to solve numerically the Lyapunov equation for the covariance matrix of the models, following the approach of [1]. The results for the temperature profile are plotted in fig. 3. As we can see, in both cases the profile is flat and given by the average temperature of the baths. Indeed, for the temperature profile to vanish, it suffices that the degeneracy of the matrix Φ is of order N . For example, if we add a pinning potential to the first and last site, the degeneracy of Φ is N − 3, that is still of order N . One can then examine the equations of motion of the system as we did above: the result is that once again we have a number of degrees of freedom of order N completely decoupled from the baths. Thus we guess that in order to have thermalization one needs to break the degeneracy of Φ to a quantity of order 1.
Another possibility is to add nonlinear forces. We performed simulations with the same method used to obtain the data plotted in fig.2, i.e. via the numerical solution of the Langevin equations of motion, now adding a term −x 3 i to the right-hand side of Eqs.31 [45]. We found that kinetic temperatures settle to the average of the temperatures of the baths. Furthermore, the long-range version of the Fermi-Pasta-Ulam-Tsingou chain was numerically studied in [33]: the results for the mean-field case show that the system does thermalize to the average temperature of the baths also in this case.  Finally, we note that one can repeat the analysis presented in this subsection even if the baths induce a coloured noise: once again we obtain that the pole on the proper frequency of the system vanishes.

Heat flux
The heat flux in the quantum case is given by [6]: Furthermore, we will work in the linear response regime: To get the heat flux for intensive couplings we plug (7) into (42) and we expand to first order in ∆T : The Green's function element G 1N is the same as in the classical case, and therefore the heat flux can be written as: where we introduced a dimensionless temperature θ as: and the function I 3 (k 1 , θ, N ) is given by: The integral I 3 (47) is plotted as a function of θ in fig. 4: as expected, the heat flux goes to zero at low temperatures (when θ is small), and saturates at high temperatures (when θ is large). From the figure it is clear that there exist a characteristic temperature scale, which we will call T N (k 1 ), that discriminates between the quantum and classical regimes. It turns out that T N (k 1 ) goes as 1/N , as we can see by a direct computation of I 3 in the large N limit. This computation is reported in the Appendix, and the final result is: where the function g(x) is given by: being Γ(z) the Euler Gamma function. The temperature T N (k 1 ) is the intrinsic temperature scale of the system, below which quantum effects are important. It is given by: To get a better picture of the crossover from the quantum to the classical regime, we consider the ratio between the quantum (48) and the classical heat flux (24), J int q /J int cl . For large N , this ratio is given by the function g(x) defined in (49). Its low and high temperature behaviors can be worked out explicitly are given by: where we used the asymptotic formulas for the digamma function. We can see that at high temperature the quantum flux correctly converges to the classical one, while at low temperature it vanishes linearly with T . In fig. 5 we plot the aforementioned ratio as a function of θ for several values of N : as N increases, the saturation to the classical value takes place at lower values of T . As a final remark, we note that the flux (48) for low temperature is: Remarkably, the quantity among parentheses is recognized to be the quantum of thermal conductance, introduced in [46] for heat transport in ballistic quantum wires. It is a universal quantity, independent of all the system parameters (the coupling constants

Temperature profile
In the quantum case, the velocity-velocity correlator is given by [6]: plugging (7) in (52) we get: For site number 1, that is, the one coupled to the left bath, the leading order term is: We get the same expression for site N , with L → R. The integral in Eq. (54) is logarithmically divergent at large frequencies. This problem is unrelated with the longrange properties of the interaction, it is present even if we couple a single oscillator to two ohmic baths (see for example [47]). The divergence stems from the implicit hypothesis that the bath is able to excite arbitrarily high frequencies, hidden in the choice Γ ∼ γω.
Physically, there has to be a cutoff at high frequencies. Now let's consider the case i = 1, N : in this case it is easy to see that in the linear response regime the term proportional to ∆T vanishes, and the leading term is of order ∆T 0 : where I 4 is given by: While we cannot solve the integral I 4 exactly, we can compute its low-temperature behaviour by approximating the hyperbolic cotangent in (56) with the following series, valid for large x: The integral of the first term can be computed exactly (note that the denominator in (56) is the same as in (19)): where R is the real part of the two roots that do not vanish in the large N limit (which is the same for both of them since they are complex conjugates). The integral of the second term in (57) can be computed fairly easily with some algebra in the large N and low-temperature limit (T T N ). In these limits, we find the following result for the temperature profile: The first term, which is non-zero even if the temperature of the external baths vanishes can be interpreted as a zero point energy of the quantum system.

Extensive coupling, classical case
We now couple the baths to a number of sites that scales as N : the left-hand bath will be coupled to a subset N L = α L N of sites starting, while the right-hand bath will be coupled to a subset of N R = α R N sites. We will also assume the condition N L + N R ≤ N , that is, a site can be coupled at most with one bath. The computation of the Green's function can be carried out in the same way as in the intensive case, the only difference being that now the Γ matrices are given by (8). The matrix that we need to invert is now: As in the intensive case, we can decompose this matrix as −Z ext = D est + uu T , where u is defined in (17), and D ext is given by: Therefore the Green's function can be computed exactly also in the case of extensive coupling to the baths. It is important to note that the results obtained in this section are formally valid also when we couple a finite number of sites, that is, a number that does not scale with N , to the baths (a case that has to be considered an intensive coupling). An analysis similar to the one of the previous section has to be performed in order to extract the proper dependence on N . The result is that the scalings do not change and only the prefactors are affected.

Heat flux
To compute the heat flux, we first substitute (8) in (12) and then we plug in the relevant matrix elements of the Green's function: where we defined I 5 as: Note that the sum over the coupled sites collapses to N R N L |G il | 2 since due to symmetry (and as can be checked by explicit calculation), G il with i = l is actually independent on i and l. As a check, we can recover the results (19) of the intensive case by putting α R = α L = 1/N in (63).
In fig. 7 we report the plot of I 5 as a function of k 1 for some fixed values of α L and α R . The qualitative behaviour is the same as in the intensive case: the heat flux vanishes for both small and strong coupling, and as the fraction of coupled sites decreases, the heat flux decreases as well, as could be expected. It is also interesting to note that in the extensive case the flux does not depend on N , in contrast with the N −1 scaling of the intensive flux (24). However, as we are going to see in the next section, the temperature profile in the bulk still goes to zero as N −1 in the thermodynamic limit for the same reasons as in the intensive case. This seems an inconsistent result, but it is actually only an apparent dichotomy that can be reconciled with the following argument. As we saw in section 3 the coupling with the baths is very weak for the particles in the bulk, so we can picture heat transport as heat flowing, at leading order in N , directly between the sites that are coupled to the baths. This means that if we increase the energy pumped into the system by a factor N -as we do by coupling an extensive number of sites to the baths-the heat flux will increase by that factor, but the temperature profile will still scale as N −1 .

Temperature profile
The velocity-velocity correlator computed via (26) by substituting (8): As in the intensive case, we get different results if i is directly coupled to a bath or not. If i is coupled to the left/right bath we get at leading order T ext cl,i = T L/R , as in the intensive case. If i is not coupled to any bath, then we have: where we introduced the integral I 6 : if we set α L = α R = 1/N we recover the intensive case (29), apart from a factor of 2N due to the different definitions of these integrals. In fig. 8 we report the dependence of I 6 on α L + α R : for (α L + α R ) → 0 I 6 diverges, as it is necessary to match with the intensive case. Also note that even in this case the temperature of the uncoupled sites scales as N −1 : indeed, the same argument used in section 3 for the intensive coupling case holds also in the extensive case, by replacing S with the sum of the positions of the sites coupled to the baths, as can be explicitly checked by solving the equations of motion.

Heat flux
The heat flux is obtained by substituting (8) in (42). In the linear response regime we get, as in the intensive case, a factor related to the derivative of the Bose function: where we introduced the function I 7 : In fig. 9a we plot I 7 as a function of θ with fixed α L , α R : as expected, the heat flux vanishes at low temperature, while it saturates at high temperature. Unfortunately, we cannot compute I 7 exactly as we did with its intensive counterpart I 3 (47), but we can obtain an estimate for the low-temperature behaviour using the following result: so that I 7 for small θ is given by: and the heat flux vanishes linearly with the temperature, as in the intensive case (48). In fig. 9b, we report the numerical exact plot of I 7 and the low-temperature approximation (70) for several values of k 1 , respectively. It is evident that the value of θ below which the linear approximation is valid decreases as a function of k 1 . This fact implies the presence of a characteristic temperature scale of the system, of which we are however unable to provide an explicit expression, since we are unable to solve (68) exactly. As in the intensive case (4.1), at low temperature the flux can be expressed in terms of the quantum of thermal conductance: where now we also have a "geometrical factor" that depends on the fraction of coupled sites. Note that for α L = α R we recover the intensive result (4.1).

Temperature profile
In the quantum case we have to plug (8) in (52): where θ L,R = 2k B T L,R / 2k/m the integral I 8 (θ, k 1 , α L , α R ) as: In the linear response regime we can Taylor expand (73) around ∆T = 0 for T L,R = T ± ∆T /2: where: Coupling Quantity Classical Quantum  (3). We only report the low-temperature behaviour of the quantum results.
In fig. 10a we report the plot of I and we notice that it goes to a nonzero constant at low temperature: indeed, I In fig. 10b we report the plot of I and (78) and we see that there is good agreement. Note that (75) entails that, if we couple the same fraction of sites to the left and to the right bath, the term linear in ∆T vanishes. This was the case for (59), the intensive counterpart of (75), and can be obtained from (75) setting α L = α R = 1/N . Moreover, since one can check numerically that I (0) 8 goes as θ 2 for small θ, we conclude that at low temperature:

Conclusions
In this work we analyzed a harmonic mean-field model in various settings. To summarize our results we refer the reader to the table 1, in which we report the scaling of the temperature profile and the heat flux. We considered both the case in which only two sites are coupled to the baths (which we call intensive coupling case), and the one in which an extensive number of sites is coupled to the baths (which we call extensive coupling case).
Let us now comment on our results, starting from the intensive case. In the classical regime, the peculiar scaling of the temperature profile with N is the result of two properties of the model: having a quadratic Hamiltonian and a degenerate matrix of interactions Φ. If we remove either of these properties, the profile flattens on the average temperature T as we show in figure 3. At low temperatures quantum effects become relevant: interestingly enough, since the scale of temperature of the system scales as T N (k 1 ) ∼ N −1 , the region where quantum effects are noticeable shrinks to a point in the thermodynamic limit. At variance, the classical heat flux scales as N −1 in the classical regime (in contrast with the short-range case, where the flux is constant for large N ). Let us now turn to the extensive coupling case starting from the classical regime. The temperature profile scales once again as N −1 , but the flux is independent of N . Indeed, we are pumping more energy into the system, and so the heat flux is larger. The fact that the temperature profile still scales as N −1 is due to the fact that, since all sites are coupled irrespectively of their distance, heat can simply flow from a site coupled to the left bath to one coupled to the right bath. This minimizes, in the large N limit, the amount of energy given to the uncoupled sites. We also note that at low temperatures, both in the intensive and extensive coupling case, the heat flux vanishes linearly with T . The prefactor is given by the quantum of thermal conductance in the intensive case, as expected from [46], and in the extensive case we get a contribution related to the fraction of coupled sites. It would be interesting to further study the quantum regime in the extensive coupling case to better understand the dependence of the temperature scale of the system with respect to the coupling constant to the baths and the fractions of coupled sites. From the analysis conducted in this paper, one concludes that for a mean-field system the coupling to an external bath essentially affects only the sites directly coupled to the bath. It would be interesting to see if and to what extent this property stays true if the role of the bath is played by a subsystem that we trace out, for example, in the computation of the entanglement entropy of the system.
Note added: During the completion of this manuscript, an interesting and related paper by L. Defaveri, C. Olivares and C. Anteneodo [48] appeared in the arXiv. The authors study heat transport in the same model, for the classical case with intensive couplings, while we also considered the quantum case and the one with extensive coupling. Our results and conclusions are in perfect agreement with theirs in the classical case with intensive coupling and equal masses. They extend their analysis to the case of graded and random masses, where the degeneracy of the model is removed and the system reaches the thermal state. In our paper we show that this also happens if one breaks the degeneracy by adding a power-law long-range interaction or a nearest-neighbor one.

Re(z)
Im(z) Figure A1: The contour Γ n used to compute the integral I 3 where we conveniently made the change of variable x = y/k 2 and a is given by: As in the classical case (19), we cannot directly take the limit N → ∞, because in this limit a = 0 and (A.1) diverges. To compute (A.1) we employ contour integration and Cauchy's theorem. Let us introduce the following function of complex variable z: The poles of f are all located on the imaginary axis, at the following positions (with the corresponding residue): z n = inπ, n ∈ Z{0}, Res n = 2iπa 2 n (a 2 − π 2 n 2 ) 2 , (A.4) z ±a = ±ia, Res ±a = ∓ ia 2 sin 2 (a) . (A.5) Consider now the contour Γ n plotted in figure A1: it is composed by a segment [−R, R] and a semicircle C n of radius R, which is such that Γ n contains the first n z k poles and the one in z a . Let now be I R the integral of f (z) over the aformentioned segment. Then, by the residue theorem, we have: Res k + Res +a . (A.6) I 3 can then be obtained by taking the limit R → ∞ of I R as follows: We now have to compute the limit R → ∞, which also entails the limit n → ∞ of the right-hand side of (A.6). For large R the integral of f over C n is: Cn dzf (z) ≈ iR π 0 e iθ dθ sinh 2 (Re iθ ) = iR(−2i coth(R)/R) → 2. (A.8) The sum over the residues becomes a series that can be resummed. We can thus finally express I 3 as: where the function g(x) is given by (49).