Thermal rectification in silicon by a graded distribution of defects

Riccardo Dettori, Claudio Melis, Riccardo Rurali, and Luciano Colombo Dipartimento di Fisica, Universit a di Cagliari, Cittadella Universitaria, I-09042 Monserrato (Ca), Italy Istituto Officina dei Materiali, CNR-IOM SLACS Cagliari, Cittadella Universitaria, Monserrato (Ca) 09042-I, Italy Institut de Ciència de Materials de Barcelona (ICMAB-CSIC), Campus UAB, 08193 Bellaterra, Barcelona, Spain Catalan Institute of Nanoscience and Nanotechnology (ICN2), CSIC, and The Barcelona Institute of Nanoscience and Nanotechnology, Campus UAB, 08193 Bellaterra, Barcelona, Spain


I. INTRODUCTION
Thermal rectification [1][2][3] occurs whenever the heat flux is affected by the actual direction of the thermal gradient applied to the system.The amount of rectification of the heat current is usually quantified by the factor R ¼ j Jfwd j=j Jrev j À 1, where Jfwd and Jrev are the heat fluxes corresponding to the forward and reverse thermal bias conditions, respectively.The situation is conceptualized in Fig. 1 for the prototypical configuration corresponding to an interface between two materials A and B. In Fig. 1, we assumed that Jfwd > Jrev , but we remark that the definition of "forward" or "reverse" bias condition is only a matter of convention.6][7] Here, the generation, control, and manipulation of lattice heat (or, equivalently, phonon flux) are the main tools to engineer devices with functionality similar to their electronic counterparts (like electrical diodes or transistors).
In most cases of practical interest, 5 thermal rectification is observed when two materials with unlike thermal conduction properties are interfaced, as indeed shown in Fig. 1.In this configuration, the role played by the sharp interface is crucial: here a temperature drop DT fwd;rev occurs in whatever bias condition, giving rise to an interface thermal resistance (ITR) quantified by the ratio jDT fwd;rev j=j Jfwd;rev j as extensively discussed in Refs.2, 3, and 8.While ITR is largely referred to in semiconductor systems, also interfacial hybrid systems like metal/superconductor 9 or crystal/polymer 10 junctions have been shown to provide thermal rectification.However, the actual rectifying properties of an interface are affected by many features (structural details, chemical contamination, and interdiffusion to name just a few) which require nontrivial nanofabrication techniques for their government, possibly resulting in a rather difficult technological task to be accomplished.
This scenario suggests that it would be interesting to observe thermal rectification without any localized (i.e., abrupt) temperature drop.This basically requires that such a sharp interface is not present in the rectifying device.In this regard, low-dimensional systems (either model or realistic) prompt several possible rectifying configurations, such as non-uniform mass loading, suited distributions of defects, and tailored shaping, [11][12][13][14][15][16][17][18][19][20]33 where the above combination is actually realized. In his work, we further exploit this concept in bulk-like systems.In particular, by nonequilibrium molecular dynamics (NEMD) simulations, we address the thermal rectification factor R in Si bulk structures containing a gradual distribution of compositional or structural defects.The reason to select Si for the present investigation is twofold, namely, (i) it has an impact on the emerging phononic nanotechnology quoted above and (ii) it has very wellknown thermal properties.

