Anomalous heat transport in classical many-body systems: overview and perspectives

In this review paper we aim at illustrating recent achievements in anomalous heat diffusion, while highlighting open problems and research perspectives. We briefly recall the main features of the phenomenon for low-dimensional classical anharmonic chains and outline some recent developments on perturbed integrable systems, and on the effect of long-range forces and magnetic fields. Some selected applications to heat transfer in material science at the nanoscale are described. In the second part, we discuss of the role of anomalous conduction on coupled transport and describe how systems with anomalous (thermal) diffusion allow a much better power-efficiency trade-off for the conversion of thermal to particle current.


I. INTRODUCTION
Anomalous diffusion is a well-established concept in statistical physics and has been invoked to describe many diverse kinetic phenomena. A very detailed insight has been achieved by generalizations of the motion of Brownian particles as done for the continuous-time random walk, Lèvy flights and walks. A formidable body of literature on the topic exists, and we refer for example to [1] as well as the present issue for an overview.
The above particle models are based on a single-particle description, whereby the single walker performs a nonstandard diffusive motion. How do these features emerge when dealing with a many-body problem? What are the conditions for a statistical system composed of many interacting particles to yield effectively anomalous diffusion of particles or quasi-particles? A further question regards how such anomalies in diffusion can be related to transport and whether they can be somehow exploited to achieve some design principle, like efficiency of energy conversion. In fact, although thermoelectric phenomena are known since centuries, it is only recently that a novel point of view on the problem has been undertaken [2]. Generally speaking, this renewed research activity is motivated also by the possibility of applications of the thermodynamics and statistical mechanics to nano and micro-sized systems with applications to molecular biology, micro-mechanics, nano-phononics etc. This requires treating systems far from thermodynamic limit, where fluctuations and interaction with the environment are essentially relevant and require to be treated in detail.
In the present contribution, we first review how anomalous energy diffusion arises in lattices of classical oscillator as a joint effect of nonlinear forces and reduced dimensionality (and in this respect we will mostly discuss onedimensional chains). In words, this amounts to say that the anomalous dynamics of energy carriers is an emergent feature stemming from correlations of the full many-body dynamics. As a consequence, Fourier's law breaks down: motion of energy carriers is so correlated that they are able to propagate faster than diffusively. In the second part of the paper, we discuss how this feature influences coupled transport and how it can be used to enhance the efficiency of thermodiffusive processes.
We conclude this review about the multifaceted problem of heat transport in classical systems with a short overview summarizing possible extensions to the quantum domain, with reference to related open problems, which certainly will deserve attention in the near future.

