Jesús
Carrete‡
a,
Miquel
López-Suárez‡
b,
Martí
Raya-Moreno
c,
Anton S.
Bochkarev
d,
Miquel
Royo
b,
Georg K. H.
Madsen
a,
Xavier
Cartoixà
c,
Natalio
Mingo
d and
Riccardo
Rurali
*b
aInstitute of Materials Chemistry, TU Wien, A-1060 Vienna, Austria
bInstitut de Ciència de Materials de Barcelona (ICMAB–CSIC), Campus de Bellaterra, 08193 Bellaterra, Barcelona, Spain. E-mail: rrurali@icmab.es
cDepartament d'Enginyeria Electrònica, Universitat Autònoma de Barcelona, 08193 Bellaterra, Barcelona, Spain
dCEA, LITEN, 17 Rue des Martyrs, 38054 Grenoble, France
First published on 12th August 2019
We combine state-of-the-art Green's-function methods and nonequilibrium molecular dynamics calculations to study phonon transport across the unconventional interfaces that make up crystal-phase and twinning superlattices in nanowires. We focus on two of their most paradigmatic building blocks: cubic (diamond/zinc blende) and hexagonal (lonsdaleite/wurtzite) polytypes of the same group-IV or III–V material. Specifically, we consider InP, GaP and Si, and both the twin boundaries between rotated cubic segments and the crystal-phase boundaries between different phases. We reveal the atomic-scale mechanisms that give rise to phonon scattering in these interfaces, quantify their thermal boundary resistance and illustrate the failure of common phenomenological models in predicting those features. In particular, we show that twin boundaries have a small but finite interface thermal resistance that can only be understood in terms of a fully atomistic picture.
Phonons are scattered at the interface between two different materials, giving rise to the phenomenon of thermal boundary resistance (TBR). The TBR is due to the breakdown of translational symmetry by the interface, and the difference in the mass of the constituent atomic species plays a prominent role. Accordingly, one would expect that phonons could pass through a crystal-phase interface, where the atomic mass does not vary, without being scattered. On the other hand, it has been recently shown that the thermal conductivity can decrease by as much as 40% from a cubic to a hexagonal polytype of the same material.6–8 Therefore, there is a need for detailed insight from atomistic calculations of heat transport across homojunctions and of their TBR. Understanding heat transport in these systems is important both for fundamental reasons and with a view to applications.9,10 On the one hand, phonon scattering at crystal-phase interfaces and, particularly, twin boundaries defies many phenomenological models used to account for interface thermal resistance; on the other hand, crystal-phase engineering can be used to design materials with tailor-made properties3,11–15 and twinning superlattices are playing an increasingly important role in nanowire science.16–19
In this work we use both state-of-the-art nonequilibrium Green's functions (NEGF) and nonequilibrium molecular dynamics (NEMD), based on carefully parameterized interatomic potentials, to calculate the TBR of a few prototypical homojunctions, compare them with the values obtained in similar, conventional heterojunctions and elucidate the atomic-scale mechanisms that originate interface scattering. These computational techniques provide complementary information: while NEGF accounts for the quantum nature of phonons in the low to mid temperature regime, NEMD naturally accounts for anharmonicity at all orders, which becomes more important at high temperature.
To address this concern in binary semiconductors, in this article we use a many-body interatomic potential that builds on the previous work of Vashishta and coworkers,25,26 and whose parameterization for GaP27 and InP28 reproduces the properties of different crystal phases, including ZB, WZ, and rock salt. As we are mostly interested in the vibrational properties, we calculated the phonon dispersions of the ZB and WZ crystals and compared them with those obtained by first-principles DFT calculations. We found that the ZB and the WZ phase of both GaP and InP are indeed described to a similar level of accuracy. Acoustic branches, whose accuracy is important in the study of interface or dislocation scattering,29 are reasonably well reproduced (see ESI†). There is a certain overestimation of the group velocity of the transverse modes in the case of GaP and InP, but it is similar for both polytypes.
In contrast, common interatomic force fields for Si, such as the Tersoff30 or the Stillinger–Weber potentials,31 give a rather poor description of vibrational properties.29 Therefore, we used machine-learning techniques to build a neural-network (NN) potential for Si that enables phonon dispersions with ab initio quality for both the cubic diamond (3C) and the hexagonal diamond (lonsdaleite, 2H) polytypes. The energy of the system is expressed as a sum of contributions from the environments of the different atoms, . An atom's environment is defined by the collection of distances, {rij}, between atom i and its neighbors, up to a cutoff radius rc. This collection of distances is transformed into an array of descriptors,
, for the environment of each atom i, where (μl, σl) are adjustable parameters for each member l of the array, and
is a smooth cutoff function. The contribution Ei is given by a neural network that has Dil as inputs, and the forces
on each atom k are the anaytical derivatives of E with respect to atom k's coordinates. The network weights, as well as the sets of parameters (μl, σl), are adjusted to minimize a cost function,
, of the differences between the network- and DFT-calculated energies and forces of all the training systems, s. This minimization is performed by stochastic gradient descent and backpropagation using tensorflow's32 GPU implementation. We used 32 different values of μ and α, and neural network with 2 hidden layers, i.e. 3 layers of coefficients, with 200 and 100 units, c = 0.01, obtaining a RMSE of 0.05 eV Å−1 for forces and 1.3 meV per atom for energies.
Even in the presence of a periodic reconstruction, two discrete translation symmetries can be identified along directions parallel to an interface between two crystalline regions. Therefore, elastic scattering of phonons at the interface is subject to the conservation of a two-dimensional parallel momentum q∥. The problem of heat-carrier reflection and transmission by an interface can thus be decoupled into an infinite set of one-dimensional scattering problems, one for each possible value of q∥ in a two-dimensional parallel Brillouin zone, BZ∥. The full solution to each of those problems consists in obtaining the detailed transmission probabilities per unit time among phonon modes on the left- and right-hand sides of the interface, constrained by the conservation of energy. For instance, for a mode incident upon the interface from the left lead with parallel and perpendicular components of the momentum (q∥, q⊥) and angular frequency ω, one first has to search for all the phonon modes in the left and the right sides compatible with both q∥ and ω, and then calculate the amplitude of each allowed scattering process. This involves exploring an effective one-dimensional Brillouin zone BZ⊥(q∥), looking for values of q⊥ and phonon branch indices yielding the same frequency as the incident phonon. The main difference with respect to a typical one-dimensional problem involving scattering by a barrier is the lack of time-reversal symmetry, which only applies to the line crossing the Γ point.
The detailed procedure described above is computationally expensive and has a tendency to numerical instability because of the need to solve the inverse problem of finding wave vectors as functions of frequency. Therefore, an alternative approach is often used, both in actual one-dimensional thermal transport problems and for interfaces: computing the total phonon transmission at each frequency as the trace of a matrix product, which is independent from the basis used for the wave functions. An example of such approach is the Caroli (or three-region) formula. The phonon transmission is enough to calculate the total thermal conductance in the Landauer–Büttiker formalism and can be used to implement semi-detailed interface scattering in the context of device simulations.33 However, it does not contain enough information to understand how phonon modes at either side are connected to one another. Fortunately, a systematic approach to the full calculation has emerged in recent years,34 which we implemented and used for the present work. In this method, the phonon modes in the leads at a given frequency are extracted from the solutions to ordinary eigenvalue problems for the so-called Bloch matrices, each of which is defined as the product of a Green's function times a hopping IFC matrix, i.e., the IFC block connecting adjacent segments of the leads. The group velocities and transmission are then computed through further matrix operations, and the propagating character of each mode is only analyzed once, so that inconsistencies between different steps of the calculation are avoided.34
We consider the transport direction to be the cubic [111] crystal axis, parallel to the hexagonal [0001] axis, and hence the layer stacking order goes from ABCABC in the ZB to ABABAB in the WZ. This crystallographic axis corresponds to the growth orientation of NWs featuring crystal-phase junctions and twin boundaries reported in the experiments. We take the z axis as perpendicular to the interface and therefore parallel to the transport direction. We study bulk junctions as an effective model of NWs with realistic diameters, which have bulk-like phonon dispersions. Phonon scattering by the sidewalls of the NW, which to a first approximation is independent on the interface scattering that we study here, is thus not explicitly considered.
The general workflow, in the transport geometry described above, is as follows. We obtain a full set of harmonic IFCs for a supercell containing two mirror copies of the interface with the Phonopy code.35 The supercell is long enough that interactions between those two copies are negligible, and the regions in between far from the interfaces can be taken as models for the bulk leads. From this set we extract all the IFC submatrices required to describe the harmonic interactions between the scattering regions and the leads, as well as the on-site and hopping matrices of the infinite leads themselves. We enforce the acoustic sum rules on these constants using a Lagrange-multiplier approach to eliminate spurious numerical effects and strictly enforce the short-sighted character of the interactions in this picture. The next step is to sample the two-dimensional BZ of parallel wave vectors using a 61 × 61 regular grid. For each q∥ we apply the procedure described in ref. 34. The ingredients are the effective complex IFCs computed by Fourier-transforming the real-space IFCs along the parallel directions for that value of q∥, the Green's function of the scattering region obtained by direct inversion, and the surface Green's functions of the leads calculated by a real-space renormalization-group method, decimation.36,37 With those ingredients we obtain the mode-by-mode energy transmission coefficients between phonons at both sides of the interface34 and, by summing over channels, the total phonon transmission for a particular q∥. An average over the two-dimensional parallel BZ is required to compute the total transmission determining the thermal conductance. The final step is to isolate the effect of the interface, whose thermal resistance per unit area is extracted using the standard approach of subtracting the contribution of the leads:
![]() | (1) |
The thermal conductances in eqn (1) are computed using the Landauer formula38,39 as
![]() | (2) |
We also perform NEMD calculations on GaP and InP boxes containing a crystal-phase interface that separates a ZB and a WZ segment of 40 nm each with the LAMMPS code.40 After an initial structural relaxation with a standard conjugate gradient algorithm we connect the ends of the computational cell to two Nosé–Hoover thermostats set at a hot and a cold temperature, TH and TC. We then evaluate the resulting heat flux from the work done by the thermostats.41,42 The total simulation time is 2.5 ns and we use a timestep of 1 fs. The heat flux and temperature are averaged over the last nanosecond, after checking that the system reached the nonequilibrium steady state prior to that. The temperature as a function of position, T(z), is computed as 2Ek/(3kB), where Ek is the average kinetic energy per atom within a slice of thickness δz of the computational cell. The TBR is calculated as the ratio between the temperature discontinuity across the interface (ΔT) and the heat current through it (J), as well as with the nonequilibrium thermodynamics formalism of ref. 43; the two approaches yield values in good agreement. We use a fixed ΔT = TH − TC of 100 K, similarly to other NEMD calculations,42,44,45 but vary TH and TC in order to sample different interface temperatures.
With the purpose to understand the difference between the two crystal phase interfaces, we analyze the mismatch in vibrational properties between the ZB and the WZ phase of the two materials. In principle, a situation where the (E) of the left and right leads are as similar as possible favors phonon propagation across the interface, since energy can be more easily conserved. To check whether this intuitive fact correlates with the thermal transport results, in Fig. 2 we plot the absolute value of the difference between the transmission of the ZB and the WZ crystal-phase,
ZB(E) and
WZ(E), for both GaP and InP. As can be seen from looking directly at the curves or at their integrals (depicted in the bottom panel of the same figure), the differences between the ZB and WZ transmissions are mostly concentrated in the low-frequency range in the case of InP. Hence they are indeed expected to result in a larger resistance, since this range is especially critical for thermal transport because of the higher group velocities and lower scattering rates of phonons with such energies. This observation agrees well with the calculated TBR of Fig. 1(c). It is also related to the significantly larger ratio between the thermal conductivities of the cubic and hexagonal phases in InP with respect to GaP, as discussed in detail in ref. 8. Nonetheless, it is already clear that this qualitative analysis will not be applicable to the case of twin boundaries where, no matter the compound, the transmissions of the perfect leads are identical by construction.
![]() | ||
Fig. 2 (Top) Absolute value of the difference between the transmission of the ZB and the WZ leads, Δ![]() |
Within the NEGF scheme, the interface is described at the harmonic level and phonon–phonon scattering processes are neglected. In a bulk material anharmonicities result in resistive collisions, i.e. Umklapp scattering, and are the main cause for the finite value of the thermal conductivity in solids. For an interface, conversely, anharmonic processes can increase the conductance by opening new transmission channels, analogously to the case of inelastic electron tunneling spectroscopy.47,48 This is easily illustrated considering the limiting case in which the leads have misaligned phonon gaps: a phonon of energy ħω1 of one of the leads cannot be transmitted across the interface if the other lead has no available states at that energy. With anharmonicities in play, on the other hand, three-phonon process such as ħω1 + ħω2 → ħω3 and ħω1 → ħω2 + ħω3 become possible at the interface and increase the probability of transmission. Our NEMD simulations, which naturally account for anharmonicities of all orders, allow us to quantify the influence of the effect just described by comparing the results with the NEGF calculations carried out using the same potentials. It must be kept in mind, however, that within NEMD the atoms move classically following Newton's equations of motion, so in the low-temperature regime where the effects of quantum statistics are relevant these results should only be taken qualitatively.
Such a comparison is presented in Fig. 1(c). The good agreement between NEGF and NEMD for the GaP crystal-phase junction indicates that in this case anharmonicities at the interface do not play an important role and the interface thermal resistance decreases only marginally when they are included. A similar agreement has been recently reported for PbTiO3 domain walls.49 The behavior of the InP crystal-phase interface, on the other hand, is quite different: at T > 300 K we find that the TBR is five times lower than the prediction of the harmonic approximation and of the same order of the one of GaP. The much larger significance of anharmonicity in InP can be probably traced to the higher degree of mismatch between the ZB and WZ phases at the harmonic level in the regions most relevant for thermal transport, as discussed in previous paragraphs based on the results presented in Fig. 2.
We find very low TBRs for both GaP and InP twin interfaces, closely resembling the results obtained for the crystal-phase interface of GaP. Specifically, the values of the TBR at T > 100 K are approximately 5 × 10−10 m2 K W−1 (see Fig. 3). Besides the relevance that they have for applications, the results obtained with the twin interfaces are important because they directly challenge common phenomenological models for the calculation of the TBR such as the acoustic mismatch model (AMM) or the diffuse mismatch model (DMM), either within the Debye approximation56 or in the full-band formalism.57,58 Both of these models derive the TBR from the mismatch of the vibrational properties of the materials involved. In a twin boundary, however, the phonon dispersions on the two sides of the interface are the same. Therefore the AMM would predict a zero TBR, in disagreement with our results, which yield a small, but finite value. Incidentally, we observe that the twin homojunction is the prototype of a flat interface, where the AMM should provide more accurate results. The DMM, on the other hand, would yield a nonzero TBR, but it would have no way to discriminate between a twin interface and a virtual interface with truly zero resistance, i.e., one separating two halves of a fully homogeneous sample considered as different regions. This is a known computational artifact of the DMM.
For a closer look at the scattering properties of the twin boundary, we have analyzed some local properties of the interface itself. In Fig. 4 we plot the average nearest-neighbor distance, the on-site and the nearest-neighbor harmonic IFCs moving across the interface. As expected, in the case of the crystal-phase interface these magnitudes change from a ZB to a WZ bulk value and the spatial extent of this transition is a good estimate of the thickness of the interface.43,59,60 In the case of the twin boundary the bulk values on either side are the same, but all these magnitudes feature a bump at the interface, an indication of the small, but non-negligible relaxation of the atomic structure. Phase field crystal methods have been shown to be able to capture the effects of a finite interface thickness and of the related strain in semiconductor systems.61 We argue that our results can be used to inform one of such phase field models, which can then be used for computationally lighter simulations.
To get additional insight into the microscopic origin of the thermal resistance at a twin interface, we also analyze the detailed transmission function in twinned InP. There is a clear one-to-one correspondence between phonon modes at both sides, and indeed the overall phonon transmissions of both leads are exactly identical. However, although the 2D Brillouin zones from which q∥ is sampled are equivalent, they are not identical because the respective mode polarization vectors are rotated by a 60° angle. This is well illustrated by Fig. 5(a), where we plot the phonon dispersion of both leads for a specific low-symmetry q∥ value, (qx, qy) = (cos(π/8), sin(π/8)). Since q∥ is conserved in the scattering process, this becomes a one-dimensional scattering problem between effective leads with different phonon spectra due to the rotation of the mode polarizations. The less-than-perfect transmission found for this value of q∥, which we present in Fig. 5(b), is a clear reflection of this fact. If we look at how modes incident on the interface from the left can be transmitted to the right lead, it becomes even more apparent that, even when there is only one vibrational mode available at each side at a particular energy, the coupling between them is far from ideal. To a certain extent this can be understood in terms of the overlap between the rotated phonon polarizations.
These two effects –the structural distortion localized at the interface and the rotation of the polarization vectors– are the physical reasons behind the TBR of the twin boundary. It is useful to compare it with the TBR of a simple stacking defect such as ABCABC|AB|ABCABC, i.e. an all-ZB system with a missing layer. In this interface there is no rotation of the polarization vectors between the leads, which are identical, so all the scattering must come only from the local distortion at the interface. The resulting TBR, plotted in Fig. 3 for the case of InP, is slightly larger than the one of the twin boundary, indicating that in this case a larger local distortion (the interface region is thicker than in the twin boundary) can have, on its own, an effect comparable to the more complex scattering at a twin boundary. Indeed, these results are a useful reminder of the fact that, while it is certainly true that homointerfaces are flat and have no chemical intermixing, they do nonetheless have an effective thickness, and that distances between atomic planes (and thus interatomic forces) are altered within a finite, if narrow, spatial region.
The transmission and TBR for a crystal-phase interface between 3C- and 2H-Si and for a twin boundary are shown in Fig. 6. We obtain a remarkably low thermal resistance associated to both the twin boundary and the crystal-phase interface. The values are only slightly lower than those obtained for GaP and InP in previous sections, which can be somewhat surprising given that the additional contribution to interface mismatch coming from the misalignment of polar In/Ga–P bonds is obviously absent in Si. It must be taken into account, however, that the level of ionicity of GaP and InP is relatively low, with those bonds being mostly covalent. The similar values allow us to reasonably speculate that the geometrical mechanism behind the TBR (interfacial readjustment and rotation of the 2D Brillouin zone plus the polarization vectors) have a certain degree of universality.
These results have been obtained with a carefully optimized NN potential that gives phonon dispersions in excellent agreement with DFT. For comparison we also calculate (E) and the TBR using the IFCs computed with the Tersoff potential.30 As can be seen, neither the
(E) of any of the leads nor that of a system with an interface are well reproduced by the Tersoff potential, which in particular overestimates the phonon bandwidth by as much as 10 meV. These results point to the need of a well tested potential that provides phonon dispersions in good agreement with DFT benchmarks in order to obtain reliable results for the phononic properties. Nevertheless, we note that, due to a cancellation of errors in the three
(E), the estimate of the TBR computed from the Tersoff potential is reasonably accurate, reinforcing the idea of some relatively universal geometric effects in play as stated above.
![]() | ||
Fig. 7 Transmission through single and double InP twin boundaries for a parallel wave vector of (qx, qy) = (cos(π/8), sin(π/8)). |
As can be seen, the addition of a second barrier does not necessarily imply a marked decrease of the transmission. This is a known behavior in the harmonic regime in the limit of large reflection coefficients: if a phonon is efficiently stopped by an interface, adding a second identical interface will not have a marked effect on the overall transmission. More importantly, however, we observe as well that several antiresonances appear throughout the transmission spectrum, an evidence that the two interfaces give rise to destructive interference between the incoming and the reflected waves. These results suggest that if we added more barriers, i.e. twin boundaries, the antiresonances would result in forbidden bandgaps, just like it is observed with heterointerfaces. In other words, on a qualitative level twinning superlattices can behave like conventional superlattices, as also hinted by the NEMD simulations of Porter and coworkers.75
However, these considerations are valid only in the purely coherent and harmonic regime, where it is assumed that phonons retain their coherence while traveling from one interface to the next. As a matter of fact, it has been demonstrated76 that, as far as thermal transport is concerned, phonon scattering by different periods of a superlattice can be treated as incoherent. This assertion is strongly supported by the accuracy and predictive value afforded by calculations based on that hypothesis. Moreover, regarding forbidden regions, the realization of phononic metamaterials for thermal transport applications is hindered by the extremely long wavelength of the relevant phonons. Whether twinning superlattices can mitigate these effects is still an open question. Thus, we argue that a more realistic description of heat transport in twinning superlattices with multiple interfaces calls for a theoretical approach where phonons are scattered elastically at the interfaces, according to the detailed transmission discussed in this work, and predominantly inelastically in between.33 Experimental measurements will also be crucial to shed light on these issues.
Footnotes |
† Electronic supplementary information (ESI) available. See DOI: 10.1039/C9NR05274G |
‡ These authors contributed equally to this work. |
This journal is © The Royal Society of Chemistry 2019 |