II. SIMULATION SETUP
The NEMD simulations have been performed on the Si tetragonal cells with section S ¼ n Â n a 2 0 (n ¼ 10 or 13, as indicated below) normal to the (001) direction of heat transport (hereafter named z) and a total length L z varying in the range of 100-350 a 0 (as indicated below).At the left and right extrema of such a simulation cell, two further slabs of thickness 10 a 0 were added and coupled to thermostats (see below).Periodic boundary conditions have been applied along the two directions normal to z.We set a 0 ¼ 0.543 nm as predicted by the Tersoff potential providing the force field for the present investigation. 21The simulation cells have been at first filled by a diamond lattice of silicon atoms.Then, a nonuniform dispersion of defects was obtained: (i) by randomly replacing Si atoms with Ge ones, up to a 20% of Ge content, which is the minimum doping corresponding to the maximum reduction of lattice thermal conductivity in a SiGe alloy 22,23 or (ii) by removing clusters of Si atoms so as to create voids with a random position, size, and shape up to a maximum 31% porosity.In any case, the two thermostatted slabs were not affected by Si !Ge replacements or void generation.This procedure generated a defect distribution characterized by different concentration profiles along the z direction, as shown in Fig. 2 (where the thermostatted slabs are not shown for the sake of clarity).A structural relaxation followed, respectively, (i) through a careful energy minimization by conjugate gradients or (ii) through a high-temperature simulated annealing at 900 K.The latter procedure (implemented over 5 Â 10 5 time steps of duration 0.5 Â 10 À15 s) was indeed required in order to allow for the full reconstruction of dangling bonds created by atom removals.In turn, the reconstruction generates a shell of amorphous matter at the void surface, providing an important source of phonon scattering in nanoporous Si. 24 In the case of a SiGe graded alloy, we have also taken care of any possible structural relaxation along the z direction due to the Si-Ge lattice mismatch by further performing constantpressure, constant-temperature MD simulations as long as 5 Â 10 5 time-steps, each lasting 1.5 Â 10 À15 s.Since we observed just a very minor variation of the cell length and, in any case, no detectable effect in the calculation of the heat flux (see below), we present here the results obtained with constant-volume cells.To this aim, fixed boundary conditions were imposed along z by adding one more plane on both sides where atomic positions have been clamped anytime during the simulations.Finally, for the sake of comparison in both cases, it has generated a sharp interface between pure Si and, respectively, a homogeneous SiGe alloy with a 20% of Ge content or a nanoporous Si sample with 31% porosity.
The desired steady state condition where to investigate possible rectification effects was generated in each system by coupling its left and right 10a 0 -thick terminal slabs (see above) to a heat reservoir, respectively, set at T h ¼ 700 K and T c ¼ 500 K by using Nos e-Hoover thermostats.These temperature values guarantee that most part of the system stays above the Si Debye temperature and, therefore, quantum effects play a very minor role.While the initial temperature of each sample was set at the average ðT h þ T c Þ=2 value, the MD simulation was aged until a steady state regime was reached for the selected thermal bias condition (assigned by the relative position of the hot and cold thermostats).Given the very small resulting thermal conductivity, this required as many 3 Â 10 6 or 7 Â 10 6 time steps of MD for the Ge-doped and nanoporous structures, respectively.In order to assess the steady state, we calculated both the works W in and W out , respectively, spent by the hot and cold thermostats, and the corresponding heat fluxes j Jin;out j ¼ ð1=SÞð@W in;out =@tÞ.A steady state condition was proclaimed only when j Jin j ¼ j Jout j within the accepted numerical error (which we have set at 1%).This further required additional 1 Â 10 6 or 3 Â 10 6 time steps, according to the system.Once reached such a condition, the steady state heat flux for the assigned bias (i.e., j Jfwd j or j Jrev j) was calculated as the average ðj Jin j þ j Jout jÞ=2.By inverting the two thermostats and repeating the calculation, the other heat flux was similarly calculated and the rectification R was eventually obtained.In all samples investigated here, we named "forward" the bias condition where the pure silicon part of the system (i.e., the left end in Fig. 2, all panels) was set at T h.All simulations have been executed by using the LAMMPS code. 25

