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

Francesco Andreucci1,2,∗, Stefano Lepri3,4,†, Stefano Ruffo1,2,3,‡, Andrea
Trombettoni5,1,2,6,# 1SISSA - Via Bonomea 265, 34136 Trieste, Italy  2INFN
sezione di Trieste, 34100 Trieste, Italy  3Consiglio Nazionale delle Ricerche,
Istituto dei Sistemi Complessi, Via Madonna del Piano 10, I-50019 Sesto
Fiorentino, Italy  4Istituto Nazionale di Fisica Nucleare, Sezione di Firenze,
Via G. Sansone 1 I-50019, Sesto Fiorentino, Italy  5Department of Physics,
University of Trieste, Strada Costiera 11, I-34151 Trieste, Italy  6CNR-IOM
DEMOCRITOS Simulation Center, Via Bonomea 265, I-34136 Trieste, Italy

###### Abstract

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 case, this implies that the steady state depends
on the initial state of bulk particles. We also show that quantum effect 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.

Email:
<EMAIL_ADDRESS>
<EMAIL_ADDRESS>
<EMAIL_ADDRESS>
<EMAIL_ADDRESS>

##  1 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 $\kappa$, which
is defined by Fourier’s law:

$\mathcal{J}=-\kappa\nabla T,$ (1)

where $\nabla T$ is the temperature gradient and $\mathcal{J}$ is the heat
flux flowing through the system. For normal (diffusive) systems thermal
transport follows the above law, that is, $\kappa$ 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
inter-particle 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 $\kappa\sim 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 $\kappa\sim N^{\alpha}$. There are
strong evidences (both numerical and theoretical) that the exponent $\alpha$
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, non-integrable, 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)\sim r^{-d-\sigma}$, 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 $\sigma$. In particular the critical exponents for $\sigma<0$ are
the mean-field ones obtained by putting $\sigma=-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
$\sigma$ 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
non-equilibrium 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 $\sigma\to-d$ of systems with spatial couplings of the
form $\sim 1/r^{d+\sigma}$, 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.

##  2 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:

$H=\frac{1}{2m}\sum_{i=1}^{N}p_{i}^{2}+\frac{k}{2N}\sum_{ij}(x_{i}-x_{j})^{2},$
(2)

where $p_{i}=m\dot{x}_{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:

$H=\frac{1}{2m}\sum_{i=1}^{N}p_{i}^{2}+\frac{1}{2}\sum_{i,j=1}^{N}\Phi_{ij}x_{i}x_{j},$
(3)

where the matrix $\Phi$ is given by:

$\Phi_{ij}=2k\left(\delta_{ij}-\frac{1}{N}\right).$ (4)

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:

$\displaystyle-m\omega^{2}x_{i}(\omega)$
$\displaystyle=-\Phi_{ij}x_{j}(\omega)+(\Gamma_{L,ij}(\omega)+\Gamma_{R,ij}(\omega))x_{j}(\omega)+\xi_{L,i}(\omega)+\xi_{R,i}(\omega),$
(5)

where $\xi(\omega)$, the noises introduced by the baths, and
$\Gamma_{L,R}(\omega)$, the Fourier transforms of the memory kernels, are
related by the fluctuation-dissipation relation:

$\Braket{\xi_{L,i}(\omega)\xi_{L,j}(\omega^{\prime})}=\text{Im}[\Gamma_{L,ij}(\omega)](1+f(\omega,T_{L,}))\frac{\hbar}{\pi}\delta(\omega+\omega^{\prime}),$
(6)

and an equivalent relation holds for the right bath. The specific form of the
$\Gamma$ 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$,

$\displaystyle\Gamma_{L,ij}(\omega)=i\gamma\omega\delta_{ij}\delta_{i1},\
\quad$ $\displaystyle\Gamma_{R,ij}=i\gamma\omega\delta_{ij}\delta_{1N},$ (7)
$\displaystyle\Gamma_{L,ij}(\omega)=i\gamma\omega\delta_{ij}\sum_{k=1}^{N_{L}}\delta_{jk},\
\quad$
$\displaystyle\Gamma_{R,ij}=i\gamma\omega\delta_{ij}\sum_{k=N-N_{R}+1}^{N}\delta_{jk},$
(8)

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:

$\Braket{\xi_{L/R,i}(\omega)\xi_{L/R,j}(\omega^{\prime})}=\frac{\gamma
k_{B}T_{L/R}}{\pi}\delta(\omega+\omega^{\prime})\delta_{ij}.$ (9)

The solution of the equations of motion (5) can be written as:

$x_{i}(\omega)=\sum_{j}G_{ij}(\omega)(\xi_{L,j}+\xi_{R,j}),$ (10)

where the matrix $G_{ij}(\omega)$ is the Green’s function, given by:

$G_{ij}(\omega)=\left[-m\omega^{2}\mathbb{I}+\Phi+\Gamma_{L}(\omega)+\Gamma_{R}(\omega)\right]_{ij}^{-1}\equiv
Z^{-1}(\omega).$ (11)

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 $\sigma>-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.

##  3 Intensive coupling, classical case

###  3.1 Heat flux

In the classical case the flux is given by:

$\mathcal{J}_{cl}=\frac{k_{B}\Delta
T}{\pi}\int_{-\infty}^{\infty}d\omega\,\Tr\left[G(\omega)\Gamma_{L}G^{\dagger}(\omega)\Gamma_{R}\right],$
(12)

where $\Delta T=T_{L}-T_{R}$. Plugging (7) in (12) we get:

$\mathcal{J}_{cl}^{int}=\frac{k_{B}\Delta
T}{\pi}\gamma^{2}\int_{-\infty}^{\infty}d\omega\,\omega^{2}|G_{1N}(\omega)|^{2}.$
(13)

Now we need to compute the Green’s function by inverting the following matrix:

$\displaystyle
Z^{int}_{ij}=\begin{cases}-m\omega^{2}-i\gamma\omega+2k(1-1/N),\quad&i=j=1,N\\\
-m\omega^{2}+2k(1-1/N),\quad&i=j\neq 1,N\\\ -2k/N,\quad&i\neq j\end{cases}.$
(14)

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
$\bm{u}$ and $\bm{v}$, the inverse of $A=M+\bm{u}\bm{v}^{T}$ is:

$A^{-1}=M^{-1}+\frac{A^{-1}\bm{u}\bm{v}^{\,T}A^{-1}}{1+\bm{v}^{\,T}A^{-1}\bm{u}}.$
(15)

It is easy to show that the $Z^{int}$ matrix (14) can be written as
$-Z^{int}=D^{int}+\bm{u}\bm{u}^{T}$, with:

$\displaystyle
D^{int}_{ij}=\begin{cases}m\omega^{2}+i\gamma\omega-2k,\quad&i=j=1,N,\\\
m\omega^{2}-2k,\quad&i=j\neq 1,N,\\\ 0,\quad&i\neq j,\end{cases},$ (16)
$\displaystyle u_{i}=\sqrt{2k/N}\,\,\forall i.$ (17)

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:

$\mathcal{J}_{cl}^{int}=\frac{k_{B}\Delta T\sqrt{2k/m}}{2\pi}I_{1}(k_{1},N),$
(18)

where the function $I_{1}$ is given by the following integral
$(y=\omega\sqrt{2k/m})$:

$I_{1}(k_{1},N)=\frac{2k_{1}^{2}}{N^{2}}\int_{-\infty}^{\infty}\frac{(y^{2}-1)^{2}}{(y^{2}-1)^{2}+k_{1}^{2}y^{2}}\frac{dy}{y^{2}(y^{2}-1)^{2}+k_{1}^{2}(y^{2}-\frac{2}{N})^{2}},$
(19)

where we introduced the dimensionless coupling constant:

$k_{1}^{2}=\frac{\gamma^{2}}{2mk}.$ (20)

Consider now the second factor of the denominator of the integrand in (19), as
a polynomial in $s=y^{2}$:

$s(s-1)^{2}+k_{1}^{2}(s-2/N)^{2}.$ (21)

As $N\rightarrow\infty$ (21) has a vanishing root:

$s_{0}(k_{1},N)=-4k_{1}^{2}/N^{2}+o(N^{-2}),$ (22)

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):

