Thermal rectification in mass-graded next-nearest-neighbor Fermi-Pasta-Ulam lattices

We study the thermal rectification efficiency, i.e., quantification of asymmetric heat flow, of a one-dimensional mass-graded anharmonic oscillator Fermi-Pasta-Ulam lattice both with nearest-neighbor (NN) and next-nearestneighbor (NNN) interactions. The system presents a maximum rectification efficiency for a very precise value of the parameter that controls the coupling strength of the NNN interactions, which also optimizes the rectification figure when its dependence on mass asymmetry and temperature differences is considered. The origin of the enhanced rectification is the asymmetric local heat flow response as the heat reservoirs are swapped when a finely tuned NNN contribution is taken into account. A simple theoretical analysis gives an estimate of the optimal NNN coupling in excellent agreement with our simulation results.


I. INTRODUCTION
While electronics has been able to control the charge flow for decades, the manipulation of heat current remains elusive more than 80 years after the first experimental observation of asymmetric heat flow in solids [1].The envisaged device capable of producing heat rectification can be generally described as a system wherein the heat current changes its magnitude as the position of the thermal reservoirs at both ends of the sample are swapped.Such a device is known as a thermal diode and is the fundamental tool for the manipulation and control of the heat flow.Despite the large amount of theoretical studies performed in the last decade [2,3], the progress in the field remains largely stalled in terms of management of thermal currents due to the lack of efficient and feasible diodes.Most of the theoretical proposals for a thermal diode are based on the sequential coupling of two or three segments, these being oscillator lattices, with different anharmonic interactions that provide the necessary asymmetry for thermal rectification to exist [4][5][6].One of the most severe problems of these devices is that they have a small rectification efficiency that rapidly decays to zero as the system size increases [7,8].Therefore, searching for alternative rectifying systems to overcome the deterioration of the rectification effect for large systems is of great current interest.
Besides the coupling of a few asymmetric segments [4][5][6] another way to produce the required asymmetry for thermal rectification to occur is to consider a mass gradient along the system length.This idea was first employed to experimentally build a rectifying device with a carbon and boron-nitride nanotube inhomogeneously mass-loaded with heavy molecules [9].This concept was later implemented in an anharmonic oscillator lattice with nearest-neighbor (NN) interactions and with a linear mass gradient along its length [10].This configuration shows good rectifying properties and is proven to be robust against the variation of various structural parameters [11] and to the change in the form of the mass distribution [12].The idea of employing a graded property has been further applied to a closed billiard model with a graded magnetic field along its length to obtain thermal rectification [13] and also to a chain of elastically colliding, asymmetrically shaped massgraded particles, where rectification is theoretically predicted when the system has a local heat flux proportional to the temperature gradient and a local conductivity dependent on temperature [14].Much theoretical insight has been obtained with the study of a graded harmonic lattice with a quartic on-site anharmonic potential, self consistent heat baths and weak interparticle interactions [15].For this model it has been recently shown that the conditions for the existence of thermal rectification are the existence of a temperature gradient in the bulk, local conductivity dependent on temperature, and a graded structure [16].Remarkably, if one further includes harmonic long-range interactions in the model, an increase in rectification efficiency is obtained [17].However, in the latter analytical study only infinitesimal temperature differences between reservoirs were considered, which in turn led to infinitesimal rectification.
The relevance of long-range interactions in order to achieve efficient thermal rectification that does not vanish with the system size has been further investigated by nonequilibrium molecular dynamics simulations in Ref. [18].In that paper a system consisting of a one-dimensional (1D) mass-graded lattice of oscillators with harmonic long-range interactions and an anharmonic on-site potential was investigated, wherewith significant rectification factors, which do not decrease with the system size, were achieved with a moderate mass asymmetry.
From the above mentioned studies one may be tempted to conclude that long-range interactions play an essential role in improving the efficiency of thermal diodes based on asymmetric lattices.The question that arises is whether one could achieve efficient rectification at all with short-range potentials and interactions in this type of system.Therefore, it would be very interesting to explore alternative short-range models that take advantage of the amplifying effect on thermal rectification of interactions beyond the NN range leading to highly effective thermal diodes even as the system size is increased.
In this paper we address the problem of improving the thermal rectification of a mass-graded anharmonic lattice [10,11] by the addition of next-nearest-neighbor (NNN) interactions.Thus this system allows to study in a systematic way the effect on thermal rectification of forces with significant magnitude beyond the NN range.An additional motivation stems from the fact that systems with such a type of interactions present unusual dynamical features.For example, a crucial modification to the multibreather properties [19] and a nontrivial intrinsic phase structure [20] of solutions have recently been observed in the Klein-Gordon lattices and in the discrete nonlinear Schrödinger lattices, respectively, when the NNN interactions are considered.Thus it can be expected that such unusual microscopic dynamics should certainly have an influence on the macroscopic properties of this system.An example of such a connection is afforded by a 1D oscillator lattice with randomly distributed NNN interactions, being thus a model of a glass, has been used to determine that the observed energy localization increases the relaxation time of the system [21].Furthermore, it has been determined that NNN interactions crucially contribute to the nonuniversality of the divergence exponent of the thermal conductivity with the system size [22,23].Therefore, it is reasonable that there might be an influence of the NNN interactions on the rectification properties of the herein considered anharmonic oscillator lattice.
This paper is organized as follows: in Sec.II the model system and methodology are presented.Our results on the dependence of rectification on the strength of the NNN interactions and other involved parameters are reported in Sec.III.The discussion of the results, as well as our conclusions, are presented in Sec.IV.

