Heat transport in oscillator chains with long-range interactions
coupled to thermal reservoirs
We investigate thermal conduction in arrays of long-range interacting rotors and Fermi-Pasta-Ulam (FPU) oscillators coupled to two reservoirs at different temperatures. The strength of the interaction between two lattice sites decays as a power of the inverse of their distance. We point out the necessity of distinguishing between energy flows towards/from the reservoirs and those within the system. We show that energy flow between the reservoirs occurs via a direct transfer induced by long-range couplings and a diffusive process through the chain. To this aim, we introduce a decomposition of the steady-state heat current that explicitly accounts for such direct transfer of energy between the reservoir. For , the direct transfer term dominates, meaning that the system can be effectively described as a set of oscillators each interacting with the thermal baths. Also, the heat current exchanged with the reservoirs depends on the size of the thermalised regions: in the case in which such size is proportional to the system size , the stationary current is independent on . For , heat transport mostly occurs through diffusion along the chain: for the rotors transport is normal, while for FPU the data are compatible with an anomalous diffusion, possibly with an -dependent characteristic exponent.
pacs:63.10.+a 05.60.-k 44.10.+i
The combination of various ingredients such as nonlinearity, reduced dimensionality, disorder and topology, may yield quite a complex scenario of transport properties in many-body systems. In a general perspective, statistical mechanics should provide us the tools for investigating this problem when a physical system is driven out of equilibrium by some external non-conservative force or when a gradient is imposed by external reservoirs, exchanging energy, momentum and mass with the system. In particular, one is interested in investigating stationary conducting states, characterized by a minimal entropy production rate, whose physical manifestation is the presence of stationary currents flowing through the system.
Within this vast, interdisciplinary research domain, insight into universal features can be achieved by studying simple paradigmatic models. Among them, systems of classical coupled oscillators are of particular interest, as they represent a large variety of different physical problems, like atomic vibrations in crystals and molecules or field modes in optics and acoustics.
For nonlinear, low-dimensional lattice models there is currently a detailed understanding of their anomalous transport properties LLP03; Basile08; DHARREV; Lepri2016 leading to the breakdown of the classical Fourier law, typically due to a superdiffusive transport mechanism.
On the other hand, transport properties in long-range interacting systems has received some attention only recently avila2015length; Olivares16; bagchi2017thermal. The problem is certainly far from trivial, because it is well known that these systems may relax to long-living metastable states and exhibit anomalous diffusion of energy Bouchet2010; RuffoRev; Campa2014. Unusual effects like the lack of thermalization upon interaction with a single external bath deBuyl2013 or the presence, in isolated systems, of non-isothermal inhomogeneous stationary states where the density and the temperature are anticorrelated 2015PhRvE..92b0101T; gupta2016surprises have been observed. Another unusual feature is the fact that perturbations may propagate with infinite velocities making this class of systems qualitatively different from their short-ranged counterparts Torcini1997; Metivier2014.
What should we expect when two or more external reservoirs drive a long-range interacting system in out-of-equilibrium conditions? Some peculiar features, pointed out in a recent paper by Avila et al. avila2015length have been confirmed and more systematically investigated in a long-range coupled rotor chain Olivares16. This model is quite interesting, because for nearest-neighbor (i.e. fully localized) interactions it exhibits normal transport properties Giardina99; gendelman2000, at variance with most models of nonlinear chains, that have been found to exhibit a diverging heat conductivity with the system size LLP03; Basile08; DHARREV; Lepri2016. The main outcome reported in Olivares16 is that Fourier law is recovered only for sufficiently short-range interactions, while the bulk conductivity seems to vanish in the mean-field limit.
An analogous investigation has been performed for the Fermi-Pasta-Ulam (FPU) -chain with long-range interactions bagchi2017thermal. The main result in this case is that exhibits a power-law divergence with the system size for any value of the long-range exponent. Surprisingly enough the divergence exponent of does not depend monotonously on the range exponent. Moreover, for the range exponent equal to 2, it takes a peculiar value close to 1, which seems to indicate a sort of ballistic transport regime, as observed in integrable models, such as the harmonic chain RLL67, or the Toda lattice Toda79; Kundu2016.
In this paper we analyze similar long-range models and we introduce a suitable setup for clarifying the basic mechanisms underlying heat transport processes (see Sec. II). The proper observables for characterizing the non equilibrium stationary regime in a system with long-range interaction are introduced in Sec. III, while in Sec. IV we present the results obtained for the rotor chain. The interpretation of the transport mechanism is outlined in terms of a two-component stationary heat current (see Sec. V). In Sec. VI we extend the same kind of analysis to the -FPU model. Conclusions and perspectives are summarized in Sec. VII.
Ii The setup
Long-range Hamiltonian models - We consider a chain of anharmonic oscillators whose dynamics is governed by the Hamiltonian
where and are canonically conjugated variables, is a two-body interaction potential and is the mass of each oscillator. In order to have a direct control on the range of the interaction, we specialize to the following form
where is the range parameter and is a generalized Kac prescription, which guarantees the extensivity of Hamiltonian (1) in the thermodynamic limit. We have also assumed that the interaction between two sites and depends only on the relative displacements through the function . For the potential in Eq. (2), is defined as
Note that for (mean field interaction) reproduces the standard Kac prescription, while for (nearest-neighbour interaction) one has . Overall, for any fixed size , is a monotonic decreasing function of , while the interval identifies the parameter region corresponding to non-additive interactions in one dimension RuffoRev.
In this paper we are going to compare two different interaction functions , namely the one of the rotor-chain Tamarit2000 (sometimes called the -XY model or -HMF vanderberg2010)
and the one of the FPU- chain
For the sake of simplicity we have eliminated any explicit dependence on any coupling parameter, that have been set to unit. In fact, this choice simplifies numerical studies, while it is not a prejudice of generality, since physical scales can be recovered by suitable rescaling of time and of the amplitudes of the canonically conjugated variables and . For what concerns boundary conditions, they will be specified throughout the paper. Here we just want to point out that in a long-range interacting system driven in an out-of-equilibrium stationary state the choice of the boundary conditions deserves some attention.
In the limit , the above models reduce to their nearest-neighbor versions, whose transport properties have been studied in great detail in the last two decades, see Refs. Giardina99; gendelman2000; Yang2005; Iubini2016 and LLP97; Lepri05; Wang2011, respectively. For finite , transport properties of the rotor model (4) have been studied recently in Ref. Olivares16. We remark that the FPU model considered in Ref. bagchi2017thermal differs from (5) as the long-range interaction is only in the nonlinear term, while here we consider the case in which the long-range terms are present already in the harmonic limit. Other variants of the model where nonlinearity is local have been considered in the literature Miloshevich2015; Miloshevich2017.
Heat transport in a temperature gradient - Nonequilibrium stationary states (NESS’s) are studied by attaching reservoirs at temperature and to the particles in two different regions of the system, A and B, extending over and lattice sites, respectively (e.g., see Fig 1).
The reservoir dynamics here implemented corresponds to the so-called Maxwellian heat bath LLP03: thermalized system particles interact by a sequence of elastic collisions with particles of mass , selected from a one-dimensional ideal gas in equilibrium at temperature . The corresponding numerical algorithm can be simplified by adopting the condition , so that each collision event reduces to exchanging the momentum of the system particle with the momentum of the selected gas particle. The sequence of time intervals between collisions is obtained by a Poissonian distribution
where is the strength of the coupling with the Maxwellian heat bath. For the numerical simulations presented in this paper, we have chosen and . Denoting with the sequence of the Poisson-distributed time intervals, labeled by the integer index , in region , the -th collision event occurs at time . The same definition holds for collision timing in region , just by changing label with . If the interaction regions with the reservoirs involve more than one lattice site (i.e. or ), independent collisions are generated on each site.
Note that during the time between two consecutive collisions the dynamics of the system is Hamiltonian. We have numerically integrated the equations of motion and using a symplectic 4-th order McLachlan-Atela algorithm mclachlan1992accuracy, with time step . This choice guarantees a sufficiently accurate sampling of trajectories in all the conditions reported throughout this paper. More precisely, the simulation protocol here adopted amounts to the following steps. First, the system is evolved to a thermal equilibrium state at temperature . Then, it is further evolved only in the presence of a temperature gradient, imposed by the heat baths at temperature and , acting in regions and . This evolution lasts over a time interval , long enough for observing a NESS, characterized by a stationary entropy production rate, together with stationary currents and temperature profiles. In general, depends on the system size and on the range parameter : in the whole range of and for the largest systems sizes explored in this paper we have checked that time units is sufficient for obtaining an effective sampling of the NESS’s.
Assuming that local equilibrium conditions set in the NESS, the temperature profile is computed as the time average of the local momentum fluctuations with , i.e.
In order to obtain a clear understanding of energy transport mechanisms in long-range models it is useful introducing suitable definitions of heat currents, that we are going to illustrate in the next section.
Iii heat currents
In the out-of-equilibrium setup sketched in Fig. 1, stationary heat currents depend on the boundary conditions imposed on the thermalized regions and , such as the boundary temperatures and and the bath coupling strength . While this scenario is well understood for short-range systems LLP03, several unusual properties arise when long-range interactions are considered. In order to analyze the NESS generated by a given thermal gradient, it is first necessary to define the heat current flowing from region to region . This task can be accomplished by computing the heat currents, and , exchanged by the system with the reservoirs at temperatures and . This amounts to measure the time average of the variation of kinetic energy of the particles in contact with the reservoirs in regions and . More precisely, we define
where is the variation of local kinetic energy produced by the collision of the Maxwellian heat bath with the particle in at time . The definition of is readily obtained by changing label with in (8).
When a NESS is reached, the energy balance condition reads
In what follows, rather than and separately, we prefer to consider the average total heat current
because it incorporates the fluctuations originated by both heat reservoirs.
In addition to this global observable, it is worth analyzing the average local heat current from site to site , i.e.
where is the force exerted on the particle at site by the one at site . This definition stems from the condition of local energy conservation LLP03. Note that, due to the symmetry of , is antisymmetric under the exchange of indices and .
The current-matrix , whose entries are the average local currents defined by Eq. (11), provides information on the relevant internal channels of the system, through which heat can flow. It should be clearly distinguished from the flow towards the the reservoirs, and and to this aim we use for it a different symbol. As shown in liu2007heat; liu2010heat, the observables defined in Eqs. (8) and (11) have proved to be useful for a clear understanding of heat transport in generic oscillator networks. In what follows, we are going to discuss how these observables depend on the range parameter .
One of the key features displayed by systems with long-range interactions is the lack of additivity RuffoRev: even in the thermodynamic limit, it is impossible to decompose an initial system at thermodynamic equilibrium into two effectively noninteracting subsystems. Equivalently, when one puts in contact two identical and independent long-range subsystems at thermal equilibrium, the whole system may display properties that do not correspond to the original equilibrium state. Therefore, when dealing with thermal conduction problems, a natural question is whether a similar indecomposability holds in the nonequilibrium setup shown in Fig. 1, where one has implicitly identified a bulk system and two thermalized leads.
To better investigate this point, we have considered two different coupling schemes that specify the system-reservoir interaction. In the first scheme (extensive coupling) we introduce a finite fraction of thermalized sites in both regions A and B, i.e. we assume that and increase linearly with the system size . In the second scheme (finite coupling), we assume that and are fixed, independently of .
Iv the rotor chain
In a first series of simulations, we have studied the NESS’s of the rotor chain (4) for different values of the range parameter and for different thermal coupling schemes.
In Fig. 2 we show the stationary temperature profiles for a chain in contact with two reservoirs at temperature and and fixed boundary conditions. The resulting NESS’s have been studied for both the extensive coupling (see panel (a)) and for the finite coupling (see panel (b)). In agreement with previous observations obtained within the finite coupling scheme Olivares16, we find that the temperature profiles are almost flat for , while they are linear for . Moreover, the profiles do not display any relevant dependence on the choice of the coupling scheme.
Quite a different behavior is observed when one analyzes the total average heat current exchanged with the external reservoirs (see Eq. (10)). Let us discuss this point starting from the case , that refers to a rotor chain with nearest-neighbor interactions. In this case, normal heat conduction was observed as a consequence of the periodicity of the interaction potential gendelman2000; Giardina99. In fact, we recover the characteristic diffusive scaling for both choices of the thermal coupling (see the black circles in both panels of Fig. 3). Such a scaling is maintained also when the range parameter is reduced down to , i.e. when additivity is preserved: in both panels of Fig. 3 we have just reported the cases and 2. These results indicate that heat transport occurs essentially via energy diffusion through the bulk of the system for and independently of the extensive- or finite-coupling scheme.
Conversely, for the current displays two different scalings, depending on the choice of the coupling scheme. More precisely, for the extensive-coupling case (see panel (a) in Fig. 3) is essentially constant with , while for the finite-coupling case (see panel (b) in Fig. 3) scales again as . We have checked that this scenario is not altered when dealing with a chain with periodic boundary conditions: in this case the chain turns into a ring of length , where also , and the maximum distance between a pair of particles is .
It is important to point out that both of the situations observed for are originated by the presence of effective “short circuits”, that allow for a direct energy transfer between the thermalized leads. In particular, the -scaling of observed in the finite-coupling case has nothing to do with energy diffusion, as we discuss in Sec.V. An intuitive argument in the mean-field limit is that any pair of thermalized sites and can interact directly through a long-range force . However, it is less evident to figure out how preferred channels of heat transport organize, when mutually interacting oscillators are considered and, more importantly, how they depend on the range parameter . In order to spotlight these peculiar energy pathways, we have computed the stationary current-matrix for different values of . The results are shown in Fig. 4 for a system with sites (periodic boundary conditions) in the presence of the extensive thermal coupling.
It is useful to recall that for a rotor chain with pure nearest-neighbor interactions , the local conservation of energy implies 111Here we assume for simplicity that the sites and are not in contact with the external reservoirs.
where Giardina99; gendelman2000. As a result is a tridiagonal matrix with vanishing elements on the diagonal (data not shown) and it reflects the existence of a channel through which energy is transported, see panel (e). For (panel (a)) maintains essentially the same properties of the nearest-neighbor case, although the nonvanishing elements are now smeared over a broader region around the diagonal. Upon reducing the range parameter from to , this structure progressively disappears in favor of a new pattern which is organized over horizontal and vertical bands, that involve the dominant contribution coming from thermalized sites (panels (b), (c), (d)). For a coexistence of the two patterns is observed, while for the original diffusive pattern discussed in the nearest-neighbor limit has completely disappeared. A closer inspection to panel (d) reveals that the two reservoirs exchange energy either via a direct transfer between two thermalized sites or via a mediated process, that involves a generic site that is not directly thermalized. This mechanism is sketched in panel (f) and provides a simple argument to explain the presence of flat temperature profiles close to , as shown in Fig. 2. In practice, each lattice site in the bulk is effectively coupled both to the hot and to the cold reservoir. As a result, its temperature settles to the average value independently of , so that the average heat current exchanged with any other site in the bulk is practically negligible.
Altogether, by comparing panels (e) and (f) in Fig. 4, we conclude that in the mean-field limit, , the system is organized as an array of oscillators connected in parallel to the reservoirs.
V Decomposition of stationary flow
Before addressing the case of a generic , let us first discuss the two extreme cases of mean-field and short range interaction, respectively. In the former case , the total heat current is constant for the extensive coupling and scales like for the finite coupling. These numerical results can be explained in terms of a two-site diffusion process between a thermalized site and a generic bulk site, whose temperature does not depend on . Indeed, let us evaluate how the current element , with representing a generic site thermalized at (the donor) and a generic site in the bulk of the chain (the acceptor) at temperature , scales with the system size . Since equivalent energy channels are available for the donor site , an energy fluctuation created on it has a probability to be transferred to the acceptor site . Moreover, due to the Kac prescription (see Eq. (3)), the interaction force scales as . As a result, the generic current element scales as . Since in any NESS the total current , defined in Eq. (8), has to be equivalent also to the following expression
we can conclude that is constant for the extensive coupling and for the finite coupling, in agreement with the numerical results.
We point out that the above mechanism does not depend on the details of the interaction potential . Indeed, a similar result is found for the long-range FPU chain (see Sec. VI). In addition, we have verified that the patterns of the current-matrix , obtained for , are stable with the system size .
In the short range limit , the above two-site (donor-acceptor) picture does not apply: for the case of the rotor chain transport is known to be normal Giardina99; gendelman2000 and energy transfer occurs through a standard diffusion process, going through all the sites of the chain. In this case, the nonvanishing current elements are proportional to the global temperature gradient and therefore scales as .
Altogether, these considerations suggest that for any finite , the stationary heat current flowing through the chain is a combination of mean-field-like and local components. This is why it is worth describing the heat transport process in long-range models in terms of a two-component flow. Let us make explicit the dependence of the current-matrix on the range parameter, i.e. . Again we want to focus our attention on the matrix elements , where and .
We assume that for large, we can decompose as follows
where and are summable quantities, namely
We assume also that and are finite quantities also in the thermodynamic limit. Upon this decomposition, the total heat current defined in Eq. (13) can be expressed in the form
where the first term accounts for the long-range contribution to heat transport, while the second one accounts for the short-range diffusion process. Note that the factor in front of the first addendum is a consequence of the assumption that the quantities and are independent of , i.e. the long-range contribution from the thermalized sites in region is the same from all of these sites.
The two contributions in Eq. (14) for different values of the range parameter are highlighted in Fig. 5, that has been obtained for chains with periodic boundary conditions. In particular, we show the current profiles obtained from numerical simulations of the rotor chain with periodic boundary conditions and extensive coupling, . In this setup, Eq. (17) predicts that long-range contributions to the heat current are constant, while short-range ones scale as . In order to separate visually the thermalized sites from the bulk sites, in these pictures we have chosen a spatial reference frame where the site labels the first particle which is not in contact with the heat bath at temperature . Accordingly, the thermalized sites, , have a negative space coordinate for any . For (see panel (a) in Fig. 5) we observe that the rescaled observable is essentially independent of , thus proving that the only relevant term in Eq. (14) is due to , while has to vanish. The flat profile indicates that each thermalized site transfers energy uniformly to all the other sites of the bulk. For (panel (b)), the rescaled observable exhibits a good data collapse in a macroscopic region region around , where the reservoir at temperature is located. This indicates that the relevant term for large in Eq. (14) is due to . On the other hand, in a limited region close to we observe a significative deviation from the data collapse (see the upper inset in panel (b)). In fact, in the lower inset of panel (b) we show that the rescaled quantity exhibits a definitely better data collapse, thus indicating that the main contribution in this region is due to . In summary, for we are facing a situation where the two contributions appearing in Eq. (14) do not vanish: the first one determines the scaling behavior of in the large limit (i.e., it is constant), while the other simply affects corrections to scaling in a limited region around . Finally, for (see panel (c)), we observe that a a good data collapse is obtain for the rescaled observable : provides the main contribution to Eq. (14) relevant, while the long-range component vanishes. Moreover, the quantity is found to decay quite rapidly as a function of the distance from the reservoir at . From a fit obtained for the chain with largest size analyzed in this paper (), we can estimate an asymptotic power-law decay as with (see the inset in panel (c)). This result confirms the summability of over , as required for the ansatz (14). We have also checked the summability of over the reservoir indices (data not reported).
We can conclude that this phenomenological description of the scaling of the total heat current with the system size in terms of a two-component stationary flow is able to capture the basic mechanisms that coexist and eventually prevail onto each other in the long-range rotor chain. Pictorially, the overall information is summarized in the form of the current matrix , shown in the upper panels of Fig. 4.
Vi the FPU chain
Here we extend our study of heat transport in long-range models to the FPU- chain, whose interaction potential is given in Eq. (5). We report just the main results, without going through a detailed description as for the rotor chain.
A first series of numerical simulations have been performed with the nonequilibrium setup described in Sec. II with and . Fig. 6 shows the scaling of the stationary total heat current as a function of the system size for the extensive thermal coupling and fixed boundary conditions. We have verified that for (i.e., nearest-neighbour interaction) numerics essentially agrees with the expected anomalous scaling consistent with previous studies of the short-range case LLP2003; Pereverzev2003; Lukkarinen2008; Wang2011. The bending of the black circles with respect to the dashed line in Fig. 6 is due to finite size effects. We want to stress that this scenario is different from the rotor chain, which exhibits normal conductivity, i.e. , in the nearest-neighbour limit.
On the other hand, for also for the FPU- chain we have evidence that does not depend on the system size , as shown in Fig. 6 for the cases and 1. The only difference with respect to the rotor chain is that the value of decreases more rapidly with . This result indicates that the basic mechanism of heat transport in this non-additive regime is the same “parallelization” process of heat channels discussed in Fig. 4 for the rotor chain. Indeed, the current-matrix , ( see panels (b),(c) and (d) in Fig. 7) exhibits patterns that are very similar to those shown in Fig. 4, with a fully parallelized structure for .
The scenario for is definitely richer than in the rotor chain. The dependence of on for various values of in this range is reported in Fig. 8. On the basis of these results we can speculate that the closer is to 1, the more finite size effects prevent the possibility of recovering the anomalous scaling of the nearest-neighbour case. In fact, all curves in this Figure exhibit a clear tendency to bend for increasing : only for the largest value of reported in this Figure, namely , we can observe an approach to the expected anomalous scaling for sufficiently large values of . Moreover, a study of the flux elements performed for the system sizes and suggests that a decomposition similar to Eq. (14) holds also for the FPU chain, provided that the diffusive term is replaced by an anomalous one of the form (data not shown). We want to point out that a deeper inspection of this regime is practically demanding, because of limitations of the available numerical resources. Anyway, we should conclude that for increasing values of there is a monotonous tendency of to recover the expected anomalous scaling over a larger and larger domain of . The results reported in Fig. 6 for (see the red squares) indicate that this is not the case: seems to be independent of . We can exclude that this result can be attributed to the parallelization process that sets in for . Indeed, the inspection of the current-matrix reveals that, for , the relevant energy channels are exclusively those around the diagonal, as shown in panel (a) of Fig. 7. At first sight, it seems that in this peculiar case we are facing a sort of ballistic regime, characterisitc of integrable models, e.g. the harmonic oscillator chain or the Toda lattice hu2000heat; LLP03. We point out that the same observation has been reported very recently by Bagchi bagchi2017thermal for a FPU- model with tunable range parameter, where only the nonlinear quartic term was put in the form of Eq. (2), while the harmonic interaction was maintained in the nearest-neighbor form. In Ref. bagchi2017thermal, the apparent ballistic transport for was attributed to the emergence of a quasi-integrable dynamics. On the other hand, recovering the same kind of peculiar behavior for in our version of the long-range FPU- chain testifies at the robustness of this effect, but makes doubtful the possibility of associating it to integrability, or quasi-integrability, whatever this means.
A better understanding of this phenomenon can be obtained by comparing the temperature profile of the stationary state for with those of the harmonic chain and of the Toda lattice. For what concerns the harmonic chain, the temperature profile can be computed analytically RLL67 and it is known to exhibit for increasing values of a flatter and flatter shape far from the boundaries around a temperature . In the absence of an analytic estimate, numerics shows that a similar shape is obtained for the Toda chain. In the inset of Fig. 6 we can observe that for model (5) with the temperature profile deviates more from flatness the larger the system size . On the basis of this observation, we can just argue that finite size effects are playing a dominant role in this case. This said, it seems quite hard to provide any convincing argument that could explain what happens for . Preliminary results (not reported here) for long-range models with different kinds of FPU-like interactions confirm the peculiarity of this case together with the important role played by finite size effects.
Heat transport in long-range models exhibits quite interesting and peculiar features, that depend on the nature of the interaction potential. The comparison between the rotor chain and the FPU- chain confirms this scenario, that was already outlined in some recent pubblications avila2015length; Olivares16; bagchi2017thermal. In this paper we have shown that in the non-additive regime corresponding to the mechanism of heat transport is dominated by the direct, mean-field like contribution, irrespectively of the kind of the interaction potential.
To rationalize the numerical results, we have proposed a decomposition of the steady-state flow, Eq.(14), that accounts for the direct energy transfer between the heat baths. Such decomposition allows us to conclude that for heat transport is dominated by local interaction mechanisms, that are expected to reproduce the same scaling properties of the total heat flux with the system size , observed in the nearest-neighbor limit . On the other hand, finite size effects have been found to have a strong influence in the study of heat transport in long-range models for any finite value of . This is the main reason why clear-cut conclusions are very hard to be drawn just relying upon accessible numerical studies. This is the case also when dealing with the peculiar situation observed in the FPU- chain for . In this case we are facing quite an unexpected behaviour, that seems to be dominated by a ballistic transport mechanism, rather than reproducing the anomalous scaling of , typical of the nearest-neighbour limit. The fact that such a peculiarity of the case seems to emerge in different versions of long-range FPU-like models with an ubiquitous manifestation of subtle finite-size effects testifies at the interest of this phenomenon, that is still waiting for a convincing explanation.