$I_{1}=\frac{2k_{1}^{2}}{N^{2}}\int_{-\infty}^{\infty}\frac{dy}{y^{2}+4k_{1}^{2}/N^{2}}=\frac{\pi
k_{1}}{N},$ (23)

and so the heat flux at leading order in $N^{-1}$ reads as:

$\mathcal{J}_{cl}^{int}=k_{B}\Delta T\sqrt{2k/m}\frac{k_{1}}{2N}.$ (24)

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. 1$a$. 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}$111It is worth noting that this dependence on
the coupling constant is qualitatively the same as the one in short-range case
[1].. To better understand the problem, we also plot $NI_{1}(k_{1},N)$ as a
function of $k_{1}$ in figures 1$a$ and 1$b$ 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\rightarrow\infty$
the maximum is virtually at $k_{1}=\infty$ so that only the linear region of
$I_{1}(k_{1},N)$ is visible.

Figure 1: In fig. $(a)$ we report the plot of $I_{1}(k_{1},N)$ for
$N=100,200,300,400,500$ from top to bottom, respectively. Notice how the flux
reaches a maximum value and then decreases. In fig. $(b)$ we report the plot
of $NI_{1}(k_{1},N)$ for $N=1000,2000,3000,4000,5000$, respectively. Note how
the agreement with the analytical prediction (represented by the dashed black
line) gets better and better as $N$ grows.

###  3.2 Temperature profile

We consider the kinetic definition of the temperature of the $i^{th}$ site:

$T_{i}=\frac{m}{k_{B}}\Braket{\dot{x}_{i}^{2}}.$ (25)

The classical velocity-velocity correlator can be expressed in terms of the
Green’s function as [6]:

$\Braket{\dot{x}_{i}^{2}}=\frac{k_{B}T_{L}}{\pi}\int_{-\infty}^{\infty}d\omega\omega\left[G(\omega)\Gamma_{L}(\omega)G^{\dagger}(\omega)\right]_{ii}+(L\rightarrow
R),$ (26)

and by substituting in (7) we get:

$\displaystyle\Braket{\dot{x}_{i}^{2}}$ $\displaystyle=\frac{\gamma
k_{B}T_{L}}{\pi}\int_{-\infty}^{\infty}d\omega\omega^{2}|G_{i1}(\omega)|^{2}+\frac{\gamma
k_{B}T_{R}}{\pi}\int_{-\infty}^{\infty}d\omega\omega^{2}|G_{iN}(\omega)|^{2}.$
(27)

For the first site, $i=1$, in the large $N$ limit we get:

$\displaystyle\Braket{\dot{x}_{1}^{2}}$ $\displaystyle=\frac{\gamma
k_{B}T_{L}}{\pi}\int_{-\infty}^{\infty}\frac{\omega^{2}}{(m\omega^{2}-2k)^{2}+\gamma^{2}\omega^{2}}+O(N^{-1})=\frac{k_{B}}{m}T_{L}+O(N^{-1}),$

where the term $O(N^{-1})$ is proportional to the heat flux. The same formula
also holds for $i=N$, with $L\rightarrow R$. For all the other sites, $i\neq
1,N$, the velocity-velocity correlator is given by:

$\displaystyle\Braket{\dot{x}_{i}^{2}}=\frac{k_{B}}{m}\frac{T_{L}+T_{R}}{2\pi}I_{2}(k_{1},k_{2},N),$
(28)

with:

$\displaystyle
I_{2}(k_{1},N)=\frac{4k_{1}}{N^{2}}\int_{0}^{+\infty}\frac{dy}{y^{2}(y^{2}-1)^{2}+k_{1}^{2}(y^{2}-\frac{2}{N})^{2}}.$
(29)

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:

$T_{cl,i}^{int}=\frac{T_{L}+T_{R}}{2\pi}\frac{I_{1}}{k_{1}}=\frac{T_{L}+T_{R}}{2N}.$
(30)

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.