II. SYSTEM DESCRIPTION
The reference Hamiltonian for the 1D model we are considering can be written as where N is the system size and V (x) is the potential, which corresponds to the Fermi-Pasta-Ulam (FPU) β model [24] characterized by V (x) = x 2 /2 + x 4 /4, {m i ,q i ,p i } N i=1 are the dimensionless mass, displacement, and momentum of the ith oscillator; fixed boundary conditions are assumed (q 0 = q N +1 = 0).The tunable parameter γ specifies the relative strength of NNN coupling compared to the NN one.We consider a mass-graded lattice wherein the mass of the ith oscillator is given by m i = M max − (i − 1)(M max − M min )/(N − 1), i.e., a mass gradient decreasing from left to right.M max and M min are the mass of the oscillators at left and right ends, respectively.In the following we will always take M min = 1.The equations of motion for each lattice oscillator can be written as qi = p i /m i and ṗi where F (x) = −∂ x V (x) is the force and ξ (i)  L/R is a Gaussian white noise with zero mean and correlation R (taken as = 0.1 in all computations hereafter reported) the coupling strength between the system and the left and right thermal reservoirs operating at temperatures T L and T R , respectively; furthermore, The equations of motion (2) were integrated with a stochastic velocity-Verlet integrator [25] with a time step of t = 0.01 for 8 × 10 7 time units after a transient of 6 × 10 7 time units.
Once the nonequilibrium stationary state is attained, the local heat flux is computed as In the stationary state the local heat flux J i becomes constant along the chain apart from thermal fluctuations, i.e., J i ∼ J , where J is the total heat flux.However, to improve the accuracy its value is calculated as the algebraic average of J i over the system bulk, i.e., J = (1/N) is the number of oscillators in the bulk, i.e., excluding the reservoirs.We use J + to denote the heat flux obtained when the high-temperature reservoir is attached to the left end of the lattice and J − when that same reservoir is attached to the right end.The thermal diode efficiency r is defined as the ratio between the heat current transmitted in the forward temperature bias configuration (J + ) and that generated when the T -bias is reversed (J − ) leading to r = |J + /J − |.A highly efficient diode corresponds to r 1.In the following the behavior of r is studied as a function of the temperature difference T = T L − T R and for a constant value of the average temperature T 0 ≡ (T L + T R )/2 of the system.

