## Abstract

We consider non-equilibrium steady-state situations for thermoelectric systems with non-local and non-linear effects. We show that the Onsager symmetry relations for effective transport coefficients break down. We also estimate the consequences of such a breakdown for the efficiency of the thermoelectric energy conversion which, under some conditions, could be higher than in the usual linear regime with Onsager symmetry.

## 1. Introduction

Coupled transport processes find a natural theoretical framework in classical thermodynamics of irreversible processes [1–3]. Among the different coupled phenomena, thermoelectric effects (i.e. the possibility of coupling thermal transport and electric current) are currently attracting much attention, since they offer a promising avenue in energy management [4–6].

The analysis of thermoelectric effects rests on the so-called *second Kelvin relation* *Π*=*ϵT*, *Π* being the Peltier coefficient, *ϵ* the Seebeck coefficient and *T* the temperature. This relation, stated by Lord Kelvin in 1854, expresses a subtle and fundamental connection between the Peltier effect (namely, the heating, or cooling, of an electrified junction) and the Seebeck effect (i.e. the conversion of temperature differences into electricity). On thermodynamic grounds, proof of the validity of that relation may be given by means of the Onsager relations (ORs) [7,8], which state the symmetry of the transport coefficients matrix and link the thermodynamic fluxes with the conjugated thermodynamic forces.

Although one may find several examples in which the ORs comply with the experimental evidence [2,3,9,10], in some situations, especially in the nonlinear regime, the ORs are no longer valid [11–13].

It is worth observing that the problem of checking the validity of the ORs goes beyond pure theoretical interests, as these relations could also have consequences in practical applications. For example, in maximizing the efficiency of thermoelectric energy conversion, it can be proved that, when the Onsager symmetry (OS) holds, the maximum efficiency depends only on the so-called figure of merit, but a more general expression is required if such symmetry is broken.

Nowadays, some generalizations of the constitutive equations describing thermoelectric effects are being explored in search of new possible strategies to improve the efficiency of thermoelectric energy conversion [13]. These generalizations are necessary when dealing with nanosystems, which are viewed as relevant candidates to improve thermoelectric efficiency. In nanosystems, special attention has to be paid to non-local and nonlinear effects [14–18], since even a small difference in temperature, or electric potential, over a very short length may generate very high gradients. Thus, in practical applications at the nanoscale, it is only natural to wonder whether the ORs still hold, or whether they suffer a break [11,19]. This is the aim of the present investigation.

The layout of the paper is as follows. In §2, we briefly review a way of obtaining the classical second Kelvin relation in the linear regime. In §3, we prove the breakdown of the OS in a mesoscopic nonlinear model of thermoelectric coupling in nanosystems. In §4, we evaluate how a breakdown of the classical relation involving the Seebeck and Peltier coefficients may influence the efficiency of thermoelectric energy conversion. Final comments are given in §5.

## 2. Thermoelectric equations in the linear regime

Consider a thermodynamic model of a rigid heat conductor accounting for both thermal and electric effects. In this case, suitable unknown fields are the internal energy *u* per unit mass and the number of electrons *z* per unit mass [5]. These quantities are governed by the local balance equations of energy and electric charge, namely,
*ρ* as the mass density, ** q** as the local heat flux,

**as the electric field, and**

*E***as the electric current density. As in this case the constitutive equation for the specific entropy**

*i**s*depends on both

*u*and

*z*, by introducing the chemical potential of the electrons

*μ*

_{e}=−(

*z*

_{e}

*T*)∂

*s*/∂

*z*, where

*z*

_{e}is the electric charge per unit mass of the material (i.e.

*z*=

*z*

_{e}

*c*

_{e}, and

*c*

_{e}is the mass fraction of electrons), and the absolute temperature

*T*

^{−1}=∂

*s*/∂

*u*, it is easy to obtain the well-known Gibbs relation

The combination of equations (2.1) and (2.2) gives

Recalling that the usual form of the local balance of the entropy is
*J*^{(s)} is the entropy flux and *σ*^{(s)} is the entropy production, from equation (2.3) we infer that the constitutive equations for *J*^{(s)} and *σ*^{(s)} are, respectively, given by

In particular, equation (2.5b) shows that the entropy-source density *Tσ*^{(s)} consists of a sum of products of the so-called thermodynamic fluxes *J*^{(α)} and thermodynamic forces *X*^{(α)} [1–3], namely,