###  3.3 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)

$\ddot{x}_{i}=-\sum_{j}\Phi_{ij}x_{j}+\delta_{i1}(\xi_{L}-\gamma\dot{x}_{i})+\delta_{i1}(\xi_{L}-\gamma\dot{x}_{i}),$
(31)

where $\Phi$ is defined in (4), and $\xi_{L,R}$ are Gaussian noises with
correlation given by (9). We now introduce the “total magnetization”
$M(t)=\sum_{i}x_{i}/N$, and $S=x_{1}+x_{N}$. The equations of motion (31) can
then be cast in the following form:

$\displaystyle\ddot{S}=-\lambda S-2kS+4kM+\xi,$ (32)
$\displaystyle\ddot{M}=(\xi-\lambda\dot{S})/N,$ (33)
$\displaystyle\ddot{x}_{i}=-2kx_{i}+2kM,\quad i=2,...,N-1.$ (34)

Switching to Fourier space, we find the following solution for $M(\omega)$ and
the position $x_{i}(\omega)$ of the uncoupled sites:

$\displaystyle
M(\omega)=\frac{2k-\omega^{2}}{N\omega^{2}(\omega^{2}+i\lambda\omega-2k)-4ik\lambda\omega}\xi(\omega),$
(35) $\displaystyle x_{i}(\omega)=\frac{2kM(\omega)}{2k-\omega^{2}}.$ (36)