A. Maximum rectification efficiency
In Fig. 1 we present the results of the dependence of the rectifying efficiency as a function of the relative strength of the NNN potential γ for various system sizes N and a fixed value of M max = 10.It is clear that, for small γ values, r tends to the value in the case when there are no NNN interactions.For large γ values the second term of the potential, i.e., the NNN contribution, dominates over the NN term and r decreases to a value close to that corresponding to the one term potential.This suggests that both contributions are needed in order to achieve the desired rectification effect.This behavior leads to the existence of a maximum for the rectification efficiency at the very precise value of γ c = 0.45 for the largest system size considered (for lower N values small system size dependent corrections appear).We have corroborated that, for greater M max values, the maximum rectification is again obtained for γ = γ c .Superficially our results bear some resemblance to those obtained with a lattice with harmonic long-range interactions [18].However, in that case the maximal rectification could be easily associated with the value of the control parameter that maximizes the strength of the long-range potential; in contrast, no such association can be made in the present case.At this point we recall that, for the very same model herein considered, the divergence exponent of the thermal conductivity with the system size changes continuously as γ increases [22], just as in the herein considered dependence of r on γ observed in Fig. 1 and which seems to indicate that the maximum rectification is constrained to be in in some narrow range of γ values.Remarkably, the existence of an optimum γ can be justified analytically under certain assumptions as will be shown in Sec.III B, seemingly being generic in this model and independent of the temperature difference and mass distribution.
In Fig. 2 we present the temperature profiles, i.e., the local temperature computed as T i = p 2 i /m i , for both configurations of the thermal reservoirs employed in the case of a lattice which has only NN interactions and other with extra NNN interactions with γ = 0.45.In both instances the asymmetry in the shape of the profiles when the thermal reservoirs are swapped becomes evident, which indicates the presence of thermal rectification.Notice that the difference between the profile curves for γ = 0 (no NNN interactions) and γ = 0.45 is far more pronounced for T L < T R , which is related to the existence of NNN interactions since the asymmetric mass distribution is the same in both instances.Thus these could be held responsible for the rectification increase observed in Fig. 1 of the lattice with NNN interactions with respect to the one with only NN ones.In Fig. 3 we plot separately the contribution of the NN term J (1)   i to the local heat flux, that of the NNN one J (2)  i , and the local heat flux J i given by Eq. ( 3) for a lattice with bulk size N = 200.It can be readily observed that, for the T L > T R case, i.e., when the applied temperature gradient is in the same direction as mass gradient, the contribution of the NN term decreases in the direction of the total heat flow along the lattice, whereas the contribution arising from the NNN interactions increases, both in a linear way.If we take the spatial average of both local heat fluxes we obtain J (1) + = 9.7 × 10 −4 and J (2)  + = 7.3 × 10 −4 , with a total heat flux of J + = 1.7 × 10 −3 .However, when the reservoirs are swapped, we observe that the aforementioned linear behavior J i FIG. 4. Same as Fig. 3, but now for a lattice of N = 2000 oscillators.(a) T L > T R and (b) T L < T R .
is not maintained along the whole length of the lattice; rather we observe a nonlinear behavior towards the heavy loaded end connected to the cold reservoir.In this case, if we again take the spatial average, we obtain J (1)  − = −2.3× 10 −4 and J (2) − = −3.0× 10 −4 , with J − = −5.3× 10 −4 and J (1)  − ∼ J (2)  − .In both instances it is clear that the NNN contribution to the local heat flux is stronger in the end connected to the cold reservoir, but for the T L < T R configuration this enhanced contribution is responsible of making the complete local heat flux J − less negative, thus diminishing its absolute magnitude and increasing the rectification efficiency of the system.We have also studied the cases with γ = 0.01 and γ = 1 (not shown) for the configuration T L > T R .In the first case the aforementioned behavior of the local heat flux associated with the NNN interactions is present but with a magnitude much less than that of the one corresponding to the NN interactions, thus making its contribution to the total heat flux J irrelevant.In this latter case both J (1)  i and J (2)  i are on average constant along the lattice and of similar magnitude; so, the rectification value has a value similar to the case wherein there are no NNN interactions present.Therefore we can conclude that γ = 0.45 represents an intermediate value wherein the mechanism detected in Fig. 3(b) is maximum, being of negligible magnitude for γ = 0.01 and completely absent for γ = 1.
In Fig. 4 we present the same results as in Fig. 3 but now for a system of N = 2000 oscillators.For the T L > T R case both contributions in Eq. ( 3) to the local heat flux at each site are almost constant along the bulk of the system, with small but nevertheless noticeable deviations at the system borders.The corresponding spatial averages are J (1) + = 1.6 × 10 −3 and J (2)  + = 8.9 × 10 −4 , with a total heat flux of J + = 2.5 × 10 −3 .For the configuration T L < T R , on the other hand, we immediately notice that the aforementioned boundary behavior is much stronger and affects the local heat flux in the bulk of the system.This boundary effect is so strong that right in the edges the NN contribution to the local heat flux in the hot end and the NNN one in the cold end become positive.In this configuration J (1)  − = −4.2× 10 −4 , J (2) − = −4.1 × 10 −4 , and J − = −8.3× 10 −4 .Therefore, just as in the N = 200 case, we have again that J (1)  − ∼ J (2) − .