III. RESULTS AND DISCUSSION
Table I summarizes the results obtained for the configurations shown in Fig. 2, reporting rectifications in the range of 2.0%-3.5% and 1.4%-3.2%for Ge and pore distributions; on average, the error in estimating rectification is about 0.45% and 0.32%, respectively.These data provide evidence that not only a rectification is indeed found but also it is ruled over by changing the distribution of Ge atoms or pores along the z direction.Overall, the predicted R is comparable to what was observed in other low-dimensional Si-based systems, [26][27][28][29][30][31] proving that rectification is indeed possible in a bulk-like system lacking sharp interfaces.Interestingly enough, rectifications as small as 3%-4% have been indeed measured in the Si-based systems; 14 an even smaller rectification of about 1%-2% has been experimentally reported, although in a rather different system as reduced graphene oxide. 18So, the rectification values predicted in this work should be within the experimental capability of measure.
FIG. 1. Left: schematic representation of the "forward" and "reverse" thermal bias conditions.Rectification occurs whenever Jfwd 6 ¼ J rev .Right: zoomed interface region where the temperature drop DT occurs (full black line), generating localized thermal resistance.
Varying the profile of the defect distribution is an effective way to control the resulting R, and the present simulations suggest that the exponential profile turns out to be the most efficient in generating different values for j Jfwd j and j Jrev j.This is not, however, the only way to tune rectification features.In fact, we calculated the dependence of the rectification on the value of the imposed temperature offset T h À T c ¼ DT, as shown in Table II for the same L z ¼ 100 a 0 sample containing a graded distribution of Ge defects with similar exponential profile.The results are along the expectations: by decreasing the temperature offset between the hot and cold thermostats, the rectification is reduced from 3.5% to 2.7%.On average, the error in estimating the rectification for these systems is 0.37%.Interestingly enough, however, such a reduction is weak: while the temperature offset was reduced by a factor 4Â, the calculated rectification is only reduced by a factor 1.3Â.We believe this is an interesting result, making clear that the predicted rectification feature is robust.Interestingly enough, the rectification is also affected by varying the absolute temperature of the two thermostats, but still preserving their offset: as a matter of fact, by repeating the calculation for an exponential profile of Ge substitutional defects with T h ¼ 900 K and T c ¼ 700 K we obtained R ¼ 4.3%.We attribute such an increased rectification to the different average interface temperatures. 32nother intriguing feature of the rectification phenomena reported here is that they are not paralleled by the onset of any interface temperature drop, possibly causing ITR effects.Fig. 3 reports the calculated temperature profiles in the steady state conditions for all graded systems shown in Fig. 2, both in the forward (open blue symbols) and reverse (red full symbols) thermal bias conditions.The profiles have been obtained by calculating the local temperature of slabs as thin as 2.7 nm aligned along the z direction, over which the atomic velocities have been averaged for 1 Â 10 6 and 3 Â 10 6 time steps, according to the system.
We fathom the present results in terms of a nonseparable dependence of the thermal conductivity upon the z-coordinate and temperature T which, in turn, defines a non-linear heat transfer regime.Such a non-linear regime is the key feature for rectification, which cannot be simply ascribed to an asymmetric scattering of the heat carriers by defects when inverting the thermal bias.As a matter of fact, the structural inhomogeneity generated by the graded distribution of compositional or structural defects makes the thermal conductivity a function of z.On the other hand, it turns out that the same quantity is explicitly a function of T as well, since the system is out of equilibrium (although in steady state).In other words, we take for granted that each single portion of the system is transmitting heat while experiencing a different local temperature with respect to the temperature of its neighbouring regions.This is tantamount to say, that in all investigated samples the function j ¼ jðz; TÞ FIG. 2. Graded distribution of substitutional Ge defects (left) and pores (right) in a Si lattice (light blue).From top to bottom, it is shown a linear, quadratic, exponential, and step-like distribution of defects.Their average concentration is shown by a yellow dotted-line.Pictures show a 4a 0 -thick longitudinal section of each sample.TABLE I. Rectification calculated for the graded distributions of Ge atoms or pores shown in Fig. 2. On average, the error in estimating rectification is about 0.45% and 0.32%, respectively.For all samples, the temperature offset between the hot and cold thermostats is set at DT ¼ 200 K and it is centred at an average temperature of 600 K. Simulation cells have a 13 Â 13 a 2 0 section and a length L z ¼ 100 a 0 .2) as a function of the temperature offset DT between the hot and cold thermostats (in all cases, the average temperature is 600 K).On average, the error in estimating rectification is 0.37%.The simulation cells have a section S ¼ 13 Â 13 a 2 0 and a length L z ¼ 100 a 0 .