II. ANOMALOUS HEAT TRANSPORT IN CLASSICAL ANHARMONIC LATTICES
The presence of a diverging heat conductivity with the systems size in a chain of coupled nonlinear oscillators was first pointed out in [3,4]. This was the beginning of a research field that, over more than two decades, has been devoted to the understanding of the mechanisms yielding anomalous transport in low-dimensional systems. Far from being a purely academic problem, it has unveiled the possibility of observing such peculiar effects in nanomaterials, like nanotubes, nanowires or graphene [5,6]. Extended review articles on this problem exist since many years [7,8], while a collection of contributions about more recent achievements is contained in [9] and in a review article of the present issue [10]. Here it is useful to provide first a short summary about the state of the art in this field, while the main effort will be devoted to some recent achievements that point to promising and challenging directions for future investigations.
In any model where anomalous transport has been observed it emerges as a hydrodynamic effect due to the conspiracy of reduced space-dimensionality and conservation laws, yielding non standard relaxation properties even in a linear response regime. As a reference we may consider the basic class of models, represented by a Hamiltonian of the following form: Typical choices for the interaction is the famous Fermi-Pasta-Ulam-Tsingou (FPUT) potential where V (x) = V F P U T (x) ≡ 1 2 x 2 + α 3 x 3 + β 4 x 4 and rotor (or Hamiltonian XY) model For what concerns conservation laws, in d = 1 anomalous transport has been generically observed in Hamiltonian models of type (1), where energy and momentum and the "stretch" variable n (q n+1 − q n ) are conserved. It is worth recalling that any approach aiming at describing out-of-equilibrium conditions, like stationary transport processes, has to be based on the hydrodynamic equations associated with such locally conserved quantities. One-dimensional oscillator models with only one (e.g. the Frenkel-Kontorova or φ 4 models [11]) or two conserved quantities like the rotor model [12,13] or the discrete nonlinear Schrödinger lattice [14,15] show instead standard diffusive transport. Intuitively, this is due to the presence of scattering sources for acoustic waves propagating through the lattice imposed by the very presence of a local nonlinear potential, which breaks translation invariance (i.e. momentum conservation). This argument does not apply to the rotor model, where stretch only is not conserved, due to the angular nature of the q n variables: anyway standard diffusion is allowed because of the boundedness of the cosine potential. For what concerns dimensionality, in d = 3 normal diffusion regimes are expected to characterize heat transport in nonlinear lattices. Only in d = 2 one can find evidence of a diverging heat conductivity according to a logarithmic dependence with the system size L [16][17][18].
The main distinctive feature of anomalous heat transport in one-dimensional Hamiltonian models of anharmonic lattices is that the finite-size heat conductivity κ(L) diverges in the limit of a large system size L → ∞ [3] as with 0 < γ ≤ 1, (the case γ = 1 corresponding to integrable models, like the Toda lattice discussed in Sec. II.B). This implies that this transport coefficient is,in the thermodynamic limit, not well-defined. In the linear response regime, this is equivalent to find that the equilibrium correlator of the energy current J(t) displays for long times t a nonintegrable power-law decay, with 0 ≤ δ < 1. Accordingly, the Green-Kubo formula yields an infinite value of the heat conductivity and allows to establish the equivalence of the exponents, i.e. γ = δ, provided sound velocity is finite [4]. In Fig.1 we illustrate two typical simulations of the FPUT model demonstrating the results above. The most basic issue of anomalous feature is related to anomalous dynamical scaling of equilibrium correlation of the hydrodynamic modes. A simple way to state this is to say that fluctuations of the conserved quantities with small wavenumber k evolve on time scales of order τ (k) ∼ |k| −z . For standard diffusion one has z = 2. Within the nonlinear fluctuating hydrodynamics approach it has been shown [19,20] that models like (1) belong generically to the universality class of the famous Kardar-Parisi-Zhang (KPZ) equation, originally formulated in the context of growing interfaces. It is well known that this equation in d = 1 is characterized by the dynamical exponent z = 3/2. The origin of this nontrivial dynamical exponent can be traced back to the nonlinear interaction of long-wavelength  (2) corresponds to a divergence ω −δ at small frequencies.
modes. This leads to the prediction γ = (2 − z)/z = 1/3 (at least in the linear response regime), a value that should be largely universal, as confirmed by many numerical experiments. The above consideration applies generically for anharmonic chains with three conservation laws [20]. There is however the possibility of having a different universality class depending on the number of conserved quantities [21] and/or on the nonlinear coupling among the hydrodynamic modes [20]. For, instance model (1) with an even potential V (x) = V (−x) should belong to the a different universality class having a different exponent γ. Actually, the precise value γ is still somehow controversial: the theoretical prediction from mode-coupling approximation of the hydrodynamic theory yields γ = 1/2 [20,22,23] while kinetic theory yields γ = 2/5 [24], a value closer to the one measured numerically [25,26] (see also [27,28] for related results on exactly solvable models). The existence of the two classes can be demonstrated either by direct measurements of the exponents [25], but also by suitable changes of the thermodynamic parameters. For instance, a nonlinear chain with a symmetric potential subject to a suitable pressure acting at its boundaries may exhibit a change from exponent γ = 1/2 to γ = 1/3 [29]. This is a relevant observation for possible experimental verification of anomalous heat transport, because the pressure or torque applied to any one-dimensional material should be taken into account for a correct comparison with theoretical predictions.
A physically intuitive way to describe anomalous heat transport is to think in terms of a Lévy walk, namely an ensemble of random walkers performing free ballistic steps with finite velocity for times that are power-law distributed [30]. This simple description accounts very well for many features of anharmonic lattices and fluids in various nonequilibrium setups [31][32][33]. For instance energy perturbations propagate superdiffusively [31,34,35]: an initially localized perturbation of the energy broadens and its variance grows in time as σ 2 ∝ t η with η > 1. These empirical observation have a theoretical justification within the nonlinear fluctuating hydrodynamics. Indeed the theory predicts a hydrodynamic "heat mode" that has the characteristic shape given by a Lévy-stable distribution, see [20,36] for details. Further support comes from mathematical results: superdiffusive behavior has been proven for one dimensional infinite chain of harmonic oscillators undergoing stochastic collisions conserving energy and momentum [37,38]. In the same spirit, the more difficult case of nonlinear oscillators with conservative noise has been discussed [39]. For exponential interactions (the Kac-van Moerbecke model), superdiffusion of energy is again demonstrated, and a lower bound on the decay of current correlation function is given [40]. In Ref. [20] it is argued also such models should belong to the KPZ class.
One related distinctive feature of anomalous transport is that the temperature profiles in nonequilibrium steady states are nonlinear, even for vanishing applied temperature gradients [32,41]. There is indeed a close connection with fractional heat equation, that has been demonstrated and discussed in the recent literature [10,42].