B. Analytical approximation to the maximum rectification efficiency
In the stationary regime and for either of the reservoir configurations considered we have J ± = J (1)  ± + 2γ J (2)  ± .This expression, together with the definition of the rectification efficiency so far employed, will allow us to perform an analytical approximation to the maximum rectification efficiency γ c computed numerically in order to gain further insights of the conditions wherein we can expect that result to hold up.
We are looking for γ = γ c such that (∂r/∂γ )| γ =γ c = 0, so we calculate which becomes zero for γ = γ c satisfying where ± , being sgn(h) the sign function and neglecting any unaccounted dependence of J ± on γ .So, within this approximation, the maximum condition (5) becomes sgn(J + )J (2)  + (J (1) From here we can solve for γ c and obtain This solution exists if and only if the condition is fulfilled, otherwise the denominator is zero and there is no maximum of r(γ ).Therefore, whenever ( 8) is fulfilled we find the unique solution Now in order to proceed any further we make the approximation |J (1)  ± | ≈ |J (2)  ± |, which is rather good for the numerical data obtained from Fig. 3 (N = 200) and less accurate for those corresponding to Fig. 4 (N = 2000).Within this approximation we arrive at which is a low order approximation to γ c but seems to be rather good when compared with the value obtained from simulation in Fig. 1.This analysis also shows that γ c may be independent of the details of the specific mass distribution employed or temperature bias, in agreement with the numerical results as well.The question of how good the employed approximation is to obtain γ c in the large system size limit remains open, but according to our numerical results (see next section and Fig. 7) the rectification efficiency saturates for large systems.This suggests that, even if deviations from the theoretical value of γ c = 1/2 might exist as we go to very large system sizes, the optimal rectification efficiency reaches its asymptotic value already for system sizes of order N = 500.

C. Rectification dependence on model parameters and system size
In the following we analyze the dependence of the rectification efficiency on other parameters that determine the behavior of the system.In Fig. 5 we plot r as a function of the temperature difference T for the case of NN and NNN interactions, fixed system size, and M max = 10.It is evident that there is an increase of in the latter case compared to the former.In Fig. 6 we plot r as a function of M max and thus corroborate that by increasing the asymmetry of the system we obtain an increase of rectification.The increment in r follows the same pattern as the NN case also depicted, but with a larger magnitude for each M max value considered; it can also be inferred that the rectification tends to saturate as the asymmetry increases, and thus, at some point, further increments of M max will not result in a significant increment in r.Now, in these last two figures it is clear that the γ c = 0.45 value optimizes the rectification efficiency of the system for each of the considered conditions.
As for the dependence of r on system size, in Fig. 7 we report our results.For the NN case rectification decays as the lattice length increases, whereas when NNN interactions are included r steadily increases up to a size around N = 500.For even larger lattices, r is almost constant with a extremely slow decay with N .At this point it is important to recall the results depicted in Fig. 4, wherein the possible origin of the increase of rectification with NNN interactions was attributed to a boundary effect affecting the bulk of the system.A possible corroboration to this conjecture is that, for N = 2000, r = 3.04, whereas for N = 200 we have r = 3.2.Since rectification has a small decrease as the system size increases, then it is reasonable to suspect that the aforementioned conjecture could be correct, since boundary conditions tend to be less important in the large system size limit.Additional simulations would be needed to further explore this behavior, although the large simulation times that are needed to attain the stationary state, which existence and uniqueness was proven in Ref. [26], would make very difficult to study the rectification in larger system sizes than those depicted in Fig. 7.

IV. CONCLUDING REMARKS
In this work we have performed the study of the rectification properties of an anharmonic mass-graded lattice with NNN interactions.For the value range of the considered parameters we determined that rectification increases in all instances compared to the NN case.The thermal rectification efficiency is maximized for a certain value of the coupling constant of the NNN interaction γ c .This optimal coupling can be estimated theoretically on the basis of the approximation that the NN and NNN contributions to the local heat flux are similar in magnitude, which seems to agree with simulations for moderate system sizes.Within this approximation we obtain the optimal coupling to be γ c = 1/2 in good agreement with simulations for the system sizes considered.Increasing the NNN coupling above this value systematically worsens the rectification efficiency of the system.According to the theoretical approximation, and corroborated by our simulations, γ c is independent of the specific masses or temperature bias.The rectification efficiency r increases quite rapidly with the system size until it reaches some plateau for large sizes.The obtained numerical data show that once this asymptotic rectification is reached it decays extremely slowly with N .
The anharmonic lattice herein considered allows a verification of the behavior of the rectification efficiency for large systems with NNN interactions.Indeed, our model would allow us to explore in a systematic way the influence of interactions of gradually increasing range (third, fourth neighbors and so on) on the thermal rectification efficiency of this and/or similar lattice models.This could be interesting if we consider that, although the employed interactions theoretical models are long-range, it could also be interpreted, perhaps more realistically, that the experimental systems are actually small with respect to the interaction range.Indeed, much theoretical work on phononic thermal diodes based on anharmonic lattices with long-range interactions was inspired by experimental results on atomic clusters of a few atoms [27][28][29].Notwithstanding those experiments, it is important to bare in mind that other proposed thermal diodes, like those based in permalloy and other magnetic materials [30], are indeed generic long-range interacting systems where heat conduction has an electronic or spintronic rather than phononic origin.

FIG. 1 .
FIG. 1. Thermal rectification r vs the relative strength of the NNN potential γ .For M max = 10 we plot r for N = 128 (circles), 256 (squares), and 512 (triangles) with T 0 = 0.1, T = 0.16, and n L = n R = 1.Void symbols depict the corresponding values for the γ = 0 case.Vertical dotted line indicates the value γ c = 0.45.Continuous lines are a guide to the eye.

FIG. 2 .
FIG. 2. Temperature profiles of a mass-graded lattice with N = 200; same M max = 10, T 0 , T , and n L,R values as in Fig. 1. Results for T L > T R (circles) and T L < T R (triangles) are shown.Solid lines correspond to a lattice with NNN interactions and γ = 0.45 and dashed lines to the case with NN interactions only (γ = 0).

( 2 )FIG. 3 .
FIG. 3. Local heat flux for a mass-graded lattice with M max = 10, N = 200, and γ = 0.45; same T 0 , T , and n L,R values as in Fig. 1.Dashed lines with triangles indicate the NN contribution, dot-dashed lines with circles the NNN one, and solid lines the result according to Eq. (3) for (a) T L > T R and (b) T L < T R .

1 FIG. 5 .
FIG.5.Dependence of thermal rectification r on temperature difference T for three values of γ : 0 (circles), 0.45 (squares), and 1 (triangles), all with N = 512.Same values of M max , T 0 , and n L,R as in Fig.1.Continuous lines are a guide to the eye.

1 FIG. 6 . 1 FIG. 7 .
FIG.6.Thermal rectification r as a function of the largest mass M max for N = 512 and various γ values.Same T 0 , T , and n L,R values as in Fig.1.Continuous lines are a guide to the eye.