Notice how the pole on the proper frequency of the system $\omega^{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:

$x_{i}(\omega)=\frac{Q(\omega)}{\omega}\xi(\omega),\quad
Q(\omega)\equiv\frac{-2k}{N}\frac{1}{\omega(\omega^{2}-2k)+i\lambda(\omega^{2}-4k/N)},$
(37)

that has no pole on the dispersion law $\omega^{2}=2k$. The mean square
velocity of the $i^{th}$ site then reads as:

$\Braket{x^{2}_{i}}=\gamma\frac{T_{L}+T_{R}}{\pi}\int d\omega|Q(\omega)|^{2},$
(38)

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:

$\ddot{z}_{i}=-2kz_{i},\quad i=2,...,N-2,$ (39)

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)=\dot{x}_{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)=\dot{x}_{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.

Figure 2: Temperature temperature obtained from simulations of the Langevin
equation of motion for $N=32$, intensive baths with $T_{L}=1.5,T_{R}=0.5$,
$\gamma=0.5$. Upper panel: zero initial conditions $x_{i}(0)=\dot{x}_{i}(0)=0$
corresponding to the choice adopted in the analytical calculations. Lower
panels: random initial conditions where $x_{i}(0)$ and $\dot{x}_{i}(0)$ are
drawn from a Gaussian distribution with zero average and variance
$(T_{L}+T_{R})/2$. Averages are over trajectories of $10^{5}$ time units.

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(\omega)/\omega$.
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 $\Phi$. In order to demonstrate that the
breaking of the matrix $\Phi$ 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:

$\displaystyle
H_{1}=\sum_{i}\frac{p_{i}^{2}}{2}+\frac{k}{2N_{\sigma}}\sum_{ij=1}^{N}\frac{(x_{i}-x_{j})^{2}}{|i-j|^{1+\sigma}},\quad
N_{\sigma}=\sum_{l=1}^{N}l^{-1-\sigma}$ (40) $\displaystyle
H_{2}=\sum_{i}\frac{p_{i}^{2}}{2}+\frac{k}{2N}\sum_{ij=1}^{N}(x_{i}-x_{j})^{2}+\frac{g}{2}\sum_{i}(x_{i+1}-x_{i})^{2}.$
(41)

Note that in the case of Hamiltonian (40) the expression of the Green’s
function is not known in the literature, to the best of our knowledge. On the
other hand, for the Hamiltonian (41), while it is possible to extend the
previous analysis for $g\neq 0$, the calculation 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 $\Phi$ is of order
$N$. For example, if we add a pinning potential to the first and last site,
the degeneracy of $\Phi$ 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 $\Phi$ 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_{i}^{3}$ 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.

Figure 3: Numerical result for the temperature profile $T_{cl,i}^{int}$ for
the Hamiltonians (40) and (41) on the top panel and on the bottom panel,
respectively. We used the following values for the parameters:$T_{L}=1.5$,
$T_{R}=0.5$ $k=g=1$, $\sigma=-0.5$.

##  4 Intensive coupling, quantum case

###  4.1 Heat flux

The heat flux in the quantum case is given by [6]:

$\mathcal{J}_{q}=\int_{-\infty}^{\infty}d\omega\,\Tr\left[G(\omega)\Gamma_{L}(\omega)G^{\dagger}(\omega)\Gamma_{R}(\omega)\right]\frac{\hbar\omega}{\pi}\left[f(\omega,T_{L})-f(\omega,T_{R})\right].$
(42)

Furthermore, we will work in the linear response regime:

$T_{L}-T_{R}=\Delta T\ll(T_{R}+T_{L})/2=T.$ (43)

To get the heat flux for intensive couplings we plug (7) into (42) and we
expand to first order in $\Delta T$:

$\mathcal{J}^{int}_{q}=\frac{\hbar\gamma^{2}\Delta
T}{\pi}\int_{-\infty}^{+\infty}d\omega\,\omega^{3}|G_{1N}(\omega)|^{2}\frac{\partial
f(\omega,T)}{\partial T}.$ (44)

The Green’s function element $G_{1N}$ is the same as in the classical case,
and therefore the heat flux can be written as:

$\mathcal{J}^{int}_{q}=\frac{k_{B}\Delta
T\sqrt{2k/m}}{4}I_{3}(k_{1},\theta;N),$ (45)

where we introduced a dimensionless temperature $\theta$ as:

$\theta=\frac{2k_{B}T}{\hbar}\sqrt{m/2k},$ (46)

and the function $I_{3}(k_{1},\theta,N)$ is given by:

$I_{3}=\frac{4}{\pi}\frac{k_{1}^{2}}{\theta^{2}N^{2}}\int_{-\infty}^{\infty}dy\,\frac{(y^{2}-1)^{2}}{(y^{2}-1)^{2}+k_{1}^{2}y^{2}}\frac{y^{2}/\sinh^{2}(y/\theta)}{y^{2}(y^{2}-1)^{2}+k_{1}^{2}(y^{2}-2/N)^{2}}$
(47)

Figure 4: In fig. $(a)$ we report $I_{3}$ as a function of $\theta$, with
$k_{1}=3$, and $N=500,1000,2000$ from top to bottom, respectively. In fig.
$(b)$ we plot the function $g(x)$ defined in (49).

The integral $I_{3}$ (47) is plotted as a function of $\theta$ in fig. 4: as
expected, the heat flux goes to zero at low temperatures (when $\theta$ is
small), and saturates at high temperatures (when $\theta$ 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:

$\mathcal{J}_{q}^{int}=\frac{k_{B}\Delta
T}{2N}k_{1}\sqrt{\frac{2k}{m}}\,g\biggr{(}\frac{T_{N}(k_{1})}{T}\biggr{)},$
(48)

where the function $g(x)$ is given by:

$\displaystyle
g(x)=\frac{x^{2}}{\pi^{2}}\left[\psi^{(1)}\left(1+\frac{x}{\pi}\right)-\psi^{(1)}\left(1-\frac{x}{\pi}\right)\right]+\frac{x^{2}}{\sin^{2}(x)}-2x,$
(49) $\displaystyle\psi^{(1)}(z)=\frac{d^{2}}{dz^{2}}\Gamma(z),$

being $\Gamma(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:

$T_{N}(k_{1})=\frac{k_{1}}{N}\frac{\hbar}{k_{B}}\sqrt{\frac{2k}{m}}.$ (50)

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), $\mathcal{J}_{q}^{int}/\mathcal{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:

$g\biggr{(}\frac{T_{N}}{T}\biggr{)}=\begin{cases}1,\quad&T\gg T_{N},\\\ \\\
\dfrac{\pi}{3}\dfrac{T}{T_{N}}\sim TN,\quad&T\ll T_{N},\end{cases}$

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 $\theta$ 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:

$\mathcal{J}_{q}^{int}=\left(\frac{\pi^{2}k_{B}^{2}T}{3h}\right)\Delta T.$
(51)

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 in our case).

Figure 5: We report the plot of the ratio
$\mathcal{J}_{q}^{int}/\mathcal{J}^{int}_{cl}$ as a function of $\theta$ for
$k_{1}=3$ and $N=500,1000,2000$ from bottom to top, respectively.

###  4.2 Temperature profile

In the quantum case, the velocity-velocity correlator is given by [6]:

$\Braket{\dot{x}_{i}^{2}}=\int_{-\infty}^{\infty}d\omega\,\omega^{2}\left[\left(G^{+}(\omega)\Gamma_{L}(\omega)G^{+\dagger}(\omega)\right)_{ii}\frac{\hbar\omega}{2}\coth\left(\frac{\hbar\omega}{2k_{B}T_{L}}\right)+(L\rightarrow
R)\right],$ (52)

plugging (7) in (52) we get:

$\Braket{\dot{x}_{i}^{2}}=\frac{\gamma}{\pi}\int_{-\infty}^{\infty}d\omega\omega^{2}\left[|G_{i1}|^{2}\frac{\hbar\omega}{2}\coth\left(\frac{\hbar\omega}{2k_{B}T_{L}}\right)+|G_{iN}|^{2}\frac{\hbar\omega}{2}\coth\left(\frac{\hbar\omega}{2k_{B}T_{L}}\right)\right],$
(53)

For site number $1$, that is, the one coupled to the left bath, the leading
order term is:

$\Braket{\dot{x}^{2}_{1}}=\frac{\hbar\gamma}{2\pi}\int_{-\infty}^{\infty}\frac{d\omega}{(m\omega^{2}-2k)^{2}+\gamma^{2}\omega^{2}}\omega^{3}\coth\left(\frac{\hbar\omega}{2k_{B}T_{L}}\right).$
(54)

We get the same expression for site $N$, with $L\rightarrow R$. The integral
in Eq. (54) is logarithmically divergent at large frequencies. This problem is
unrelated with the long-range 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
$\Gamma\sim\gamma\omega$. Physically, there has to be a cutoff at high
frequencies.

Now let’s consider the case $i\neq 1,N$: in this case it is easy to see that
in the linear response regime the term proportional to $\Delta T$ vanishes,
and the leading term is of order $\Delta T^{0}$:

$\Braket{\dot{x}_{i}^{2}}=\frac{\hbar\gamma}{2m^{2}}I_{4},$ (55)

where $I_{4}$ is given by:

$I_{4}=\frac{1}{N^{2}}\int_{0}^{\infty}dy\frac{y\coth(y/\theta)}{y^{2}(y^{2}-1)^{2}+k_{1}^{2}(y^{2}-2/N)^{2}}.$
(56)

Figure 6: Plot of $I_{4}$ as a function of $\theta$ for $k_{1}=3$ and
$N=1000$. In the inset we plot $I_{4}$ at low temperatures.

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$:

$\coth(x)=1+2\sum_{k=0}^{\infty}e^{-2(k+1)x}.$ (57)

The integral of the first term can be computed exactly (note that the
denominator in (56) is the same as in (19)):

$I_{4}(\theta=0)=\frac{1}{2N^{2}}\left[2\ln\left(\frac{N}{2k_{1}}\right)+\frac{R}{\sqrt{1-R^{2}}}\left(\frac{\pi}{2}+\arctan\left(\frac{R}{\sqrt{1-R^{2}}}\right)\right)\right],$
(58)

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\ll
T_{N}$). In these limits, we find the following result for the temperature
profile:

$T_{q,i}^{int}=\frac{\hbar\gamma
k_{B}}{2m}\left[I_{4}(\theta=0)+\frac{\pi^{2}}{48}\frac{\theta^{2}}{k_{1}^{2}}\right].$
(59)

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.

##  5 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}=\alpha_{L}N$ of sites starting, while
the right-hand bath will be coupled to a subset of $N_{R}=\alpha_{R}N$ sites.
We will also assume the condition $N_{L}+N_{R}\leq 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 $\Gamma$ matrices are given by (8). The matrix that we need
to invert is now:

$\displaystyle
Z^{ext}_{ij}=\begin{cases}-m\omega^{2}-i\gamma\omega+2k(1-1/N),\quad&i=j=1...N_{L},N-N_{R}+1...N,\\\
-m\omega^{2}+2k(1-1/N),\quad&i=j,\,\,\text{otherwise},\\\ -2k/N,\quad&i\neq
j.\end{cases}$ (60)

As in the intensive case, we can decompose this matrix as
$-Z^{ext}=D^{est}+uu^{T}$, where $\bm{u}$ is defined in (17), and $D^{ext}$ is
given by:

$\displaystyle
D^{ext}_{ij}=\begin{cases}m\omega^{2}+i\gamma\omega-2k,\quad&i=j=1...N_{L},N-N_{R}+1...N,\\\
m\omega^{2}-2k,\quad&i=j,\text{otherwise},\\\ 0,\quad&i\neq j.\end{cases}$
(61)

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.

###  5.1 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:

$\mathcal{J}_{cl}^{ext}=\frac{k_{B}\Delta
T}{\pi}\int_{-\infty}^{+\infty}d\omega(\gamma\omega)^{2}\sum_{i=N-N_{R}+1}^{N}\sum_{l=1}^{N_{L}}|G_{il}|^{2}=\frac{k_{B}\Delta
T\sqrt{2k/m}}{2\pi}I_{5}(k_{1},\alpha_{L},\alpha_{R}),$ (62)

where we defined $I_{5}$ as:

$I_{5}=k_{1}^{2}(2\alpha_{L}\alpha_{R})\int_{-\infty}^{\infty}dy\frac{(y^{2}-1)^{2}}{\left[(y^{2}-1)^{2}+k_{1}^{2}y^{2}\right]\left[y^{2}(y^{2}-1)^{2}+k_{1}^{2}(y^{2}-(\alpha_{R}+\alpha_{L}))^{2}\right]}.$
(63)

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\neq l$ is actually independent on $i$ and $l$. As a check, we
can recover the results (19) of the intensive case by putting
$\alpha_{R}=\alpha_{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 $\alpha_{L}$ and $\alpha_{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}$.

Figure 7: Plot of $I_{5}$ as a function of $k_{1}$ with
$\alpha_{L}=1/6,\alpha_{R}=1/10$ and $\alpha_{L}=\alpha_{R}=1/2$. The first
choice of parameters corresponds to the curve with larger maximum.

###  5.2 Temperature profile

The velocity-velocity correlator computed via (26) by substituting (8):

$\Braket{\dot{x}_{i}^{2}}=\frac{k_{B}\gamma}{\pi}\left[T_{L}\int_{-\infty}^{\infty}d\omega\omega\sum_{i=1}^{N_{L}}|G_{ik}(\omega)|^{2}+T_{R}\int_{-\infty}^{\infty}d\omega\omega\sum_{i=N-N_{R}+1}^{N}|G_{ik}(\omega)|^{2}\right].$
(64)

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_{cl,i}^{ext}=T_{L/R}$, as in the intensive case. If $i$ is not
coupled to any bath, then we have:

$T_{cl,i}^{ext}=\frac{\alpha_{L}T_{L}+\alpha_{R}T_{R}}{\pi
N}I_{6}(k_{1},\alpha_{L},\alpha_{R}),$ (65)

where we introduced the integral $I_{6}$:

$I_{6}=k_{1}\int_{-\infty}^{\infty}dy\frac{1}{\left[(y^{2}-1)^{2}+k_{1}^{2}y^{2}\right]\left[y^{2}(y^{2}-1)^{2}+k_{1}^{2}(y^{2}-\alpha_{L}-\alpha_{R})^{2}\right]},$
(66)

if we set $\alpha_{L}=\alpha_{R}=1/N$ we recover the intensive case (29),
apart from a factor of $2N$ due to the different definitions of these
integrals.

Figure 8: The plot of $I_{6}$ as a function of $\alpha_{L}+\alpha_{R}$ for
$k_{1}=3$.

In fig. 8 we report the dependence of $I_{6}$ on $\alpha_{L}+\alpha_{R}$: for
$(\alpha_{L}+\alpha_{R})\rightarrow 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.

##  6 Extensive coupling, quantum case

###  6.1 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:

$\mathcal{J}_{q}^{ext}=\frac{k_{B}\Delta
T}{\pi}\int_{-\infty}^{+\infty}d\omega(\gamma\omega)^{2}\sum_{i=N-N_{R}+1}^{N}\sum_{l=1}^{N_{L}}|G_{il}|^{2}\frac{\partial
f}{\partial T}=\frac{k_{B}\Delta
T\sqrt{2k/m}}{\pi}I_{7}(k_{1},\theta,\alpha_{L},\alpha_{R}),$ (67)

where we introduced the function $I_{7}$:

$I_{7}=\alpha_{L}\alpha_{R}\frac{k_{1}^{2}}{\theta^{2}}\int_{-\infty}^{\infty}dy\frac{\left(y^{2}-1\right)^{2}}{\left(y^{2}-1\right)^{2}+k_{1}^{2}y^{2}}\frac{y^{2}/\sinh(y/\theta)}{y^{2}\left(y^{2}-1\right)^{2}+k_{1}^{2}\left(y^{2}-(\alpha_{L}+\alpha_{R})\right)^{2}}.$
(68)

Figure 9:  In fig. $(a)$ we plot $I_{7}$ (68) as a function of $\theta$ for
$k_{1}=3$, $\alpha_{L}=1/20,1/3,1/2$ and $\alpha_{R}=1/20,1/5,1/2$ from bottom
to top, respectively. In fig. $(b)$ we plot $I_{7}$ for low $\theta$ with
$\alpha_{L}=\alpha_{R}=1/3$ and $k_{1}=5,10,15$ from bottom to top,
respectively. The black dashed line is the linear approximation (70).

In fig. 9$a$ we plot $I_{7}$ as a function of $\theta$ with fixed
$\alpha_{L},\alpha_{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:

$\lim_{\theta\rightarrow
0}\frac{3}{\pi^{2}\theta^{3}}\frac{y^{2}}{\sinh^{2}(y/\theta)}=\delta(y),$
(69)

so that $I_{7}$ for small $\theta$ is given by:

$I_{7}=\frac{\pi^{2}}{3}\frac{\alpha_{R}\alpha_{L}}{(\alpha_{R}+\alpha_{L})^{2}}\theta,$
(70)

and the heat flux vanishes linearly with the temperature, as in the intensive
case (48). In fig. 9$b$, 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 $\theta$ 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 (51), at low temperature the flux
can be expressed in terms of the quantum of thermal conductance:

$\mathcal{J}_{q}^{ext}=\left(\frac{\pi^{2}k_{B}^{2}T}{3h}\right)\frac{4\alpha_{R}\alpha_{L}}{(\alpha_{R}+\alpha_{L})^{2}}\Delta
T,$ (71)

where now we also have a “geometrical factor” that depends on the fraction of
coupled sites. Note that for $\alpha_{L}=\alpha_{R}$ we recover the intensive
result (51).

###  6.2 Temperature profile

In the quantum case we have to plug (8) in (52):

$\displaystyle\Braket{\dot{x}_{i}^{2}}=\gamma\biggr{[}$
$\displaystyle\int_{-\infty}^{\infty}d\omega\omega\sum_{i=k}^{N_{L}}|G_{ik}(\omega)|^{2}\frac{\hbar\omega}{\pi}\coth\left(\frac{\hbar\omega}{2k_{B}T_{L}}\right)$
$\displaystyle+\int_{-\infty}^{\infty}d\omega\omega\sum_{i=N-N_{R}+1}^{N}|G_{ik}(\omega)|^{2}\frac{\hbar\omega}{\pi}\coth\left(\frac{\hbar\omega}{2k_{B}T_{R}}\right)\biggr{]}.$
(72)

For the sites directly coupled to the baths , i.e. $i=1,...,N_{L}$ and
$i=N-N_{R}+1,...,N$, the correlator diverges due to the contribution of the
unphysical high frequencies. On the other hand, if $i=N_{L}+1,...N-N_{R}$,
then we get:

$T_{q,i}^{ext}=\frac{\hbar\gamma}{2\pi
Nmk_{B}}\left[\alpha_{L}I_{8}(\theta^{L},k_{1},\alpha_{L,R})+\alpha_{R}I_{8}(\theta^{R},k_{1},\alpha_{L},\alpha_{R})\right],$
(73)

where $\theta^{L,R}=2k_{B}T_{L,R}/\hbar\sqrt{2k/m}$ the integral
$I_{8}(\theta,k_{1},\alpha_{L},\alpha_{R})$ as:

$I_{8}=\int_{-\infty}^{\infty}dy\frac{y\coth(y/\theta)}{y^{2}(y^{2}-1)^{2}+k_{1}^{2}(y^{2}-\alpha_{L}-\alpha_{R})^{2}}.$
(74)

In the linear response regime we can Taylor expand (73) around $\Delta T=0$
for $T_{L,R}=T\pm\Delta T/2$:

$T_{q,i}^{ext}=\frac{\hbar\gamma}{2\pi
Nmk_{B}}\left[(\alpha_{L}+\alpha_{R})I_{8}^{(0)}+\frac{k_{B}\Delta
T}{\hbar\sqrt{2k/m}}(\alpha_{L}-\alpha_{R})I_{8}^{(1)}\right],$ (75)

where:

$\displaystyle
I_{8}^{(0)}=\int_{-\infty}^{\infty}dy\frac{(y^{2}-1)^{2}y\coth(y/\theta)}{\left[y^{2}(y^{2}-1)^{2}+k_{1}^{2}(y^{2}-\alpha_{R}-\alpha_{L})^{2}\right]},$
(76) $\displaystyle
I_{8}^{(1)}=\theta^{-2}\int_{-\infty}^{\infty}dy\frac{y^{2}/\sinh^{2}(y/\theta)}{\left[y^{2}(y^{2}-1)^{2}+k_{1}^{2}(y^{2}-\alpha_{R}-\alpha_{L})^{2}\right]}.$
(77)

In fig. 10$a$ we report the plot of $I^{(0)}_{8}$ and we notice that it goes
to a nonzero constant at low temperature: indeed, $I_{8}^{(0)}(\theta=0)$ can
be interpreted as the contribution to the temperature of the zero-point energy
of the particles. We can get an analytical estimate of $I_{8}^{(1)}$ at low
temperature using once again (69):

$I_{8}^{(1)}=\frac{\pi^{2}}{3k_{1}^{2}(\alpha_{L}+\alpha_{R})^{2}}\theta.$
(78)

In fig. 10$b$ we report the plot of $I_{8}^{(1)}$ 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
$\Delta T$ vanishes. This was the case for (59), the intensive counterpart of
(75), and can be obtained from (75) setting $\alpha_{L}=\alpha_{R}=1/N$.
Moreover, since one can check numerically that $I_{8}^{(0)}$ goes as
$\theta^{2}$ for small $\theta$, we conclude that at low temperature:

$T_{q,i}^{ext}=\frac{\hbar\gamma}{2\pi
Nmk_{B}}\left[(\alpha_{L}+\alpha_{R})I_{8}^{(0)}(\theta=0)+\frac{k_{B}\Delta
T}{\hbar\sqrt{2k/m}}\frac{\pi^{2}}{3k_{1}^{2}}\frac{(\alpha_{L}-\alpha_{R})}{(\alpha_{L}+\alpha_{R})^{2}}\theta\right].$
(79)

Figure 10: In fig. $(a)$ we report the plot of $I_{8}^{(0)}$ as a function of
$\theta$ with $k_{1}=3$ and $\alpha_{L}+\alpha_{R}=0.4,0.8$ from respectively.
In fig. $(b)$ we report the plot of $I_{8}^{(1)}$ as a function of $\theta$
with $k_{1}=3$ and $\alpha_{L}+\alpha_{R}=0.4$. The black dashed line is the
linear approximation (78)

##  7 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).

Coupling | Quantity | Classical | Quantum
---|---|---|---
Intensive | $\mathcal{J}^{int}$ | $\mathcal{J}_{cl}^{int}\sim\Delta T/N$ | $\mathcal{J}_{q}^{int}\sim T\Delta T$
$T_{i}^{int}$ | $T_{cl,i}^{int}\sim T/N$ | $T_{q,i}^{int}\sim T^{2}+O(\ln N/N^{2})$
Extensive | $\mathcal{J}^{ext}$ | $\mathcal{J}_{cl}^{ext}\sim\Delta TN^{0}$ | $\mathcal{J}_{q}^{ext}\sim T\Delta T$
$T_{i}^{ext}$ | $T_{cl,i}^{ext}\sim N^{-1}$ | $T_{q,i}^{ext}\sim const/N+T/N$
Table 1: Summary of main results about thermal transport in the mean-field
harmonic model (3). We only report the low-temperature behaviour of the
quantum results.

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 $\Phi$. 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})\sim 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.

SL acknowledges support from the program Collaborations of excellence in
research and education granted by SISSA (Trieste, Italy) where this work has
been initiated. This work is part of the MIUR-PRIN2017 project Coarse-grained
description for non-equilibrium systems and transport phenomena (CO-NEST) No.
201798CZL.

## Appendix: Calculation of $I_{3}$

Since the denominator of the integrand in (47) is the same as the one in (19)
we can once again exploit the presence of the vanishing root. In the large $N$
limit, the dominant contribution to $I_{3}$ is:

$I_{3}=\frac{2}{\pi}\frac{k_{1}a}{N}\int_{-\infty}^{\infty}\frac{dx}{x^{2}+a^{2}}\frac{x^{2}}{\sinh^{2}(x)},$
(80)

where we conveniently made the change of variable $x=y/k_{2}$ and $a$ is given
by:

$a=\frac{2k_{1}}{k_{2}N}=T_{N}/T,\quad T_{N}=\frac{\hbar
k_{1}\sqrt{2k/m}}{k_{B}N}.$ (81)

As in the classical case (19), we cannot directly take the limit
$N\rightarrow\infty$, because in this limit $a=0$ and (80) diverges.

Figure 11: The contour $\Gamma_{n}$ used to compute the integral $I_{3}$

To compute (80) we employ contour integration and Cauchy’s theorem. Let us
introduce the following function of complex variable $z$:

$f(z)=\frac{z^{2}}{z^{2}+a^{2}}\frac{1}{\sinh^{2}(z)},$ (82)

The poles of $f$ are all located on the imaginary axis, at the following
positions (with the corresponding residue):

$\displaystyle z_{n}=in\pi,\quad
n\in\mathbb{Z}\\{0\\},\,\quad\text{Res}_{n}=\frac{2i\pi
a^{2}n}{\left(a^{2}-\pi^{2}n^{2}\right)^{2}},$ (83) $\displaystyle z_{\pm
a}=\pm ia,\qquad\qquad\quad\quad\text{Res}_{\pm
a}=\mp\frac{ia}{2\sin^{2}(a)}.$ (84)

Consider now the contour $\Gamma_{n}$ plotted in figure 11: it is composed by
a segment $[-R,R]$ and a semicircle $\mathcal{C}_{n}$ of radius $R$, which is
such that $\Gamma_{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:

$I_{R}=\int_{-R}^{R}dxf(x)=-\int_{\mathcal{C}_{n}}dzf(z)+2\pi
i\left[\sum_{k}^{n}\text{Res}_{k}+\text{Res}_{+a}\right].$ (85)

$I_{3}$ can then be obtained by taking the limit $R\rightarrow\infty$ of
$I_{R}$ as follows:

$I_{3}=\frac{2}{\pi}\frac{k_{1}a}{N}\lim_{R\rightarrow\infty}I_{R}.$ (86)

We now have to compute the limit $R\rightarrow\infty$, which also entails the
limit $n\rightarrow\infty$ of the right-hand side of (85). For large $R$ the
integral of $f$ over $\mathcal{C}_{n}$ is:

$\int_{\mathcal{C}_{n}}dzf(z)\approx
iR\int_{0}^{\pi}\frac{e^{i\theta}d\theta}{\sinh^{2}(Re^{i\theta})}=iR(-2i\coth(R)/R)\rightarrow
2.$ (87)

The sum over the residues becomes a series that can be resummed. We can thus
finally express $I_{3}$ as:

$I_{3}=\frac{2k_{1}}{N}g(a),$ (88)

where the function $g(x)$ is given by (49).

## References

## References

  * [1] Rieder Z, Lebowitz J L and Lieb E 1967 Journal of Mathematical Physics 8 1073–1078 URL https://doi.org/10.1063/1.1705319
  * [2] Lepri S, Livi R and Politi A 2003 Physics Reports 377(1) 1–80
  * [3] Dhar A 2008 Adv. Phys. 57 457–537
  * [4] Dhar A and Saito K 2016 Heat transport in harmonic systems Thermal Transport in Low Dimensions (Springer) pp 39–105
  * [5] Zürcher U and Talkner P 1990 Physical Review A 42 3278
  * [6] Dhar A and Roy D 2006) Journal of Statistical Physics 125 801
  * [7] Asadian A, Manzano D, Tiersch M and Briegel H 2013 Physical Review E 87 012109
  * [8] Freitas N and Paz J P 2014 Physical Review E 90 042128
  * [9] Saito K, Takesue S and Miyashita S 2000 Physical Review E 61 2397
  * [10] Lepri S (ed) 2016 Thermal transport in low dimensions: from statistical physics to nanoscale heat transfer (Lect. Notes Phys vol 921) (Springer-Verlag, Berlin Heidelberg)
  * [11] Benenti G, Lepri S and Livi R 2020 Frontiers in Physics 8 292 ISSN 2296-424X URL https://www.frontiersin.org/article/10.3389/fphy.2020.00292
  * [12] Lepri S and Politi A 2011 Phys. Rev. E 83 030107
  * [13] Kundu A, Bernardin C, Saito K, Kundu A and Dhar A 2019 J. Stat. Mech: Theory Exp. 2019 013205 URL http://stacks.iop.org/1742-5468/2019/i=1/a=013205
  * [14] Roy D 2012 Phys. Rev. E 86(4) 041102 URL https://link.aps.org/doi/10.1103/PhysRevE.86.041102
  * [15] Gendelman O V and Savin A V 2000 Phys. Rev. Lett. 84(11) 2381–2384 URL https://link.aps.org/doi/10.1103/PhysRevLett.84.2381
  * [16] Giardinà C, Livi R, Politi A and Vassalli M 2000 Phys. Rev. Lett. 84(10) 2144–2147 URL https://link.aps.org/doi/10.1103/PhysRevLett.84.2144
  * [17] Toda M 2012 Theory of nonlinear lattices vol 20 (Springer Science & Business Media)
  * [18] Lepri S, Livi R and Politi A 2020 Phys. Rev. Lett. 125(4) 040604 URL https://link.aps.org/doi/10.1103/PhysRevLett.125.040604
  * [19] Politi A 2011 J. Stat. Mech.: Theory and Experiment P03028
  * [20] Prosen T and Žunkovič B 2013 Phys. Rev. Lett. 111(4) 040602
  * [21] Gupta S and Ruffo S 2017 International Journal of Modern Physics A 32 1741018 (Preprint https://doi.org/10.1142/S0217751X17410184) URL https://doi.org/10.1142/S0217751X17410184
  * [22] Defenu N, Donner T, Macrì T, Pagano G, Ruffo S and Trombettoni A 2021 (Preprint 2109.01063)
  * [23] Bouchet F, Gupta S and Mukamel D 2010 Physica A: Statistical Mechanics and its Applications 389 4389–4405
  * [24] Campa A, Dauxois T and Ruffo S 2009 Phys. Rep. 480 57–159
  * [25] Campa A, Dauxois T, Fanelli D and Ruffo S 2014 Physics of long-range interacting systems (OUP Oxford)
  * [26] de Buyl P, De Ninno G, Fanelli D, Nardini C, Patelli A, Piazza F and Yamaguchi Y Y 2013 Physical Review E 87 042110
  * [27] Torcini A and Lepri S 1997 Phys. Rev. E 55 R3805
  * [28] Métivier D, Bachelard R and Kastner M 2014 Phys. Rev. Lett. 112 210601
  * [29] Ávila R R, Pereira E and Teixeira D L 2015 Physica A: Statistical Mechanics and its Applications 423 51–60
  * [30] Olivares C and Anteneodo C 2016 Phys. Rev. E 94(4) 042117
  * [31] Bagchi D 2017 Phys. Rev. E 95 032102
  * [32] Bagchi D 2017 Phys. Rev. E 96 042121
  * [33] Iubini S, Di Cintio P, Lepri S, Livi R and Casetti L 2018 Phys. Rev. E 97(3) 032102 URL https://link.aps.org/doi/1
  * [34] Wang J, Dmitriev S V and Xiong D 2020 Physical Review Research 2 013179
  * [35] Di Cintio P, Iubini S, Lepri S and Livi R 2019 Journal of Physics A: Mathematical and Theoretical 52 274001
  * [36] Bagchi D 2021 Physical Review E 104 054108
  * [37] Jin T, Ferreira J S, Filippone M and Giamarchi T 2021 Exact description of quantum stochastic models as quantum resistors (Preprint %****␣draft2.bbl␣Line␣150␣****2106.14765)
  * [38] Tamaki S and Saito K 2020 Phys. Rev. E 101(4) 042118 URL https://link.aps.org/doi/10.1103/PhysRevE.101.042118
  * [39] Joshi M K, Kranzl F, Schuckert A, Lovas I, Maier C, Blatt R, Knap M and Roos C F 2021 Observing emergent hydrodynamics in a long-range quantum magnet (Preprint 2107.00033)
  * [40] Schuckert A, Lovas I and Knap M 2020 Phys. Rev. B 101(2) 020416 URL https://link.aps.org/doi/10.1103/PhysRevB.101.020416
  * [41] Spagnolli G, Semeghini G, Masi L, Ferioli G, Trenkwalder A, Coop S, Landini M, Pezzè L, Modugno G, Inguscio M, Smerzi A and Fattori M 2017 Phys. Rev. Lett. 118(23) 230403 URL https://link.aps.org/doi/10.1103/PhysRevLett.118.230403
  * [42] Keller T, Jäger S B and Morigi G 2017 Journal of Statistical Mechanics: Theory and Experiment 2017 064002 URL https://doi.org/10.1088/1742-5468/aa71d7
  * [43] Lepri S, Mejía-Monasterio C and Politi A 2010 Journal of Physics A: Mathematical and Theoretical 43 065002
  * [44] Sherman J and Morrison W J 1950 The Annals of Mathematical Statistics 21 124 – 127 URL https://doi.org/10.1214/aoms/1177729893
  * [45] Dauxois T, Lepri S and Ruffo S 2003 Communications in Nonlinear Science and Numerical Simulation 8 375–387
  * [46] Rego L G C and Kirczenow G 1998 Phys. Rev. Lett. 81(1) 232–235 URL https://link.aps.org/doi/10.1103/PhysRevLett.81.232
  * [47] Segal D, Nitzan A and Hänggi P 2003 The Journal of Chemical Physics 119 6840–6855 (Preprint https://doi.org/10.1063/1.1603211) URL https://doi.org/10.1063/1.1603211
  * [48] Defaveri L, Olivares C and Anteneodo C 2021 Heat conduction in chains of non-locally coupled harmonic oscillators: mean-field limit (Preprint 2111.01332)