A. The importance of being small
As mentioned in the previous section, theoretical inferences on the problem of heat transport in anharmonic chains stems from the basic assumption that one should compute any relevant quantity in the limits L → ∞ and t → ∞, performed just in this order. On the other hand, in any numerical simulation or in real low-dimensional heat conductors, such as nanowires, carbon nanotubes, polymers or even thin fibers, we have to deal with finite size and finite time corrections. These can be taken under control in a linear response regime if the mean-free path of propagating excitations, λ, and their mean interaction time, τ , are such that λ L and τ t. It is a matter of fact that, when dealing with models of anharmonic chains, such a control is often not granted, mainly because of nonlinear effects. This is may be a very relevant problem also for interpreting possible experimental verifications of anomalous transport in real systems and also for designing nanomaterials, exhibiting deviations from the standard diffusive conductivity.
As a matter of fact, very severe finite-size effects invariably arise when trying to check predictions numerically. Very often, the estimate of the relevant exponents γ or δ systematically deviates from the expected values, and sometimes are even claimed to depend on parameters [43][44][45][46]. If universality is (as we do believe) to hold, this effects should be due to subleading corrective terms to the asymptotics that are still relevant on the scales accessible to simulations. Besides those issues, other unexpected effects arise. For instance, for the FPUT [47], Toda [43] and the Kac-van Moerbecke [38] chains perturbed by conservative noise, the exponent γ increases with the noise strength. Besides the problem of evaluating the precise exponents, this is pretty surprising since it suggest that a larger stochasticity in the model makes the system more diffusive, at least for finite systems.
Another instance where finite size corrections are "amplified" by nonlinear effects is the case of anharmonic chains with asymmetric potential , i.e. with V (x) = V (−x) as the FPUT model with α = 0. As already shown in Fig.1, both equilibrium and out of equilibrium measurements of the heat conductivity in the presence of an applied thermal gradient are usually consistent with KPZ scaling. However, in other temperature regimes, Fourier's law appears to hold, i.e. thermal conductivity is constant over a large range of sizes [48]. This was traced back to the relatively long relaxation time of mass inhomogeneities induced by the asymmetry of the interaction potential and acting as scatterers of phonons [48]. Actually, it was later shown in [49][50][51] that this is a strong finite-size effect, since it persists for relatively large values of L and t. Yet, the expected theoretical prediction of a diverging heat conductivity can be recovered in simulations performed for sufficiently large values of L and t. It should be pointed out that all of these speculations rely on numerical results, while a theoretical approach capable to provide estimates of the combination of nonlinear with finite-size corrections to the hydrodynamics would be useful. Indeed, fluctuating hydrodynamics provide some prediction: subleading corrections to the leading asymptotic decay Eq. (2) can be very large and decay very slowly [20].
Finite size-effects can have other facets. A remarkable example is the Discrete Nonlinear Schroedinger equation, a well-known model for atomic condensates in periodic optical lattices. The model has two conserved quantities (energy and number density) and display normal diffusive transport [14]. However, at very low temperatures there appears a further almost conserved quantity (the phase difference among oscillators) and, for a finite chain and long times the dynamics is the the same of a generic anharmonic model, leading to KPZ scaling of correlations and anomalous transport [52,53]. Further unexpected features have been reported also in [54]. In this paper the authors study this problem for a the FPUT-β model (i.e. Eq. (1) with V = V F P U T and α = 0) with an additional local, also called "pinning", potential of the form This term breaks translational invariance making energy the sole conserved quantity. By varying the nonlinear coupling β one observes a crossover from a ballistic transport, typical of an integrable model, to an anomalous diffusive regime ruled by and exponent of the time correlation function, which corresponds to a value of γ ∼ 0.2. The crossover occurs in the parameter region 0.1 < β < 1. Numerical simulations performed for a chain of a few thousands of oscillators show that further increasing β seems to yield an increasing γ. The overall outcome challenges the basic theoretical argument, which predicts that an anharmonic chain equipped by a local potential should exhibit normal diffusion. For the sake of completeness it is worth mentioning that in this paper the model under scrutiny is compared with the so-called φ 4 -model, where the nonlinear term in Hamiltonian (1), i.e. β(q n+1 − q n ) 4 , is substituted with βq 4 n . Also for this model one observes, in the same parameter region, a crossover from a ballistic regime to an anomalous diffusive regime, but for β > 1 one eventually obtains numerical estimates yielding γ ∼ 0, i.e. the expected diffusive behavior is recovered.
All of these results have a logical interpretation only if we invoke, once again, the role of finite size corrections combined with nonlinearity. Actually, in the φ 4 -model there is no way to argue that a ballistic regime should be observed for any finite, even if small, value of β. The ballistic behavior observed in both models for β < 0.1 seems to suggest that for small nonlinearities one needs to explore definitely much larger chains and integrate the dynamics over much longer times, than those employed in [54], before in both chains phonon-like waves may experience the scattering effects due to the local potential. Moreover, the weaker quadratic pinning potential of the original model seems to be still affected by finite size corrections, even in the region β > 1. A problem that one should investigate systematically is the dependence on β of the chain length and of the integration time necessary to recover standard diffusive transport, at least in the crossover region 0.1 < β < 1, where one can expect to perform proper numerical analysis in an accessible computational time.