The simplest way to avoid any possible violation of the second law of thermodynamics (stating that *σ*^{(s)} has to be non-negative along any admissible thermodynamic process) is to assume that each *J*^{(α)} depends linearly on the whole set of *X*^{(α)} through suitable phenomenological coefficients [1–3] satisfying thermodynamic restrictions.

Once coupled with equation (2.6), this observation, also suggested by experimental evidence, for isotropic systems leads to the generalized thermoelectric equations
*L*_{11}, *L*_{12}, *L*_{21} and *L*_{22} can be identified by simple comparisons. In fact, equations (2.7) reduce to the classical thermoelectric equations
*σ*_{e} the electric conductivity, once one identifies *L*_{11}=λ*T*+*σ*_{e}*ϵΠT*, *L*_{12}=−*σ*_{e}*Π*, *L*_{21}=−*ϵσ*_{e}*T* and *L*_{22}=*σ*_{e}. If, by virtue of the ORs [7,8], in equations (2.7) one assumes *L*_{12}=*L*_{21}, the relation *Π*=*ϵT* directly follows from the identifications above.

Indeed, in the literature one can find several proposals for a generalization of the ORs in the nonlinear domain. In particular, we recall the proposal of Hurley and Garrod [20] through statistical methods, as well as that of Verhás [21] with pure macroscopic methods.

The generalizations proposed by Li [22] and Gyarmati [23,24] are also very interesting. These authors suppose that the thermodynamic fluxes *J*^{(α)} arise from a dissipation potential *Ξ*(*X*^{(α)}), in such a way that *J*^{(α)}=∂*Ξ*/∂*X*^{(α)}, and the equality of the mixed second-order derivatives of *Ξ* with respect to the thermodynamic forces leads to
*L*_{αβ} is symmetric, and the ORs hold. The dissipation potential *Ξ* is a convex functional in a neighbourhood of *X*^{(α)}≡**0**, takes a minimum at *X*^{(α)}≡**0**, and it is such that *Ξ*(**0**)≡0. Its Taylor expansion in the neighbourhood of *X*^{(α)}≡**0**, up to the third-order approximation yields
*C*_{αβγ} is such that
*Ξ* with respect to *X*^{(α)}, and of the convexity of *Ξ* in a neighbourhood of *X*^{(α)}=**0**, *C*_{αβγ} is symmetric and *L*_{αβ}(**0**) is symmetric and positive semidefinite.

It is easily seen that the linear theory is recovered by neglecting the third-order terms in equation (2.9). In such a case, both the OS and the condition *σ*^{(s)}≥0 are guaranteed. If, instead, the third-order terms in equation (2.9) are not negligible, then the additional transport matrix *a priori*, but it depends on the forces, i.e. on the thermodynamic process.

The considerations above allow us to conclude that, whenever a dissipation potential exists, the ORs hold, in both the linear and nonlinear cases.

If a dissipation potential does not exist, in the genuinely nonlinear case the matrix of the transport coefficients depends on the thermodynamic forces. In such a case, the OS is conserved if one considers the Taylor expansion of the transport coefficients up to arbitrary order of approximation but, as will be shown in the next section, it can break down in the general case, i.e. for arbitrary forms of the transport matrix.

As far as equations (2.7) are concerned, for negligible values of *μ*_{e}/*z*_{e}, by taking the Taylor expansion of the transport matrix, up to the first-order approximation the procedure would lead to the following phenomenological equations [25,13]:
*α*,*β*=1,2) are constant second-rank tensors, and *α*,*β*,*γ*=1,2) are constant third-rank tensors. Moreover, in the equations above the symbol ⊗ means the product between a third-rank tensor and a vector (i.e. *A*_{ijk}*a*_{k}) and ⊙ is the usual product between a matrix and a vector, (i.e. *B*_{ij}*b*_{j}).

In such a case, the transport matrices entering the nonlinear terms in equations (2.10), i.e.
*T*/*T* and ** E**. The admissible processes are those, and only those, in which the tensors above are positive semidefinite. As a consequence, the second law of thermodynamics becomes a restriction on the thermodynamic processes and not on the constitutive equations.

In this paper, we exclude the existence of a dissipation potential, and do not approximate the transport matrices by a Taylor series. Under such hypotheses, we analyse the consequences of a possible nonlinearity of the phenomenological laws caused by the nonlinearity of the system of equations governing the evolution of the dissipative fluxes. We find that the most important consequences are both the breakdown of the OS and a drastic modification of the efficiency of thermoelectric energy conversion.