Linear
is non-separable.Let us now assume that, contrary to the above conclusion, the thermal conductivity is separable, i.e., we can write jðz; TÞ ¼ f ðzÞgðTÞ, where f(z) and g(T) are known functions.In the steady state condition (whatever thermal bias) investigated here, the heat equation for onedimensional transport along the z direction, can be easily integrated by variable separation since the onedimensional heat flux J z is a constant and, therefore, we get where T l and T r are the temperatures of left and right terminal ends of the system, respectively, located at positions z l and z r .Equivalently, Eq. ( 2) can be cast in the form By inverting the thermal bias condition, we just overturn the upper and lower limits in the temperature integral: this will only affect the sign of the heat flux, leaving unaffected its absolute value.This implies a null rectification, i.e., R ¼ 0 since j Jfwd j ¼ j Jrev j.Therefore, the assumption that kðz; TÞ ¼ f ðzÞgðTÞ is separable has in fact defined a sufficient condition for no rectification. 33his is the key concept that allows to understand our results, providing a rationale for them.Sure enough, we can logically invert the above statement and say, that a nonseparable j ¼ jðz; TÞ form of the thermal conductivity does represent the necessary condition for rectification.This is precisely what is exploited by the combination of (i) a graded distribution of defects and (ii) a thermal bias condition.Therefore, the bulk structures here investigated must rectify a thermal current: as a matter of fact, their thermal conductivity is a complicated and non-linear convolution given by a z-dependence of the temperature which, in turn, is a function of the local stoichiometry or porosity.
The present picture on rectification is robust since it does not qualitatively depend neither on the nature of the defects (compositional or structural) distributed in the bulk structure nor on their actual distribution profile.It is also found for no matter what thermostatting condition is set: by simulating a graded exponential distribution of Ge defects under two different conditions (namely, T h ¼ 700 K with T c ¼ 500 K and T h ¼ 900 K with T c ¼ 700 K) in both forward and reverse bias, a smooth continuous temperature profile was found in all the samples, similar to that shown in Fig. 3.
When the same analysis is applied to the systems characterized by a step-like distribution of defects (Fig. 2, bottom FIG. 3. Temperature profiles calculated for the graded distributions of Ge atoms (left) and pores (right) shown in the same order of Fig. 2 (the step-like profile is omitted here for sake of clarity).For all systems, it has been set T h ¼ 700 K and T c ¼ 500 K.The forward and reverse thermal bias conditions correspond to the empty (blue) and full (red) symbols, respectively.Errors are indicated by vertical bars.FIG. 4. Zoomed temperature profiles nearby the interface (positioned at z ¼ 50 a 0 ) calculated for a step-like distribution of Ge atoms (left) and pores (right) in the T h ¼ 700 K and T c ¼ 500 K thermostatting condition.The resulting temperature drop DT is shown for both the forward (empty blue symbols) and reverse (full red symbols) thermal bias conditions.Errors are indicated by vertical bars.panels), we have indeed found a small (but definitely non-vanishing) and abrupt temperature offset at the interface, as reported in Fig. 4. In the case of Ge doping, we evaluated such an interface temperature drop as large as DT fwd ¼ 18.1 6 2.3 K and DT rev ¼ 21.6 6 2.6 K for the forward and reverse bias situation, respectively.Similarly, for a step-like distribution of pores we calculated DT fwd ¼ 11.9 6 4.2 K and DT rev ¼ 17.8 6 3.0 K, resulting in a more relevant difference between the two thermal bias conditions.
We now investigate another important issue, namely, how the predicted rectification is affected by the gradient of the defect distribution.To this aim, we have once again selected an exponential profile of Ge substitutional defects which, as in the previous cases, was varied from the minimum 0% to the maximum 20% content over a sample length L z .However, in this case we considered four different samples with increasing thickness along the direction of heat transport corresponding to L z ¼ 100, 200, 300, and 350 a 0 , respectively.The resulting defect profiles are shown in Fig. 5.The two largest increased lengths correspond to quite a big simulation cell, and therefore, in order to keep the corresponding computational workload sustainable, we reduced the sample cross section to 10 Â 10 a 2 0 .This reduction makes unfair the direct comparison with the previously investigated sample with same L z but larger S, and therefore, we recalculated the rectification for the new section.Results are shown in Table III, indicating that the rectification is predicted to decrease from 4.3% (L z ¼ 100 a 0 ) to 3.3% (L z ¼ 350 a 0 ), i.e., a 3.5 Â increase of the sample length has reduced rectification by only a factor 1.3Â.On average, the error in estimating rectification is 0.75%.So, as expected, there is indeed a reduction in the rectification features, but even in this case the dependence on the sample length is weak.
It is hard to assess whether this reduction is due to the increased (because of the increased sample length) scattering of heat carriers or, rather, to the decreased gradient in the defect distribution.In order to further substantiate the argument that the interplay between the sample length and the gradient of the defect distribution is complex, we considered one more configuration with S ¼ 13 Â 13 a 2 0 ; L z ¼ 100 a 0 and an exponential profile of Ge substitutional defects which, however, was now varied from 0% to 60% of Ge-content.Interestingly enough in this case, we found a smaller rectification than reported in Table I, a value in fact very similar to the rectification calculated for the step-like profile.A systematic set of simulations exploring various combinations of length and gradient effects would be needed to fully clarify this issue.We leave them to a following investigation.
Finally, by increasing the maximum Ge content of the doped region up to 60%, or by enlarging the ðT h À T c Þ difference up to 400 K, or even by reducing the thermostats temperature to T h ¼ 400 K and T c ¼ 200 K (and neglecting possible quantum effects), we have found full confirmation of the picture outlined above, namely, rectification up to $5% is always observed in graded bulk structures without any ITR (once again, because of the missing sharp temperature drop anywhere in the system).