B. Chimeras of ballistic regimes
In the light of what discussed in the previous section one should not be surprised to encounter further nonlinear chain models, equipped with a pinning potential, which exhibit a regime of ballisitic transport of energy, compatible with a linearly divergent heat conductivity, κ(L) ∼ L. Again, one could conjecture that this is due to the puzzling combination of nonlinearity and finite size effects, although, as we are going to discuss, the emerging scenario is more intricate and interesting than the former statement could foretell.
As a preliminary remark we want to recall that ballistic transport is the typical situation of an integrable Hamiltonian chain, whose prototypical example is the harmonic lattice, with V (x) = 1 2 x 2 in Hamiltonian (1). It is worth pointing out that the addition of the harmonic pinning term (3) again keeps the harmonic chain integrable. In this perspective, it seems reasonable to find that for sufficiently small nonlinearities both the FPUT-β and the φ 4 chains, as discussed in the Section II.B, may exhibit a seemingly ballistic regime for β < 0.1 also in the presence of the pinning potential. Recovering the expected diffusive transport regime is a matter of simulating exceedingly large chains over exceedingly long times.
Anyway, the peculiar role played by the quadratic pinning potential (3) has emerged also in a recent study of heat transport in the Toda chain [55]. It is worth recalling that the unpinned Toda chain is an integrable Hamiltonian model of the form (1) with V (x) = e −x + x − 1: in this model heat transport is ballistic due to the finite-speed propagation of solitons (rather than phonons, as in the harmonic chain). Toda solitons are localized nonlinear excitations which are known to interact with each other by a non-dissipative diffusion mechanism: the soliton experiences a random sequence of spatial shifts while it moves through the lattice and interacts with other excitations without exchanging momentum [56]. In fact, the calculation of the transport coefficients by the Green-Kubo formula indicates the presence of a finite Onsager coefficient, which corresponds to a diffusive process on top of the dominant ballistic one [57,58].
When the pinning term (3) is at work the Toda chain becomes chaotic, as one can easily conclude by measuring the spectrum of Lyapunov exponents [55]. This notwithstanding, not only the energy, but also the "center of mass" is a conserved quantity. The peculiarity of the quadratic pinning potential (3) is testified by the observation that if one turns it into a quartic one, i.e. U (q i ) = 1 4 L n=1 q 4 n , the quantity h c is no more conserved. Nonequilibrium simulations of the Toda chain with the addition of (3), where heat reservoirs at different temperature, T 1 > T 2 act at its boundaries, yield a "flat" temperature profile at T = (T 1 + T 2 )/2 in the bulk of the chain. This scenario would be expected for the pure Toda chain, but it is inconsistent with the basic consideration that the presence of (3) breaks translation invariance and total momentum is no more conserved. This notwithstanding, in order to observe a temperature profile in the form of a linear interpolation between T 1 and T 2 (i.e., Fourier's law) one has to simulate the dynamics of very large chains over very long times: typically L ∼ O(10 4 ) and t ∼ O(10 6 ), when all the parameters of the model are set to unit.
Equilibrium measurements based on the Green-Kubo relation, i.e. on the behavior of the energy current correlator (2), provide further interesting facets of this scenario. By comparing the Toda chain with quadratic and quartic pinning potentials one observes in the latter case clear indications of a diffusive regime, i.e. a finite heat conductivity, and a practically negligible influence of finite size corrections, while in the former case the power spectrum (i.e. the Fourier transform of (2)) is found to exhibit a peculiar scaling regime (with a power −5/3), before eventually reaching a plateau, that indicates a standard diffusion. In the same region of the spectrum the FPUT model (where the parameters α and β have been chosen in such a way to correspond to a Taylor series expansion of the Toda chain) with the addition of (3) is found to converge to a plateau, in the absence of any precursor of a power-law scaling.
Further details about the unexpected transport regimes encountered in the Toda chain equipped with the quadratic pinning can be found in [59]. One should point out that many of them are still waiting for a convincing theoretical interpretation.
C. Anomalous transport in the presence of a magnetic field Recent contributions [60,61] have tackled the important problem of heat transport in chains of charged oscillators in the presence of a magnetic field B. The model at hand is a one-dimensional polymer, that allows for transverse motion of the oscillators interacting via a harmonic potential. For B = 0 the exponent of the energy current correlator is δ = 1 2 , thus indicating the presence of anomalous transport and a divergent heat conductivity κ(L) ∼ L 1 2 . By switching on B the first basic consequence is the breaking of translation invariance, so that the total momentum is no more conserved. This notwithstanding, the total pseudo-momentum is conserved, but the hydrodynamics of the model is definitely modified. Actually, numerical and analytic estimates indicate that the exponent δ may turn to a value different from 1 2 . In particular, in [60] two different cases were considered: the one where oscillators have the same charge and the one where oscillators have alternate charges of sign (−1) n , n being the integer index numbering oscillators along the chain. It can be easily shown that in the former case the sound velocity is null and the energy correlator exhibits a thermal peak centered at the origin and spreading in time. Conversely, in the latter case the sound velocity has a finite, B-dependent value and the thermal peak of the energy correlator is coupled to sound modes propagating through the chain. In the case of finite sound velocity (alternate charges) the exponent ruling the divergence of the heat conductivity with the system size is found to remain the same obtained for B = 0, i.e. γ = 1 2 . This is not surprising, since also for B = 0 the sound velocity in the model is finite. In the case of equally charged bids, on the other hand, it appears a new exponent γ = 3 8 , which corresponds to a universality class, different from all the others encountered in anomalous transport in nonlinear chains of oscillators. An important remark on this new exponent is that, in the absence of a finite sound velocity, the identification of the exponents δ and γ, introduced in Sec.II, is no more correct. In fact, in this case the value of the exponent δ is found to be very close to 3 4 . Rigorous estimates of all of these exponents and also for the d = 2 and d = 3 versions of the charged polymer model have been obtained through the asymptotics of the corresponding Green-Kubo integrals, where the deterministic dynamics has been substituted with a stochastic version that conserves the same quantities [61] . For what concerns the different one-dimensional cases, these rigorous estimates agree with the previous findings, while in d = 2 and in d = 3 the expected logarithmic divergence and a finite heat conductivity have been singled out, respectively.

D. The case of long-range interactions
Long-range forces, slowly decaying in the relative distance between particles are well-studied in statistical mechanics. They characterize a wide domain of physical situations, e.g. self-gravitating systems, plasmas, interacting vortices in fluids, capillary effects of colloids at an interface, chemo-attractant dynamics, cold atoms in optical lattices, colloidal active particles etc. Several unusual features are known: ensemble inequivalence, long-living metastable states and anomalous energy diffusion [62,63], inhomogeneous stationary states [64], lack of thermalization on interaction with a single external bath [65], etc. Moreover perturbations can spread with infinite velocities, leading to qualitative difference with respect to their short-ranged counterparts [66,67].
The study of heat transport in chains with long-range interactions has been tackled only recently [68][69][70][71][72]. The main question is to what extent are the anomalous properties changed when the spatial range of interaction among oscillators increases. In two recent papers [73,74] this problem has been investigated for Hamiltonian chains with a long range potential of the form where the generalized Kac factor For what concerns the rotor chain, nonequilibrium measurements, with thermal reservoirs at different temperature, T 1 > T 2 , acting at the chain ends, show that when α > 1, the resulting temperature profile interpolates linearly between T 1 and T 2 . Despite the long range nature of the interaction this is a strong indication that a standard diffusive process still governs energy transport through the rotor chain, as in the limit α → +∞. Conversely, when α < 1 the temperature profile progressively flattens, until reaching a constant bulk temperature T = (T 1 + T 2 )/2 in the "mean-field" limit α → 0 + . It is important to point out that, at variance with the situation of chimeras of integrable models discussed in Section II.B, such flat temperature profiles have nothing to do with integrability, but they are driven by the dominance of a "parallel" energy transport mechanism, that connects the heat baths at the chain boundaries with each other directly through each individual rotors in the bulk of the chain. Energy transport along the chain is practically immaterial and the overall process is mediated by rotors, that have to compromise between the two different temperatures imposed by the reservoirs. A sketch of this mechanism is presented in figure  2. For small α each lattice site in the bulk is directly coupled to both thermal baths and its temperature sets to the average (T 1 + T 2 )/2 independently of the system size. Moreover, the average heat current exchanged with any other site in negligibly small.
At least for α < 1, a similar scenario seems to characterize the FPUT model : flat temperature profiles are observed also in this case and one can verify that the same parallel transport mechanism described for rotors is at work. On the other hand, the scenario is definitely more complicated for α > 1. Careful numerical studies, exploring finite size effects, indicate an overall scenario where an anomalous diffusion mechanism sets in, characterized by an exponent γ, that is expected to increase up to the one of the quartic FPUT model in the limit α → +∞. But one is faced with a first surprise for α = 2, where a flat temperature profile is restored, although, as numerics clearly indicates, the mechanism of transport along the chain certainly dominates over the parallel transport process. In the light of what discussed in the Section II.B this appear as a possible manifestation of a chimera ballistic regime, although there is no simple argument allowing us to invoke a relation of this special case with an integrable approximation, if any. The fact that the case α = 2 is characterized by a somehow "weaker nonintegrability" is confirmed also for a related model [70,72]. This can be traced back to the fact that in this case the lattice supports a special type of free-tail localized excitations (travelling discrete breathers) that enhance energy transfer [72].
A complementary approach is the analysis of space-time scaling of equilibrium correlations, that in the short-range case yields useful information by the dynamical exponent z [20]. A numerical study of the structure factors of the FPUT model [74], shows that for α > 1 the dynamical exponent z depends on α in a way that is certainly different form the one that could be expected on the basis of the theory of Lévy processes. Moreover, upon adding the cubic term 1 3 (q i − q j ) 3 to the potential one recovers the same dependence of z on α, up to α = 5. This is again a surprise, because in the limit α → +∞ the cubic and quartic versions of the FPUT-model should converge to different values of z. It is a matter of fact that no theoretical explanation is, at present, available to cope with this challenging scenario: in particular no hydrodynamic description is, to the best of our knowledge, available.

E. Anomalous transport via the Multiparticle Collision Method
So far we have discussed the case of lattice models. To test the generality of the results and their universality, it is important to consider more general low-dimensional many-body systems like interacting fluids or even plasmas. Although molecular dynamics would be the natural choice, it is computationally convenient to consider effective stochastic processes capable to mimic particle interaction through random collisions. A prominent example is the Multi-Particle-Collision (MPC) simulation scheme [75]. which has been proposed to simulate the mesoscopic dynamics of polymers in solution, as well as colloidal and complex fluids. Another application is the modeling parallel heat transport in edge tokamak plasma [76]. Indeed, in the regimes of interest of magnetic fusion devices, large temperature gradients will build up along the field line joining the hot plasma region (hot source), and the colder one close to the wall (that act as a sink). Besides this motivation, we mention here some results that are pertinent to the problem of anomalous transport.
In brief, the MPC method consists in partitioning the set of N p point-particles into N c disjoint cells. Within each cell, the local center of mass coordinates and velocity are computed and a rotation of particle velocities around a random axis in the cell's center of mass frame is performed. The rotation angles are fixed by imposing that the conserved quantities (energy, momentum) are locally preserved. Then, all particles are then propagated freely. Physical details of the interaction can be easily included phenomenologically, for instance introducing energy-dependent collision rates [18]. Interaction with external reservoirs can be also included by imposing Maxwellian distributions of velocity and chemical potentials on the thermostatted cells [77].
For the case of a one-dimensional MPC fluid, since the conservation laws are the same as say, the FPUT model, we expect it to belong to the same KPZ universality class of anomalous transport. Indeed, numerical measurements of dynamical scaling fully confirm this prediction [78]. This result is rather robust. The same type of anomalies have been shown to hold also for quasi-one-dimensional MPC dynamics, namely in the case of a fluid confined in a box with a large enough aspect ratio [18].

F. Anomalous heat transport in material science
The impact of the discovery of anomalous heat transport in anharmonic chains has triggered the search for such an important physical effect in real low-dimensional materials. There is nowadays a vast literature, where these kind of phenomena has been predicted and experimentally observed and part of this growing research field has been sketched in [9]. Here we just want to illustrate two recent contributions which allow the reader to appreciate the relevance of this phenomenon for nanowires and polymers. We want to point out that an overview on the recent literature in both fields is contained in the bibliography of these recent contributions.
In [79] the problem of the lattice thermal conductivity in Si-Ge nanowires has been tackled by solving the Boltzmann transport equation. More precisely, the authors used a Monte Carlo algorithm for sampling the phonon mean free path, combined with phenomenological ingredients concerning a suitable representation of realistic boundary conditions. It is quite remarkable that they find evidence of a heat transport mechanism ruled by a Lévy walk dynamics of phonon flights through the lattice structure. In particular, the phonon mean free paths are found to be characterized by a heavy-tailed distribution, which is associated to an anomalous diffusive behaviour, characterized by a size-divergent heat conductivity κ(L) ∼ L 0.33 . This behavior has been checked for system sizes in the range 10nm < L < 10µm. Note that the phonon mean free path is orders of magnitude smaller than this size range. It is important to point out that this scenario is robust for different alloy compositions, where the Ge component varies in the range [6%, 86%]. All of these results fully agree with the theoretical expectation that anharmonic chains with leading cubic nonlinearity should exhibit a divergent heat conductivity with an exponent γ = 1 3 . In [80] atomistic simulations have been performed for poly (3,4-ethylenedioxythiophene), termed PEDOT, a conjugated polymer which is considered to be of interest in view of its tunable and large electrical conductivity, transparency and air stability [81]. The authors simulated this polymer model in d = 1 and in d = 3, both in equilibrium and nonequilibrium setups. More precisely, equilibrium measurements have been performed by estimating the dependence of the heat conductivity κ on the system size L by the Green-Kubo formula, where one has to estimate the asymptotic decay in time of the correlator of the total energy current (2). The outcome of this analysis has been compared with the numerics obtained in nonequilibrium conditions. The setup adopted in this case is based on a transient measurement of the effective heat diffusivityκ. The two halves of the system are initially prepared in two thermalized states at different temperatures T 1 and T 2 . By running molecular dynamics simulation one can measurē κ as a function of L during the transient evolution to thermal equilibrium state at temperature (T 1 + T 2 )/2. More precisely, the estimate ofκ relies on the fit of the time-dependent temperature difference in the two regions (for details see Eq.(9) in [80]). Finally, the thermal conductivity is obtained by the formula κ =κ ρ c V where ρ is the polymer mass density and c V its specific heat. The authors have obtained consistent results showing that for the polymer chain anomalous diffusion is obeserved, κ(L) ∼ L γ with γ 1 2 , while for the polymer crystal standard diffusion, i.e. a size-independent finite thermal conductivity κ, is recovered. These results are quite remarkable, because they provide a very clear confirmation of the role played by space dimension in determining anomalous transport effects. On the other hand the authors point out that the exponent γ 1 2 does not agree with the expected theoretical value 1 3 , since the phenomenological AMBER nonlinear potential adopted for the PEDOT model is certainly asymmetric, i.e. dominated by a leading cubic nonlinearity. Simulations of the polymer chain have been performed for quite large system sizes, namely 0, 376µm < L < 7.526µm. This notwithstanding, we cannot exclude that, as discussed in the previous section, the combination of finite size effects and nonlinearity might be at work also in this case, yielding a power-law divergence of κ(L) that should be compatible with a symmetric phenomenological potential.
These findings should be compared with those obtained for simpler models: for instance, in the case mentioned in Sec. II.C, the exponent γ = 1 2 has been found in a chain model with a quadratic interaction potential between the beads [60] (notice that the possibility of displacements in both the horizontal and the vertical directions make this model non-integrable, at variance with the harmonic chain). On the other hand, a three-dimensional anharmonic chain with cubic and quartic interaction has been shown to belong instead to the KPZ class with γ = 1/3 [82]. We could speculate that in any polymer model the basic scaling associated to anomalous transport is determined just by the quadratic term, but this is certainly a conjecture that would deserve further theoretical and numerical investigations.
To conclude, we briefly mention some experimental investigations. Thermal properties of nanosized objects have an intrinsic technological interest for nanoscale thermal management. In this general context, nanowires and single-walled nanotubes have been analyzed to look for deviation from the standard Fourier's law [5]. Some experimental evidence of anomalous transport in very long carbon nanotubes has been reported [83] although the results appear controversial [84]. Experiments demonstrating a non-trivial length dependence of thermal conductance for molecular chains have also been reported [85]. Admittedly, such experimental evidences of anomalous transport are rather limited, and not exempt from criticism. Heat transfer measurements on nano-sized objects are notoriously difficult but may possibly be undertaken in the future, possibly guided by the theoretical insight here summarized.

III. COUPLED TRANSPORT
The purpose of this section is to discuss the relevance of anomalous diffusion in coupled transport. In particular, we shall focus on steady-state transport and use for concreteness the language of thermoelectricty [2,86], where the coupled flows are charge and heat flow (other examples where the flow coupled to heat is particle or magnetization flow could be discussed similarly). Moreover, we shall limit our discussion to power production, even though many of the results and open problems highlighted below could be readily extended to refrigeration. Thermoelectricity is a steady state heat engine. Relevant quantities to characterize the performance of a generic heat engine, operating between a hot reservoir at temperature T h and a cold one at T c , are: • The efficiency η = W/Q h , where W is the output work and Q h the heat extracted from the hot reservoir. For cyclic as well as for steady-state heat engines, the Carnot efficiency η C = 1 − T c /T h is an upper bound for the efficiency η.
• The output power P . It is common belief that an engine reaching the Carnot efficiency would require a quasistatic transformation, i.e. infinite cycle-time, implying vanishing power. For steady-state engines this argument is replaced by the one that finite currents would imply dissipation, thus forbidding Carnot efficiency for non-zero power. Hence, it is important to consider the power-efficiency trade-off. This is a key problem in the field of finite-time thermodynamics [87], in relation to the fundamental thermodynamic bounds on the performance of heat engines as well as to the practical purpose of designing engines that, for a given output power, work at the maximum possible efficiency. For classical cyclic heat engines, whose interactions with heat bath can be described by a Markov process, it was proved [88] that the mean power P has an upper bound where A is a system-specific prefactor [see also [89] for an analogous linear-response result within the framework of stochastic thermodynamics [90]]. While at first sight this bound implies that P → 0 when η → η c (and of course when η → 0), so that an engine of that kind with finite power never attains the Carnot efficiency, one cannot exclude that the amplitude A diverges as the efficiency approaches the Carnot value [91].
• The fluctuations in the power output around its mean value P . Indeed, large fluctuations render heat engines unreliable. Especially for heat engines at the nanoscale, one expects power fluctuations due, e.g., to thermal noise, which are not negligible in comparison with the mean output power. In general, one would like to obtain high efficiency (as close as possible to the Carnot efficiency), large power, and small fluctuation. However, a trade-off between these three quantities has been proved [92] for a broad class of steady-state heat engines (including machines described by suitable rate equations or modeled with an overdamped Langevin dynamics): where the (steady-state) power fluctuations are given by where P (t) is the mean delivered power up to time t. Since P (t) converges for t → ∞ to P as 1/ √ t, an additional factor of t in (8) is needed to obtain a finite limit for ∆ P . Eq. (7) tells us that efficiency close to Carnot and high power entail large fluctuations. We note that the bound (8) has been recently generalized to periodically driven systems [93].

A. Linear response
Within linear response, the relationship between currents and generalized forces is linear [94,95]. In particular, for thermoelectric transport we have where j e is the electric current density, j u is the energy current density, and the conjugated generalized forces are F e = −∇(µ/eT ) and F h = ∇(1/T ), where µ is the electrochemical potential and e is the electron charge. The coefficients L ab (a, b = e, u) are known as kinetic or Onsager coefficients; we will denote L the Onsager matrix with matrix elements L ab . Note that the (total) energy current j u = j h + (µ/e)j e is the sum of the heat current j h and the electrochemical potential energy current (µ/e)j e .
The Onsager coefficients must satisfy two fundamental constraints. First, the second law of thermodynamics, that is, the positivity of the entropy production rate,ṡ = F e j e + F u j u ≥ 0, implies Second, for systems with time-reversal symmetry, Onsager derived fundamental relations, known as Onsager reciprocal relations: L eu = L ue . The kinetic coefficients L ab are related to the familiar thermoelectric transport coefficients: the electrical conductivity σ, the thermal conductivity κ, the thermopower (or Seebeck coefficient) S, and the Peltier coefficient Π: For systems with time reversal symmetry, due to the Onsager reciprocal relations Π = T S.
The thermoelectric performance is governed by the thermoelectric figure of merit Thermodynamics imposes a lower bound on the figure of merit: ZT ≥ 0. Moreover, the thermoelectric conversion efficiency is a monotonous increasing function of ZT , with η = 0 at ZT = 0 and η → η C in the limit ZT → ∞. Nowadays, most efficient thermoelectric devices operate at around ZT ≈ 1. On the other hand, it is generally accepted that ZT > 3 − 5 is the target value for efficient, commercially competing, thermoelectric technology. It is a great challenge to increase the thermoelectric efficiency, since the transport coefficients S, σ, κ are generally interdependent. For instance in metals σ and κ are proportional according to the Wiedemann-Franz law, and the thermopower is small, and these facts make metals poor thermoelectric materials. It is therefore of great importance to understand tha physical mechanisms that might allow to independently control the above transport coefficients.

B. Anomalous transport and efficiency
The theoretical discussion of the role of anomalous (thermal) diffusion in thermoelectric transport is based on the Green-Kubo formula. Such a formula expresses the Onsager coefficients in terms of dynamic correlation functions of the corresponding currents, computed at thermodynamic equilibrium. When the current-current correlations j a (0)j b (t) ( . denotes the canonical average at a given temperature T ) do not decay after time averaging, then by definition the corresponding Drude weight is different from zero. Here Ω is the system's volume and Λ the system's size along the direction of the currents. It has been shown [96][97][98][99] that a non-zero Drude weight D ab is a signature of ballistic transport, namely in the thermodynamic limit the corresponding kinetic coefficient L ab diverges linearly with the system size. Non-zero Drude weights can be related to the existence of relevant conserved quantities, which determine a lower bound to D ab [100,101]. By definition, a constant of motion Q is relevant if it is not orthogonal to the currents under consideration, in thermoelectricity j e Q = 0 and j u Q = 0. With regard to thermoelectric efficiency, a theoretical argument [102] predicts that, for systems with a single relevant conserved quantity, as it is the case for non-integrable systems with elastic collisions (momentum-conserving systems), the figure of merit ZT diverges at the thermodynamics limit, so that the Carnot efficiency is attained in that limit. Indeed, for systems where total momentum along the direction of the currents is the only relevant constant of motion, as a consequence of ballistic transport the Onsager coefficients L ab ∝ Λ. Therefore, the electrical current is ballistic, σ ∼ L ee ∼ Λ, while the thermopower is asymptotically size-independent, S ∼ L eu /L ee ∼ Λ 0 . On the other hand, for such systems the ballistic contribution to det L is expected to vanish [102]. Hence the thermal conductivity κ ∼ det L/L ee grows sub-ballistically, κ ∼ Λ γ , with γ < 1. Since σ ∼ Λ and S ∼ Λ 0 , we can conclude that ZT ∼ Λ 1−γ . That is, ZT diverges in the thermodynamic limit Λ → ∞.
This result has been illustrated in several models: in a diatomic chain of hard-point colliding particles [102] (see Fig. 3), in a two-dimensional system [77], with the dynamics simulated by the multiparticle collision dynamics method mentioned in Section II.E [75], and in a one-dimensional gas of particles with screened (nearest neighbour) Coulomb interaction [104]. In all these (classical) models, collisions are elastic and the only relevant constant of motion is the component of momentum along the direction of the charge and heat flows. In the numerical simulations, openings connect the system with two electrochemical reservoirs. The left (L) and right (R) reservoirs are modeled as ideal gases, at temperature T γ and electrochemical potential µ γ (γ = L, R). A stochastic model of the reservoirs [105,106] is used: whenever a particle of the system crosses the opening, which separates the system from the left or right reservoir, it is removed. Particles are injected into the system through the openings, with rates and energy distribution determined by temperature and electrochemical potential (see, e.g., [2]).
We now show that systems with anomalous (thermal) diffusion allow a much better power-efficiency trade-off than achievable by means of non-interacting systems or more generally by any system that can be described by the scattering theory.
We first briefly discuss noninteracting systems. In this case, we can write the charge current following the Landauer-Büttiker scattering theory [107], adapted to classical physics: where f γ ( ) = e −βγ ( −µγ ) is the Maxwell-Boltzmann distribution function for reservoir γ and T ( ) is the transmission probability for a particle with energy to go from one reservoir to the other: 0 ≤ T ( ) ≤ 1. Similarly, we obtain the heat current from reservoir γ as For a given output power P = (∆µ/e)J e (∆µ = µ R − µ L > 0 and we set for the sake of simplicity µ L = 0), the transmission function that maximizes the efficiency of the heat engine, η(P ) = P/J h,L (P, J h,L > 0, T L > T R ) was determined in [103], closely following the method developed for the quantum case in [108,109]. The optimal transmission function is a boxcar function, T ( ) = 1 for 0 < < 1 and T ( ) = 0 otherwise. Here 0 = ∆µ/η C is obtained from the condition f L ( 0 ) = f R ( 0 ), and corresponds to the special value of energy for which the flow of particles from left to right is the same as the flow from right to left. Thus, if particles only flow at energy 0 , the flow can be considered as "reversible", in a thermodynamic sense. The energy 1 and ∆µ are determined numerically in the optimization procedure [103,108,109]. The maximum achievable power is obtained when 1 → ∞: where ∆T = T L −T R , A ≈ 0.0373, and the superscript st reminds us that the obtained results are within the scattering theory approach. At small output power, P/P where B ≈ 0.493. Note that in the limit P → 0 the upper bound on efficiency achieves the Carnot efficiency and the energy window for transmission, δ = 1 − 0 → 0. That is, we recover the celebrated delta-energy filtering mechanism for Carnot efficiency [110][111][112]. Hence, the Carnot limit corresponds to the above mentioned reversible, zero power flow of particles. It is clear that selecting transmission over a small energy window reduces power production. We can then expect that a different mechanism to reach Carnot efficiency might be more favorable for power production. Such expectations are borne out by numerical data for the above described interacting momentum-conserving systems. For these systems, the Carnot efficiency can be reached without delta-energy filtering [113], and the power-efficiency trade-off can be improved. In Fig. 4, we show, for a given ∆T and different system sizes, η/η C as a function of P/P max . These curves have two branches. Indeed, they are obtained by increasing ∆µ from zero, where trivially P = 0, up to the stopping value, where again P = 0. In the latter case, power vanishes because the electrochemical potential difference becomes too high to be overcome by the temperature difference. Power first increases with ∆µ, up to its maximum value P = P max , and then decreases, leading to a two-branch curve. Note that, despite the relatively high value of ∆T /T = 0.2, the numerical results are in rather good agreement with the universal linear curves, which only depends on the figure of merit ZT [103]. Not surprisingly, such agreement improves with increasing the system size, since |∇T | = ∆T /N decreases when N increases. In Fig. 4, we also show the limiting curve corresponding to ZT = ∞, obtained in momentum-conserving models in the thermodynamic limit N → ∞. The upper branch of this curve is the universal linear response upper bound to efficiency for a given power P . For P/P max 1, this bound reads as follows: which is a much less restrictive bound than bound (20), obtained above from scattering theory. Note that, using the linear response result P max ∝ (∆T ) 2 , from (21) we obtain P ∝ ∆T (η C − η). Accordingly, when η ≈ η C ∝ ∆T , we find the same dependence as in bound (6), that was obtained in a rather different context.

C. Open problems
Several open questions about the role of anomalous transport in coupled transport remain, notably: • As discussed in Sec. II.D, heat transport in the presence of long-range interactions has been recently investigated. However, the role of the range of interactions on coupled transport, and in particular on the power-efficiency trade-off, is unknown.
• While momentum-conserving systems greatly improve the power-efficiency trade-off with respect to the noninteracting case, it is not known how they would behave with respect to the bound (7), which simultaneously involves efficiency, power, and fluctuations. That bound was obtained within the framework of stochastic thermodynamics [90], with the transition rates between the system's states which obey the local detailed balance principle, without precise modeling of the underlying particle-particle interactions. On the other hand, in the models described in this short review paper stochasticity is confined to the baths and the internal system's dynamics plays a crucial role. This allowed us to discuss the impact of constants of motion and anomalous transport on the efficiency of heat-to-work conversion and could be relevant also in relation to the bound (7).
• While the results of this section have been corroborated by numerical simulations of several classical oneand two-dimensional systems, their extension to the quantum case remains as a challenging task for future investigations.
• The discussion of this section neglected phonons. However, besides their fundamental interest, the results presented here could be of practical relevance in very clean systems where the elastic mean free path of the conducting particles is much longer than the length scale associated with elastic particle-particle collisions, for instance in high-mobility two-dimensional electron gases at very low temperatures. Phonon-free thermoelectricity (more precisely, thermodiffusion) has been experimentally realized in the context of cold atoms, first for weakly interacting particles [114] and more recently in a regime with strong interactions [115]. In this latter case, a strong violation of the Wiedemann-Franz law has been observed. Such violations could not be explained by the Landauer-Büttiker scattering theory. It would be interesting to investigate whether in such systems, where a high thermoelectric efficiency has been observed, the non-interacting bound on efficiency for a given power could be outperformed.

IV. OVERVIEW
In spite of the significant progress made over the last decades, the study of anomalous heat transport in nonlinear systems still remains a challenging research field. While this review focused on some promising directions in classical systems, a main avenue for future investigations should undoubtedly be sought in the quantum domain. At the quantum level, anomalous heat transport is considerably less understood than in the classical case due to both conceptual and practical difficulties. The same definition of thermodynamic observables like temperature, heat, and work, or of the concept of local equilibrium, is problematic in nanoscale systems. For instance, in solid-state nanodevices we can have structures smaller than the length scale over which electrons relax to a local equilibrium due to electron-electron or electron-phonon interactions. Consequently, quantum interference effects, quantum correlations, and quantum fluctuation effects should be taken into account [2]. In particular, many-body localization provides a mechanism for thermalization to fail in strongly disordered systems, with anomalous transport in the vicinity of the transition between many-body localized and ergodic phases [116,117].
From a practical viewpoint, one has to face the computational complexity of simulating many-body open quantum systems, with the size of the Hilbert space growing exponentially with the number of particles. Notwithstanding these difficulties, a time-dependent density matrix renormalization group method allows the computation of transport properties of integrable and nonintegrable quantum spin chains driven by local (Lindblad) operators acting close to their boundaries [118]. Sizes up to n ∼ 100 spins can be simulated, much larger than the n ∼ 20 spins achievable with other methods like Monte Carlo wavefunction approaches [119]. The obtained results confirm the relevance of constants of motion on transport properties, with integrable systems that exhibit ballistic heat transport, while for quantum chaotic systems heat transport is normal (Fourier's law, see [120,121]). In passing, we note that for magnetization transport in some integrable models like XXZ one can obtain diffusive behavior (Fick's law, see [118]).
From a thermodynamic perspective, the use of local Lindblad operators is problematic. With the exception of quantum chaotic systems, such operators do not drive the system to a grand-canonical state [122]. Furthermore, the use of local Lindblad baths may result in apparent violations of the second law of thermodynamics [123]. Global Lindblad dissipators are free from such problems and can be used to simulate heat transport, see e.g. [124], but are not practical, in that limited to very small system sizes. Furthermore, it is crucial for the description of quantum heat engines in the extreme case where the working medium may even consist of a single two-level quantum system, to take into account medium-reservoir quantum correlations as well as non-Markovian effects, not included in the standard, weak-coupling Lindblad description of quantum open systems. For first steps in this challenging direction, see [125][126][127]. The investigation of anomalous heat transport in such regime is terra incognita.