To achieve that task, we start from enhanced thermoelectric equations and investigate their implications at the nanoscale. In this sense, the result we find for the breakdown of the OS is less general with respect to other theoretical proposals, but it has the advantage of being more pragmatic, as it is strictly related to practical applications.

## 3. Non-local and nonlinear thermoelectric equations and breakdown of the Onsager symmetry

To cope with nonlinear effects, in recent papers [26–28] nonlinear thermoelectric equations have been introduced on a mesoscopic level. Whenever terms proportional to ∇(*μ*_{e}/*z*_{e}) are neglected, those equations read
**q**^{(p)} and **q**^{(e)} are the partial heat fluxes due to phonons and electrons [26–28], respectively, *c*_{v} is the specific heat, λ_{p} and λ_{e} are the thermal conductivities due to phonons and electrons, respectively, and *τ*_{p}, *τ*_{e} and *τ*_{i} are the relaxation times, which play an important role in the model described by equations (3.1), since they introduce not only relaxational effects but also nonlinear effects. In equations (3.1), all the thermophysical quantities are supposed to be constant, so that the nonlinearity is only due to the product of the partial heat fluxes and of the electric current with their gradients. However, in the most general case the material functions are temperature dependent, and this fact introduces a further nonlinearity in the system of equations. Here, such a situation is not considered.

The thermodynamic model described by equations (3.1) is based on the hypothesis that **q** may be split into the two different contributions **q**^{(p)} and **q**^{(e)}. Then, equation (2.1a) is the sum of the following two different evolution equations:
*u*^{(p)} and *u*^{(e)} are the partial internal energies, such that *u*=*u*^{(p)}+*u*^{(e)}.

In steady states, and whenever **q** and ** i** vary only along the longitudinal direction

*ζ*of the thermoelectric device, from equation (2.1b) one has ∇

_{ζ}

*i*=0 (where ∇

_{ζ}is the longitudinal component of the ∇ operator along

*ζ*), whereas from equations (3.2) it follows that ∇

_{ζ}

*q*

^{(p)}=0, and ∇

_{ζ}

*q*

^{(e)}=

*Ei*. Once these relations have been taken into account, it is easy to see that equations (3.1) read

Note that in the particular physical situation considered here the flux of phonons is ruled by the classical Fourier law. Thus, the phonon relaxation time does not appear in the equations above and so it cannot influence the form of the second Kelvin relation.

For practical applications, in equations (3.3), we introduce the following effective coefficients:
*Ω* is the cross section of the system, and the superscripts (T) and (E) mean the values of the corresponding fields when thermal or electric effects are neglected, respectively, that is, when ∇_{ζ}*T*=0, or *E*=0, in equations (3.3). Moreover, in equations (3.4) Δ*T*/*L* is the constant value of ∇_{ζ}*T* measured on the corresponding cross section. Then, if we suppose that each field *q*^{(E)}, *i*^{(E)}, *E*^{(T)} and *q*^{(T)} assumes a constant value in each cross section, the coupling of equations (3.3) and (3.4) yields

Equation (3.5d) shows a non-classical part of the effective Peltier coefficient related to nonlinear effects due to the coupling between **q**^{(e)} and ** i** in equations (3.2b) and (3.3b) and suggests that at nanoscale the effective Peltier coefficient, which is related to the experimental measurements, may be different from the bulk value. Since a similar coupling is not present in equation (3.2a), and, moreover,

*τ*

_{p}is lacking in equation (3.3a), no physical quantities related only to phonons may contribute to determine

*Π*

^{eff}. From the physical point of view, the fact that

*Π*

^{eff}only depends on

*τ*

_{e}seems rather natural, since the effective Peltier coefficient has been calculated by considering the heat flux when the thermal drag due to the gradient of temperature is neglected; namely, when the only dragging force is the electric field, which, of course, cannot drag the phonons.

Equations (3.5) also allow us to point out the influence of nonlinear effects on the breakdown of the OS. In fact, if one assumes that *Π*−*ϵT*=0 (i.e. ORs holding for the bulk coefficients), then from equations (3.5c) and (3.5d) it follows that
*x**=1 (i.e. when the term *τ*_{e}*Ei*/(*c*_{v}*T*) vanishes), the usual ORs are recovered. Thus, at a given temperature, the degree of violation of the OS depends on the material functions *τ*_{e} and *c*_{v}, as well as on the intensity of the electric field and of the electric current; namely, the higher the product *Ei*, the higher the deviation *Π*^{eff}−*ϵ*^{eff}*T* from zero.

## 4. Efficiency of thermoelectric energy conversion with broken Onsager relations

From the theoretical point of view, in the literature several criteria have been introduced to optimize heat transport [11,29]. In a simple and practical application, here we are faced with an analogous problem; namely, we focus our attention on the optimization of the thermoelectric energy conversion. The generation of electric energy from thermal energy by thermoelectric devices, in fact, is one of the main goals of future energy management. Many research groups are focusing their attention on the search for different strategies to enhance the efficiency in thermoelectric energy conversion. There are many reports in the literature of engineering studies that are seeking the most suitable structure for a thermoelectric device to achieve better thermoelectric conversion of energy. Here we analyse the consequences of accounting for the breakdown of the ORs on the efficiency *η* in the thermoelectric energy conversion, defined as
*P*_{el} being the electric-power output, and *Q*_{tot} being the total heat supplied per unit time. In particular, here we account for the possibility of a breakdown of the ORs by assuming that the second Kelvin relation holds for the bulk values of *ϵ* and *Π* (i.e. *Π*=*ϵT*), but not necessarily for their effective values *ϵ*^{eff} and *Π*^{eff} (i.e. *Π*^{eff}≠*ϵ*^{eff}*T*). In fact, we observe that, at the nanoscale, the material functions may display values which are rather different with respect to the corresponding bulk values on some occasions [27].

Along with equation (3.7), here we assume that the breakdown of the OS is endorsed by the non-dimensional parameter *x**=*Π*^{eff}/(*ϵ*^{eff}*T*).

To derive the expression for thermoelectric efficiency, we use the simplest and most common situation, i.e. we consider a single one-dimensional thermoelectric element of length *L* at steady state, the sides of which are kept at two different temperatures, *T*_{h} (the hottest temperature) and *T*_{c} (the coldest temperature). We assume that both an electric current and a quantity of heat per unit time enter uniformly into the hot side of the element. Recalling that in estimating *P*_{el} one has to take into account both the Seebeck contribution and the Joule dissipation [5], in a first-order approximation we have
*σ*^{eff}_{e} stands for the effective value of the electric conductivity.

The total heat supplied per unit time *Q*_{tot}, instead, is the sum of the purely conductive contribution (given by Fourier's Law) and of the Peltier one [5], that is,
^{eff} is an effective value of the total thermal conductivity of the device.

Introducing equations (4.2) and (4.3) into equation (4.1), by straightforward calculations, we have
*η*_{C}=1−*T*_{c}/*T*_{h} is the usual Carnot efficiency, and
*η*_{C}. On the other hand, from equation (4.4) it is clear that the higher *η*_{r}, the higher *η*.

At this stage, it is convenient to express equation (4.5) as a function of the ratio *ξ*=*iL*/ [λ^{eff}(*T*_{h}−*T*_{c})]. In this way, it becomes
*x**=1. Note that in deriving equation (4.9) we used the relation *x**=*Π*^{eff}/*Π*, the validity of which is stated by equations (3.5*c*,*d*) and (3.7).

From equation (4.9), we recover the well-known result that the higher *Z*^{eff}*T*, the higher *η*. In fact, the search for materials with high values of *Z* is a major topic of research [5,27,30], and the figure of merit of many materials has been explored [31,32].

Indeed, equation (4.9) suggests another way to reach higher values of *η*. In fact, for a given value of the figure of merit, from this equation we also infer that the smaller *x**, the higher *η*_{max}. This means that the thermoelectric efficiency also depends on the degree of the breakdown of the second Kelvin relation for the effective Seebeck and Peltier coefficients. This result is illustrated in figure 1, wherein we plot the thermoelectric efficiency as a function of *x**, for several different values of *Z*^{eff}*T*. It is worth noticing that the effects of the breakdown of the OS on the efficiency is not negligible, and it is higher for higher values of *Z*^{eff}. Thus, in order to enhance the thermoelectric efficiency, it would be desirable to act on both *Z*^{eff}*T* and *x**, making the first as high as possible and the latter as small as possible. This last task can be achieved for moderate values of the product *Ei* and using materials for which the ratio 8*τ*_{e}/(*c*_{v}*T*) is as small as possible, in such a way that the values of the effective and of the bulk Peltier coefficients are sensibly close. This could be important in practical applications, wherein the principal strategy followed by many research groups is to find new materials with high values of the figure of merit by assuming linear constitutive equations.

Since in the case of the model equations (3.1) the non-dimensional parameter *x** is greater than 1, according to figure 1 the breakdown of the ORs arising from nonlinear effects reduces the maximum efficiency.

A result analogous to equation (4.9) was found by Benenti *et al.* [12], under the assumption that the time-reversal symmetry breaks down. In that paper, the authors wrote the Onsager matrix assuming that the material functions depend on the external magnetic field **B**. They assumed that *ϵ*(**B**)≠*ϵ*(−**B**), in contrast to the equality required by Onsager–Casimir symmetry relations [33], whereas λ(**B**)=λ(−**B**) and *σ*_{e}(**B**)=*σ*_{e}(−**B**). The maximum value of the thermoelectric efficiency was calculated as a function of the non-dimensional ratio
*x*=1/*x** and *y*=*Z*^{eff}*x**. However, our analysis is more general, because it does not assume that the breakdown of the OS is related to the magnetic field, but to nonlinear effects, as we will see below.

## 5. Conclusion

In the early 1930s, Onsager proved the symmetry (or reciprocity) relations [7,8], in the absence of a magnetic field. Later on, his conclusion was generalized by Casimir [33] in the presence of a magnetic field. Both the proofs are based on the so-called *principle of microscopic reversibility*, and on suitable hypotheses on the behaviour of the constitutive functions under *time reversal*. The proofs combine the mathematical analysis with *a priori* assumptions based on the physical intuition of the authors.

Some irreversible processes must be described by nonlinear phenomenological laws. This is evident especially at the nanoscale, where the thermodynamic forces may be very high. Thus, in these situations, it is only natural to wonder whether the ORs still hold or whether they suffer a break, and what consequences this breakdown may have for the efficiency of thermoelectric conversion.

In this paper, we have presented a thermodynamic model predicting that the second Kelvin relation for the effective transport coefficients no longer holds at the nanoscale. Indeed, this seems rather natural when dealing with a genuinely nonlinear model, so that the validity of the classical linear relation between the thermodynamic forces and their conjugated thermodynamic fluxes could be questioned. Our proposal here is that, in nonlinear heat conduction at the nanoscale, the relation (3.6), namely
**q**^{(p)}, **q**^{(e)} and ** i**, which are able to comply with nonlinear effects at the nanoscale, in the absence of a dissipation potential.

We have further shown that the search for the limit of validity of the ORs, which at first sight may have only an academic motivation, yields indeed also possible consequences for practical applications. Our theoretical prediction (4.9), in fact, suggests that to enhance the thermoelectric energy-conversion efficiency, one should pay attention not only to the figure of merit (which remains the main factor responsible for high thermoelectric performances), but also to the degree of the breakdown of the second Kelvin relation, as shown by figure 1. The main parameter accounting for the aforementioned breakdown is the non-dimensional ratio *x**=*Π*^{eff}/(*ϵ*^{eff}*T*), whose form is strictly related to the model equations. In the case of the model equations (3.1), it is given by equation (3.7), i.e. it is related to well-known physical parameters. As equation (3.7) predicts that *x** is greater than 1, we concluded that the breakdown of the ORs, following from the effects considered in this paper, yields a reduction in the maximum efficiency of the thermoelectric energy conversion. However, this does not exclude that other kinds of nonlinear effects (such as those expressed in equations (2.10)) could lead to a higher efficiency. It would also be interesting in future research to investigate the possible breakdown of the second Kelvin relation if a different heat conduction theory, such as thermomass theory [34,35], is applied.

We conclude with more about the model equations (3.1). They represent a simplified version of more general model equations, which have the form [28]
_{p}, ℓ_{e} and ℓ_{i} are the mean-free paths of phonons, electrons and of the electric-current density, respectively, while *T*^{(p)} and *T*^{(e)} are the phonon temperature and the electron temperature. The second-order spatial derivatives in equations (5.1) account for non-local effects, which are important to obtain the size dependency of the figure of merit [27]. The model with two temperatures rests on the physical evidence that in thermoelectric coupling it is possible to run into situations characterized by different values of the phonon and electron temperatures [28]. In fact, as the electron mean-free path ℓ_{e} is usually shorter than the phonon mean-free path ℓ_{p}, in the nanowire considered here, when the longitudinal distance *ζ* between heat carriers is such that ℓ_{e}<*ζ*<ℓ_{p}, it is expected, for example, that there will be a very high number of electron collisions, and only scant phonon collisions. This yields that the electron temperature may reach its local equilibrium value, whereas the phonon temperature is still far from its own local equilibrium value. Conversely, whenever the electron mean-free path corresponding to the electron–phonon collisions is large, one may have the so-called ‘hot electrons’, namely, a population of electrons whose average kinetic energy (i.e. the kinetic temperature) is considerably higher than that of the phonons.

The compatibility of equations (5.1) with the second law of thermodynamics has been discussed by Jou *et al*. [28] by means of the classical Coleman–Noll procedure [36] in such a way that necessary and sufficient conditions for the physical compatibility are fulfilled. They also derived both the form of the specific entropy, and that of the entropy flux, which, in the case of a single temperature, respectively, read as follows:

Equation (5.2a), which is in accordance with the principle of maximum entropy at equilibrium, shows that the specific entropy is quadratic in the dissipative fluxes. Such a property makes it invariant with respect to the reversible time reversal. As a consequence of the thermodynamic relation *s*=−∂*ψ*/∂*T*, the same property also holds for the Helmholtz free energy *ψ*=*u*−*sT*. On the other hand, by observing that −*u*=−*ψ*+*T*∂*ψ*/∂*T*, we may conclude that −*ψ* and −*u* are related by a Legendre transformation, with −*s*=∂*ψ*/∂*T* playing the role of the conjugated moment derived by −*ψ*. This mechanical analogy, together with the invariance of *s* under time reversal, leads to the conclusion that, analogous to the mechanical case, in our model the thermodynamic potentials do not change under time reversible time evolution.

Equation (5.2b), instead, points out that *J*^{(s)} contains also extra-flux terms. Jou *et al.* [28] also observed that non-local and nonlinear terms in equations (5.1) change the entropy-source density in
*f*(*Σ*) is a well-defined function of the state space *Σ*, that is, in equation (5.3) there are no uncontrollable quantities.

In the case of a single temperature, and when the nonlinear terms are negligible in the system above, taking into account the interactions between the different heat carriers and the walls, in Sellitto *et al.* [27] a possible size-dependent enhancement of the figure of merit has been discussed. Here, instead, still for a system with a single temperature, we have supposed that the non-local effects are negligible, emphasizing the role of nonlinear effects.

In the future, it would also be interesting to analyse the consequences of accounting for two different temperatures for the heat carriers on the figure of merit.

Finally, let us observe that, if the relaxation effects are negligible, equations (3.1) governing the evolution of the non-equilibrium fluxes **q**^{(p)}, **q**^{(e)} and ** i** reduce to equations (2.8), which are the classical constitutive equations of thermoelectricity for the quasi-stationary fluxes

**q**and

**. Note that the partial heat fluxes are related to the total heat flux by the relation**

*i***q**=

**q**

^{(p)}+

**q**

^{(e)}, while the partial and total heat conductivities satisfy the relation λ=λ

_{p}+λ

_{e}. In such a case, the validity of the second Kelvin relation, which has been confirmed by plenty of experimental results over the last century, proves that the OS holds. From that, we infer that a solution for the entropy inequality in terms of dissipation potential should exist.

## Funding statement

V.A.C. thanks the University of Basilicata for financial support and for funding the research subject in mathematical physics Equazioni costitutive per la conduzione del calore nei nanosistemi.

## Acknowledgements

V.A.C. thanks the Italian Gruppo Nazionale per la Fisica Matematica - GNFM. D.J. acknowledges the financial support from the Dirección General de Investigación of the Spanish Ministry of Economics and Competitiveness under grant FIS no. 2012-33099, the Consolider Project NanoTherm (grant no. CSD-2010-00044) and the Direcció General de Recerca of the Generalitat of Catalonia under grant no. 2009-SGR-00164. A.S. acknowledges both the University of Basilicata for funding the research project Modeling heat and electric transport in nanosystems in the presence of memory, nonlocal and nonlinear effects, and the Italian Gruppo Nazionale per la Fisica Matematica - GNFM for financial support under grant Progetto Giovani 2012.

- Received April 2, 2014.
- Accepted July 10, 2014.

- © 2014 The Author(s) Published by the Royal Society. All rights reserved.