\section*{\contentsname}%
  \@starttoc{toc}%
}
\catcode`@=1

\begin{document} 
\title[]{Classical and quantum harmonic mean-field models
coupled intensively and extensively with external baths}


\vspace{.5cm}

\author{Francesco Andreucci$^{1,2,*}$, Stefano Lepri$^{3,4,\dagger}$, Stefano Ruffo$^{1,2,3,\ddagger}$, Andrea Trombettoni$^{5,1,2,6,\#}$}
\address{$^{1}$SISSA - Via Bonomea 265, 34136 Trieste, Italy}
\address{$^{2}$INFN sezione di Trieste, 34100 Trieste, Italy}
\address{$^{3}$Consiglio Nazionale delle Ricerche, Istituto dei Sistemi Complessi, Via Madonna del Piano 10, I-50019 Sesto Fiorentino, Italy}
\address{$^{4}$Istituto Nazionale di Fisica Nucleare, Sezione di Firenze, Via G. Sansone 1 I-50019, Sesto Fiorentino, Italy}
\address{$^{5}$Department of Physics, University of Trieste, Strada Costiera 11, I-34151 Trieste, Italy}
\address{$^{6}$CNR-IOM DEMOCRITOS Simulation Center, Via Bonomea 265, I-34136 Trieste, Italy}
\begin{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.
\end{abstract}
Email:\\
*\href{mailto:fandreuc@sissa.it}{fandreuc@sissa.it}\\
$^{\dagger}$\href{mailto:stefano.lepri@isc.cnr.it}{stefano.lepri@isc.cnr.it}\\
$^{\ddagger}$\href{mailto:ruffo@sissa.it}{ruffo@sissa.it}\\
$^{\#}$\href{mailto:andreatr@sissa.it}{andreatr@sissa.it}\\
\vspace{.5cm}


\maketitle



\section{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:
\begin{equation} \label{kappa_def}
    \mathcal{J} = -\kappa \nabla T,
\end{equation}
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.~(\ref{kappa_def}): transport is ballistic since
quasi-particles travel undisturbed throughout the system bulk. 
This was demonstrated by Rieder, Lebowitz and Lieb \cite{RLL} 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 \cite{Lepri_review,DHARREV} and also
\cite{dhar2016heat} 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 \cite{Zuercher1990}, see 
e.g. \cite{dhar2006} and \cite{asadian2013heat, freitas2014analytic, saito2000energy}.


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 \cite{Lepri_review,DHARREV,Lepri2016,benenti2020}.
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 \cite{benenti2020}. 
In those cases, Eq.~(\ref{kappa_def}) should be replaced 
by its fractional version, yielding nonlinear temperature profiles across
the system \cite{Lepri2011b,Kundu2019}.  
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 \cite{roy2012}. 
If one conservation law is broken, as in the case of coupled rotors \cite{gendelman2000,giardina2000}, standard
diffusion is restored.
Finally, nonlinear integrable
models like for instance the Toda chain \cite{toda2012} should generically display ballistic transport 
mediated by solitons (see \cite{Kundu2019,lepri2020tooclose} and references
therein for some recent results). However, there may be exceptions to this 
intuition, as in the cases of hard-point particles \cite{Politi2011} and   
the classically integrable Landau-Lifshitz spin chain \cite{Prosen2013}.


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 \cite{ruffo2017} and \cite{defenu2021}, 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$ \cite{Bouchet2010,RuffoRev,Campa2014} and even lack of thermalization upon interaction with a single external bath \cite{deBuyl2013}.

Moreover, since perturbations may propagate with infinite velocities\cite{Torcini1997,Metivier2014} anomalous transport properties are expected. Indeed, numerical studies of energy transport in classical long-range interacting systems support this expectation \cite{avila2015length,Olivares2016,Bagchi2017,bagchi2017energy,iubini2018,wang2020thermal} 
The long-range version of the XY model \cite{Olivares2016,iubini2018} and the 
Fermi-Pasta-Ulam-Tsingou chain were analyzed in \cite{Bagchi2017,iubini2018,dicintio2019equilibrium,bagchi2021heat}, 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 \cite{jin2021};
we also refer once again to \cite{defenu2021} and references therein for further results. From an analytical perspective, a classical model with stochastic momentum exchange was studied in \cite{tamaki2020}, while a hydrodynamic approach for transport in long-range quantum magnets was developed in \cite{joshi2021} and \cite{schuckert2020}.

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 \cite{defenu2021} and for
  cold atoms in double wells described by the Lipkin-Meshkov-Glick model
  \cite{spagnolli}, with the stationary state of the Hamiltonian Mean Field model
  simulated with cold atoms in a cavity \cite{Keller2017}. Recent progress in such systems
  allows for the study
  of dynamical and transport properties, despite the effective inherent absence of spatial distance \cite{defenu2021}.

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:
        \begin{itemize}
        \item[$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.
        \item[$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.
        \end{itemize}
        In Sec. \ref{section2} we introduce the model, and the tecniques we are going to use to compute the flux and the temperature profile. In Sec.\ref{section3}, we study the classical case with intensive couplings, and we discuss the peculiarities of the temperature distribution, while in Sec. \ref{section4} we tackle the quantum case. In Sec. \ref{section5} and \ref{section6} 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.

        

\section{The model and the coupling to external baths} \label{section2}

We consider a Hamiltonian system describing a network of $N$ fully-connected harmonic oscillators with displacements $x_i$. The Hamiltonian of the model is:
        \begin{equation} 
                H=\frac{1}{2m}\sum_{i=1}^{N} p_{i}^{2}+\frac{k}{2N}\sum_{ij}(x_{i}-x_{j})^{2},
        \end{equation}
        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:
        \begin{equation} \label{ham_all2all}
                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},
        \end{equation}
where the matrix $\Phi$ is given by:
        \begin{equation} \label{phi_matrix}
                \Phi_{ij}=2k\left(\delta_{ij} -\frac{1}{N}\right).
        \end{equation}


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 \cite{dhar2006}. 
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 \cite{dhar2006}).
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:
\begin{align} \label{eom_fourier}
    -m\omega^{2} x_{i}(\omega) &= -\Phi_{ij}x_{j}(\omega)+ (\Gamma_{L,ij}(\omega)+\Gamma_{R,ij}(\omega))x_{j}(\omega)+\xi_{L,i}(\omega)+\xi_{R,i}(\omega),
\end{align}
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:
    \begin{equation} \label{fluctuation_dissipation}
        \Braket{\xi_{L,i}(\omega)\xi_{L, j}(\omega')} = \text{Im}[\Gamma_{L,ij}(\omega)] (1+f(\omega, T_{L,}))\frac{\hbar}{\pi}\delta(\omega+\omega'), 
    \end{equation}
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$,  
    \begin{align}
        &\Gamma_{L,ij}(\omega)=i\gamma\omega\delta_{ij}\delta_{i1},\ \quad &&\Gamma_{R,ij}=i\gamma\omega\delta_{ij}\delta_{1N},\label{gamma_intensive}\\
        &\Gamma_{L,ij}(\omega)=i\gamma\omega\delta_{ij}\sum_{k=1}^{N_{L}}\delta_{jk},\ \quad &&\Gamma_{R,ij}=i\gamma\omega\delta_{ij}\sum_{k=N-N_{R}+1}^{N}\delta_{jk}, \label{gamma_extensive}
    \end{align}
Eq. \eqref{gamma_intensive} defines the intensive coupling case, whereas \eqref{gamma_extensive} refers to the extensive coupling case
when both  $N_L$ and $N_R$ are taken to be proportional to $N$. 
By plugging \eqref{gamma_intensive} in \eqref{fluctuation_dissipation} we recover the  familiar fluctuation-dissipation relation, for example the classical limit of \eqref{fluctuation_dissipation} becomes:
    \begin{equation} \label{fluctuation_dissipation_standard}
        \Braket{\xi_{L/R,i}(\omega)\xi_{L/R,j}(\omega')} = \frac{\gamma k_{B}T_{L/R}}{\pi}\delta(\omega+\omega')\delta_{ij}.  
    \end{equation}
The solution of the equations of motion \eqref{eom_fourier} can be written as:
    \begin{equation} \label{eom_sol_fourier}
        x_{i}(\omega) = \sum_{j}G_{ij}(\omega) (\xi_{L,j}+\xi_{R,j}),
    \end{equation}
where the matrix $G_{ij}(\omega)$ is the Green's function, given by:
    \begin{equation} \label{green_def}
        G_{ij}(\omega) = \left[ -m\omega^{2}\mathbb{I}+\Phi+\Gamma_{L}(\omega)+\Gamma_{R}(\omega) \right]_{ij}^{-1}\equiv Z^{-1}(\omega).
    \end{equation}
Once $G_{ij}$ is computed, we can use the noise correlation \eqref{fluctuation_dissipation} 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 \eqref{green_def}.

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 \cite{RLL}, the covariances 
satisfy a Lyapunov equation, which can be solved numerically \cite{lepri2010nonequilibrium}. 
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. 


\section{Intensive coupling, classical case} \label{section3}
\subsection{Heat flux}
In the classical case the flux is given by:
	\begin{equation} \label{heat_flux_intensive_classic_def} 
		\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], 
        \end{equation}
where $\Delta T=T_{L}-T_{R}$. Plugging \eqref{gamma_intensive} in \eqref{heat_flux_intensive_classic_def} we get:
        \begin{equation} \label{heat_flux_intensive_classic}
		\mathcal{J}_{cl}^{int}  = \frac{k_{B}\Delta T}{\pi}\gamma^{2}\int_{-\infty}^{\infty}d\omega\, \omega^{2} |G_{1N}(\omega)|^{2}.
        \end{equation}
Now we need to compute the Green's function by inverting the following matrix:
        \begin{align} \label{z_matrix}
                 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}.
        \end{align}
In order to do it we can employ the Sherman-Morrison formula \cite{sherman50}, 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:
        \begin{equation} \label{sherman}
                 A^{-1} = M^{-1} +\frac{A^{-1}\bm{u}\bm{v}^{\,T}A^{-1}}{1+\bm{v}^{\,T}A^{-1}\bm{u}}.
        \end{equation}
It is easy to show that the $Z^{int}$ matrix \eqref{z_matrix} can be written as $-Z^{int} = D^{int} +\bm{u}\bm{u}^{T}$, with:
	\begin{eqnarray}
		&D^{int}_{ij} = \begin{cases} m\omega^{2}+i\gamma \omega -2k, \quad &i=j=1,N, \\ 
                                                             m\omega^{2}-2k, \quad &i=j\neq1,N, \\
                                                             0, \quad &i\neq j,
                           \end{cases}, \\ \nonumber \\
				\quad &u_{i} = \sqrt{2k/N} \,\, \forall i. \label{u_vector}
        \end{eqnarray}
        Using \eqref{sherman} we can then compute the Green's function \eqref{green_def} exactly. In particular we can compute $G_{1N}$ and plug the result in \eqref{heat_flux_intensive_classic_def} to obtain the heat flux:
        \begin{equation}
		\mathcal{J}_{cl}^{int}=\frac{k_{B}\Delta T\sqrt{2k/m}}{2\pi} I_{1}(k_{1},N),
        \end{equation}
where the function $I_{1}$ is given by the following integral $(y=\omega\sqrt{2k/m})$:
 	\begin{equation} \label{i1_def}
                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}},
        \end{equation}
where we introduced the dimensionless coupling constant:
        \begin{equation}
		k_{1}^{2}=\frac{\gamma^{2}}{2mk}.
        \end{equation}
Consider now the second factor of the denominator of the integrand in \eqref{i1_def}, as a polynomial in $s=y^2$:
    \begin{equation} \label{s_denomin}
    s( s-1)^{2}+k_{1}^{2}(s-2/N)^{2}.
    \end{equation}
As $N\rightarrow\infty$ \eqref{s_denomin} has a vanishing root:
    \begin{equation} \label{vanishing_root}
        s_{0}(k_{1},N)=-4k_{1}^{2}/N^{2}+o(N^{-2}),
    \end{equation}
where $o(x)$ indicates a quantity that goes to zero faster than $x$.
By decomposing the integrand in \eqref{i1_def} in partial fractions, it is easy to see that the dominant contribution for large $N$ is the one coming from \eqref{vanishing_root}:
        \begin{equation} \label{i1_large_n}
		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},
        \end{equation}
and so the heat flux at leading order in $N^{-1}$ reads as:
	\begin{equation} \label{heat_flux_intensive_classical_final}
		\mathcal{J}_{cl}^{int} = k_{B}\Delta T\sqrt{2k/m} \frac{k_{1}}{2N}.
	\end{equation}
 
 
        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. \ref{i1_num_plot}$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}$\footnote{It is worth noting that this dependence on the coupling constant is qualitatively the same as the one in short-range case \cite{RLL}.}. To better understand the problem,
          we also plot $NI_{1}(k_{1},N)$ as a function of $k_{1}$ in figures \ref{i1_num_plot}$a$ and \ref{i1_num_plot}$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 \ref{i1_large_n}) 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 \ref{i1_num_plot} 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. 

\begin{figure}
\centering
\subfigure
{\includegraphics[scale=0.30]{fig1a.pdf}}
\subfigure
{\includegraphics[scale=0.30]{fig1b.pdf}}        
\caption{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.}
\label{i1_num_plot}
\end{figure}


\subsection{Temperature profile}
We consider the kinetic definition of the temperature of the $i^{th}$ site:
	  \begin{equation} \label{temp_kin_def}
                  T_{i}  = \frac{m}{k_{B}} \Braket{\dot{x}_{i}^{2}}.
          \end{equation}
The classical velocity-velocity correlator can be expressed in terms of the Green's function as \cite{dhar2006}: 
\begin{equation} \label{v2_classical_def}
		\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),
	\end{equation}
and by substituting in \eqref{gamma_intensive} we get:
          \begin{align}
		  \Braket{\dot{x}_{i}^{2}} &= \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.
          \end{align}
For the first site, $i=1$, in the large $N$ limit  we get:
          \begin{align}
                  \Braket{\dot{x}_{1}^{2}} &= \frac{\gamma k_{B} T_{L}}{\pi} \int_{-\infty}^{\infty} \frac{\omega^2}{(m\omega^2-2k)^{2}+\gamma^2\omega^2}\nonumber +O(N^{-1}) =\frac{k_{B}}{m}T_{L} +O(N^{-1}),
          \end{align}
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:
          \begin{align} \label{v2_int_classic}
                  \Braket{\dot{x}_{i}^{2}} = \frac{k_{B}}{m}\frac{T_{L}+T_{R}}{2\pi} I_{2}(k_{1},k_{2},N),
          \end{align}
with:
		  \begin{align} \label{i2_def}
                          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}}.
                  \end{align}
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:
			  \begin{equation} \label{temp_profile_int_classic}
				  T_{cl,i}^{int} = \frac{T_{L}+T_{R}}{2\pi} \frac{I_{1}}{k_{1}} =\frac{T_{L}+T_{R}}{2  N}.
                          \end{equation}

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.

\begin{subsection}{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) 
    \begin{equation} \label{eom_argument}
        \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}),
    \end{equation}
where $\Phi$ is defined in \eqref{phi_matrix}, and $\xi_{L,R}$ are Gaussian noises with correlation given by \eqref{fluctuation_dissipation_standard}.  
We now introduce the ``total magnetization" $M(t) = \sum_{i} x_{i}/N$, and $S=x_{1}+x_{N}$. The equations of motion \eqref{eom_argument} can then be cast in the following form:
    \begin{align}
            &\ddot{S}=-\lambda S-2kS+4kM+\xi, \\
            &\ddot{M}=(\xi-\lambda \dot{S})/N, \\
            &\ddot{x}_{i} =-2k x_{i}+2kM, \quad i=2,...,N-1. \label{eom_xi}
    \end{align}
Switching to Fourier space, we find the following solution for $M(\omega)$ and the position $x_{i}(\omega)$ of the uncoupled sites:
    \begin{align}
        &M(\omega)=\frac{2k-\omega^{2}}{N\omega^{2}(\omega^{2}+i\lambda\omega-2k)-4ik\lambda\omega} \xi(\omega),\\
        &x_{i}(\omega)= \frac{2kM(\omega)}{2k-\omega^2}. \label{x_sol_eom}
    \end{align}
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. \eqref{x_sol_eom} as:
    \begin{equation} \label{x_sol_all2all}
            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)},
    \end{equation}
that has no pole on the dispersion law $\omega^{2}=2k$. The mean square velocity of the $i^{th}$ site then reads as:
    \begin{equation}
        \Braket{x^{2}_{i}} = \gamma  \frac{T_{L}+T_{R}}{\pi} \int d\omega |Q(\omega)|^{2},
    \end{equation}
which reproduces exactly formula \eqref{v2_int_classic}.\\
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 \eqref{eom_xi} entails that:
    \begin{equation}
        \ddot{z}_{i} = -2k z_{i},\quad i=2,...,N-2,
    \end{equation}
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. \ref{langprof} 
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 \eqref{temp_profile_int_classic}. 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.\\ 
\begin{figure} \centering
\includegraphics[scale=0.5]{fig2.pdf}
\caption{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. } 
\label{langprof}
\end{figure}
    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 (\ref{phi_matrix}). The Hamiltonians corresponding to these two choices are, respectively:
    \begin{align}
        &H_{1} = \sum_{i} \frac{p_{i}^{2}}{2} + \frac{k}{2 N_{\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} \label{H1} \\
        &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}. \label{H2}
    \end{align}
Note that in the case of Hamiltonian \eqref{H1} 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 \eqref{H2}, 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 \cite{RLL}. The results for the temperature profile are plotted in fig. \ref{temp_nondeg}. 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.\ref{langprof}, 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.\ref{eom_argument} \cite{dauxois2003clustering}. 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 \cite{iubini2018}: 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.

\begin{figure}
\centering
\subfigure
{\includegraphics[scale=0.8]{fig3a.pdf}}
\subfigure
{\includegraphics[scale=0.8]{fig3b.pdf}}
\caption{Numerical result for the temperature profile $T_{cl,i}^{int}$ for the Hamiltonians \eqref{H1} and \eqref{H2} 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$.}
\label{temp_nondeg}
\end{figure}
\end{subsection}    




\section{Intensive coupling, quantum case} \label{section4}
\subsection{Heat flux}
The heat flux in the quantum case is given by \cite{dhar2006}:
\begin{equation} \label{heat_flux_quantum_def}
		\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].
	\end{equation}
Furthermore, we will work in the linear response regime:
		\begin{equation}
			T_{L}-T_{R}=\Delta T \ll (T_{R}+T_{L})/2=T.
		\end{equation}
To get the heat flux for intensive couplings we plug \eqref{gamma_intensive} into \eqref{heat_flux_quantum_def} and we expand to first order in $\Delta T$:
		\begin{equation} 
			\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}.
		\end{equation}
The Green's function element $G_{1N}$ is the same as in the classical case, and therefore the heat flux can be written as:
		\begin{equation} \label{heat_flux_intensive_quantum}
			\mathcal{J}^{int}_{q} =  \frac{k_{B}\Delta T \sqrt{2k/m}}{4}  I_{3}(k_{1},\theta;N),
		\end{equation}
where we introduced a dimensionless temperature $\theta$ as:
    \begin{equation}
\theta=\frac{2k_{B}T}{\hbar}\sqrt{m/2k}, \label{k2_def}
    \end{equation}
and the function $I_{3}(k_{1},\theta,N)$ is given by:
		\begin{equation}  
			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}}       \label{i3_def}
                \end{equation}
		\begin{figure}  
            \centering
            \subfigure
            {\includegraphics[scale=0.3]{fig4a.pdf}}
            \subfigure
            {\includegraphics[scale=0.3]{fig4b.pdf}} 
			\caption{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 \eqref{g}.} \label{i3_num_plot}
                \end{figure}
The integral $I_{3}$ \eqref{i3_def} is plotted as a function of $\theta$ in fig. \ref{i3_num_plot}: 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:
	\begin{equation} \label{heat_flux_intensive_quantum_final}
		\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),
	\end{equation}
where the function $g(x)$ is given by:
\begin{align} 
 & 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, \label{g}\\
 &  \psi^{(1)}(z) = \frac{d^{2}}{dz^{2}} \Gamma(z), \nonumber
\end{align}
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:
    \begin{equation}
        T_{N}(k_{1}) = \frac{k_{1}}{N} \frac{\hbar}{k_{B}} \sqrt{\frac{2k}{m}}.
    \end{equation}
To get a better picture of the crossover from the quantum to 
the classical regime, we consider the ratio between the quantum \eqref{heat_flux_intensive_quantum_final}  and the classical heat flux \eqref{heat_flux_intensive_classical_final}, $\mathcal{J}_{q}^{int}/\mathcal{J}^{int}_{cl}$. For large $N$, this ratio is given by the function $g(x)$ defined in \eqref{g}. Its low and high temperature behaviors can be worked out explicitly are given by:
	\begin{equation} \label{g_asymptotics}
	g\biggr( \frac{T_{N}}{T} \biggr) =
		 \begin{cases} 
				1,  \quad &T\gg T_{N}, \\ \nonumber \\ \dfrac{\pi}{3}\dfrac{T}{T_{N}} \sim T N , \quad &T\ll T_{N},
			\end{cases} 
	\end{equation}
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. \ref{r.pdf} 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 \eqref{heat_flux_intensive_quantum_final} 
for low temperature is:
    \begin{equation}
        \mathcal{J}_{q}^{int} = \left( \frac{ \pi^{2} k_{B}^{2} T}{3 h} \right) \Delta T.
    \end{equation} \label{heat_flux_quantum_low_int}


Remarkably, the quantity among parentheses is recognized to be the quantum of thermal conductance, introduced in \cite{rego98} for heat transport in 
ballistic quantum wires. It is a universal quantity, independent of all the 
system parameters (the coupling constants in our case).
    
	\begin{figure} \centering
        \subfigure
            {\includegraphics[scale=0.3]{fig5.pdf}}
        \caption{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.} \label{r.pdf}
	\end{figure}
	
	
\subsection{Temperature profile}
In the quantum case, the velocity-velocity correlator is given by \cite{dhar2006}:
	\begin{equation} \label{v2_quantum_def}
		\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],
	\end{equation}
plugging \eqref{gamma_intensive} in \eqref{v2_quantum_def} we get: 
	\begin{equation}
         \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],
	\end{equation}
For site number $1$, that is, the one coupled to the left bath, the leading order term is:
	\begin{equation} \label{correlator_velocity_quantum_1n}
		\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).
	\end{equation}
 We get the same expression for site $N$, with $L\rightarrow R$.	
The integral in Eq. (\ref{correlator_velocity_quantum_1n}) 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 \cite{dvira2003}). 
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\neq1,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}$:
	\begin{equation}
		\Braket{\dot{x}_{i}^{2}}=\frac{\hbar \gamma}{2m^{2}}I_{4},
	\end{equation}
where $I_{4}$ is given by:
	\begin{equation} \label{i4_def}
		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}}.
	\end{equation}
	\begin{figure} \centering
		\includegraphics[scale=0.3]{fig6.pdf} 
		\caption{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.} \label{i4.pdf}
	\end{figure} 
While we cannot solve the integral $I_{4}$ exactly, we can compute its low-temperature behaviour by approximating the hyperbolic cotangent in \eqref{i4_def} with the following series, valid for large $x$:
    \begin{equation} \label{coth_series}
        \coth(x) = 1+2\sum_{k=0}^{\infty} e^{-2(k+1)x}.
    \end{equation}
The integral of the first term can be computed exactly (note that the denominator in \eqref{i4_def} is the same as in \eqref{i1_def}):
    \begin{equation}
        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],
    \end{equation}
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 \eqref{coth_series} 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:
    \begin{equation} \label{t_profile_quantum_int}
        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].
    \end{equation}
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.

\section{Extensive coupling, classical case} \label{section5}
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 \eqref{gamma_extensive}. The matrix that we need to invert is now:
		\begin{align} \label{z_ext_matrix}                                                          
                        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}  
		\end{align}
As in the intensive case, we can decompose this matrix as $-Z^{ext}=D^{est}+uu^{T}$, where $\bm{u}$ is defined in \eqref{u_vector}, and $D^{ext}$ is given by:
                        \begin{align} 
                                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}                                                                                       
			\end{align}
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 \textit{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. 

\subsection{Heat flux}
To compute the heat flux, we first substitute \eqref{gamma_extensive} in \eqref{heat_flux_intensive_classic_def} and then we plug in the relevant matrix elements of the Green's function:
	\begin{equation} \label{heat_flux_extensive_classic_preliminary}
		\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}),
	\end{equation}
where we defined $I_{5}$ as:
	\begin{equation} \label{i5_def}
		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] }.
	\end{equation}
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 \eqref{i1_def} of the intensive case by putting $\alpha_{R}=\alpha_{L}=1/N$ in \eqref{i5_def}. \\
In fig. \ref{i5_k1_alpha.pdf} 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 \eqref{heat_flux_intensive_classical_final}. 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 \ref{section3} 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}$.
	\begin{figure} \centering
		\includegraphics[scale=0.3]{fig7.pdf}
		\caption{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.} \label{i5_k1_alpha.pdf}
	\end{figure}

\subsection{Temperature profile}
The velocity-velocity correlator computed via \eqref{v2_classical_def} by substituting \eqref{gamma_extensive}:
	\begin{equation} \label{v2_correlator_classic_extensive}
		\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].
	\end{equation}
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:
	\begin{equation}
		T_{cl,i}^{ext}=\frac{\alpha_{L}T_{L}+\alpha_{R}T_{R}}{\pi N} I_{6}(k_{1},\alpha_{L},\alpha_{R}),
	\end{equation}
where we introduced the integral $I_{6}$:
	\begin{equation}
		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]},
	\end{equation}
if we set $\alpha_{L}=\alpha_{R}=1/N$ we recover the intensive case \eqref{i2_def}, apart from a factor of $2N$ due to the different definitions of these integrals.
	\begin{figure}  \centering
		\includegraphics[scale=0.3]{fig8.pdf}
		\caption{The plot of $I_{6}$ as a function of $\alpha_{L}+\alpha_{R}$ for $k_{1}=3$.} \label{i6.pdf}
	\end{figure}
In fig. \ref{i6.pdf} 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 \ref{section3} 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.
\section{Extensive coupling, quantum case} \label{section6}
\subsection{Heat flux}
The heat flux is obtained by substituting \eqref{gamma_extensive} in \eqref{heat_flux_quantum_def}. In the linear response regime we get, as in the intensive case, a factor related to the derivative of the Bose function:
	\begin{equation}
		\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}),
	\end{equation}
where we introduced the function $I_{7}$:
	\begin{equation} \label{i7_def}
		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}}.  
	\end{equation}
	\begin{figure} 
        \centering
        \subfigure
        {\includegraphics[scale=0.30]{fig9a.pdf}}
        \subfigure
        {\includegraphics[scale=0.30]{fig9b.pdf}}
        \caption{ In fig. $(a)$ we plot $I_{7}$ \eqref{i7_def} 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 \eqref{i7_low_temp}.}\label{i7_k1_alpha.pdf}
	\end{figure}
In fig. \ref{i7_k1_alpha.pdf}$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}$ \eqref{i3_def}, but we can obtain an estimate for the low-temperature behaviour using the following result:
	\begin{equation}
		\lim_{\theta \rightarrow 0} \frac{3}{\pi^{2} \theta^{3}} \frac{y^{2}}{\sinh^{2}(y/\theta)} = \delta(y), \label{nascent_delta}
	\end{equation}
so that $I_{7}$ for small $\theta$ is given by:
	\begin{equation} \label{i7_low_temp}
	 I_{7} = \frac{\pi^{2}}{3} \frac{\alpha_{R}\alpha_{L}}{(\alpha_{R}+\alpha_{L})^{2}}\theta,
	\end{equation}
and the heat flux vanishes linearly with the temperature, as in the intensive case \eqref{heat_flux_intensive_quantum_final}. In fig. \ref{i7_k1_alpha.pdf}$b$, we report the numerical exact plot of $I_{7}$ and the low-temperature approximation \eqref{i7_low_temp} 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 \eqref{i7_def} exactly. As in the intensive case \eqref{heat_flux_quantum_low_int}, at low temperature the flux can be expressed in terms of the quantum of thermal conductance:
    \begin{equation}
        \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,
    \end{equation}
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 \eqref{heat_flux_quantum_low_int}.

\subsection{Temperature profile}
In the quantum case we have to plug \eqref{gamma_extensive} in \eqref{v2_quantum_def}:
	\begin{align}
		\Braket{\dot{x}_{i}^{2}} = \gamma \biggr[ &\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) \nonumber \\
		 &+ \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].
	\end{align}
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:
	\begin{equation} \label{temp_quantum_extensive}
		T_{q,i}^{ext}=\frac{\hbar\gamma}{2\pi N m k_{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],
	\end{equation}
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:
	\begin{equation} \label{i8_def}
		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} }.
	\end{equation}
In the linear response regime we can Taylor expand \eqref{temp_quantum_extensive} around $\Delta T=0$ for $T_{L,R} = T \pm \Delta T/2$:
	\begin{equation} \label{temp_quantum_ext}
		 T_{q,i}^{ext} =\frac{\hbar\gamma}{2\pi N m k_{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],
	\end{equation}
where:
	\begin{eqnarray}
		&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]
                                }, \\
                                &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]
                                }.
	\end{eqnarray}
In fig. \ref{i8}$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 \eqref{nascent_delta}:
	\begin{equation} \label{i8_lin}
		I_{8}^{(1)} = \frac{\pi^2}{3 k_{1}^{2} (\alpha_{L}+\alpha_{R})^{2}}\theta.
	\end{equation}
In fig. \ref{i8}$b$ we report the plot of $I_{8}^{(1)}$ and \eqref{i8_lin} and we see that there is good agreement. Note that \eqref{temp_quantum_ext}  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 \eqref{t_profile_quantum_int}, the intensive counterpart of \eqref{temp_quantum_ext}, and can be obtained from \eqref{temp_quantum_ext} 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:
\begin{equation}
    T_{q,i}^{ext} = \frac{\hbar\gamma}{2\pi N m k_{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].
\end{equation}
\begin{figure}
    \centering
        \subfigure
        {\includegraphics[scale = 0.30]{fig10a.pdf}}
        \subfigure
        {\includegraphics[scale = 0.30]{fig10b.pdf}}
        \caption{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 \eqref{i8_lin}} \label{i8}
\end{figure}

\section{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   \ref{table}, 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). \\
\begin{table}
    \centering
    \begin{tabular}{|c|c|c|c|}
    \hline
 Coupling    & Quantity    & Classical    & Quantum \\
 \hline
 \multirow{4}*{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})$ \\ & & &\\

 \hline
 \multirow{4}*{Extensive}
                         &$\mathcal{J}^{ext}$ & $ \mathcal{J}_{cl}^{ext} \sim \Delta T N^{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$ \\  
                         & & &\\

\hline 
    \end{tabular} 
    \caption{Summary of main results about thermal transport in the mean-field harmonic model \eqref{ham_all2all}. We only report the low-temperature behaviour of the quantum results.}  \label{table}
    \end{table}
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 \ref{temp_nondeg}. 
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 \cite{rego98}, 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.
\\ \\
\textit{Note added:} During the completion of this manuscript, an interesting and related paper by L. Defaveri, C. Olivares and C. Anteneodo \cite{defaveri2021} 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.

\ack

SL acknowledges support from the program 
\textit{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 \textit{Coarse-grained description for non-equilibrium systems and transport phenomena (CO-NEST)} No. 201798CZL.


\vspace{40px} 