IV. CONCLUSIONS
In conclusion, we have shown that the thermal rectification can be obtained in a bulk-like silicon system missing of any interface by exploring two different configurations where a graded distribution of Ge substitutional defects or nanovoids, respectively, is arranged along the direction of an applied thermal gradient.This result is consistent with the general argument that the rectification is necessarily generated in any system where the thermal conductivity is a nonseparable function of both temperature and position, 33 indeed a situation found in any structure is investigated here.For both graded distributions, we proved that the rectification is obtained still preserving a very smooth temperature profile throughout the system, i.e., without setting up any localized temperature drop.Furthermore, it is found that the resulting value of R depends on the structural features of the graded distribution, as well as on the applied thermal bias and the actual gradient of the defect distribution.This, in turn, suggests a possible way for tuning thermal rectification by defect engineering.Finally, we remark that for any system investigated here the predicted R is comparable with the rectification observed in low-dimensional Si-based systems and lies within the experimental resolution.

ACKNOWLEDGMENTS
We acknowledge the financial support by the Spanish MINECO under Grant Nos.FEDER-FIS2012-37549-C05-02, FEDER-MAT2013-40581-P, TEC2012-31330, and TEC2015-67462-C2-1-R, the Generalitat de Catalunya under Grant Nos.2014 SGR 301 and 2014 SGR 384, and the FIG. 5. Local content of Ge substitutional defects along the direction z of heat transport for four samples with different length L z .Symbols (connected by thick lines) represent the actual Ge content; thin lines represent a guide to the eye, corresponding to an ideal exponential profile.TABLE III.Rectification calculated for the graded distributions of Ge atoms with an exponential profile (see Fig. 2) as a function of the length L z of the simulation cell.On average, the error in estimating rectification is 0.75%.In this case, the cross section was reduced to 10 Â 10 a 2 0 for computational convenience.The temperature offset between the hot and cold thermostats is set at DT ¼ 200 K and centred at an average temperature of 600 K. L z ¼ 100 a 0 L z ¼ 200 a 0 L z ¼ 300 a 0 L z ¼ 350 a 0 4.3% 4.1% 3.6% 3.3%