Nanoscale This article is licensed under a Creative Commons Attribution-NonCommercial 3.0 Unported Licence. View Article Online PAPER View Journal | View Issue Phonon transport across crystal-phase interfaces and twin boundaries in semiconducting Open Access Article. Published on 12 August 2019. Downloaded on 7/21/2021 11:30:42 AM. Cite this: Nanoscale, 2019, 11, 16007 nanowires† 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 Mingod and Riccardo Rurali *b We combine state-of-the-art Green’s-function methods and nonequilibrium molecular dynamics calcu- lations 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 Received 21st June 2019, mechanisms that give rise to phonon scattering in these interfaces, quantify their thermal boundary resis- Accepted 9th August 2019 tance and illustrate the failure of common phenomenological models in predicting those features. In par- DOI: 10.1039/c9nr05274g ticular, we show that twin boundaries have a small but finite interface thermal resistance that can only be rsc.li/nanoscale understood in terms of a fully atomistic picture. Recent advances in nanowire (NW) growth have provided small, and (iv) they cannot feature any chemical intermixing. access to many crystal phases of semiconductor compounds Therefore, there is a fundamental interest in the properties of that either are precluded in bulk form or require extreme these novel interfaces, particularly concerning the role they temperature and pressure conditions to become stable.1 The play in transport phenomena. high degree of control of the growth parameters that favor one Phonons are scattered at the interface between two polytype over another allows creating homojunctions, most different materials, giving rise to the phenomenon of thermal often between a zinc blende (ZB) and a wurtzite (WZ) crystal. boundary resistance (TBR). The TBR is due to the breakdown These homojunctions form the basis of crystal-phase superlat- of translational symmetry by the interface, and the difference tices, a novel type of structure where different crystal phases of in the mass of the constituent atomic species plays a promi- the same material, rather than different materials, are alter- nent role. Accordingly, one would expect that phonons could nated with a given periodicity.2,3 An even subtler type of homo- pass through a crystal-phase interface, where the atomic mass junction is the twin boundary, which results from a 60° does not vary, without being scattered. On the other hand, it rotation of part of the ZB crystal structure that, when arranged has been recently shown that the thermal conductivity can periodically, gives rise to the so-called twinning decrease by as much as 40% from a cubic to a hexagonal poly- superlattices.4,5 These homojunctions are special in several type of the same material.6–8 Therefore, there is a need for aspects: (i) they are atomically flat and do not exhibit rough- detailed insight from atomistic calculations of heat transport ness, (ii) they are defect-free, (iii) their lattice mismatch is very 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, a phonon scattering at crystal-phase interfaces and, particularly, Institute of Materials Chemistry, TU Wien, A-1060 Vienna, Austria b Institut de Ciència de Materials de Barcelona (ICMAB–CSIC), Campus de Bellaterra, twin boundaries defies many phenomenological models used 08193 Bellaterra, Barcelona, Spain. E-mail: rrurali@icmab.es to account for interface thermal resistance; on the other hand, c Departament d’Enginyeria Electrònica, Universitat Autònoma de Barcelona, 08193 crystal-phase engineering can be used to design materials with Bellaterra, Barcelona, Spain tailor-made properties3,11–15 and twinning superlattices are d CEA, LITEN, 17 Rue des Martyrs, 38054 Grenoble, France playing an increasingly important role in nanowire science.16–19 † Electronic supplementary information (ESI) available. See DOI: 10.1039/ C9NR05274G In this work we use both state-of-the-art nonequilibrium ‡ These authors contributed equally to this work. Green’s functions (NEGF) and nonequilibrium molecular This journal is © The Royal Society of Chemistry 2019 Nanoscale, 2019, 11, 16007–16016 | 16007 View Article Online Paper Nanoscale dynamics (NEMD), based on carefully parameterized intera- system is expressed as a sum of contributions from the P tomic potentials, to calculate the TBR of a few prototypical environments of the different atoms, E ¼ Ei . An atom’s i homojunctions, compare them with the values obtained in environment is defined by the collection of distances, {rij}, similar, conventional heterojunctions and elucidate the This article is licensed under a Creative Commons Attribution-NonCommercial 3.0 Unported Licence. between atom i and its neighbors, up to a cutoff radius rc. This atomic-scale mechanisms that originate interface scattering. collection of distances is transformed into an array of descrip- These computational techniques provide complementary P   2 information: while NEGF accounts for the quantum nature of tors, Dil ¼ f rij ; rc eσl jrij μl j , for the environment of each j phonons in the low to mid temperature regime, NEMD natu- atom i, where (μl, σl ) are adjustable parameters for each rally accounts for anharmonicity at all orders, which becomes    1 πr more important at high temperature. member l of the array, and f ðr; rc Þ ¼ 1 þ cos is a Open Access Article. Published on 12 August 2019. Downloaded on 7/21/2021 11:30:42 AM. 2 rc smooth cutoff function. The contribution Ei is given by a neural network that has Dil as inputs, and the forces ~ Fk on Computational methods each atom k are the anaytical derivatives of E with respect to atom k’s coordinates. The network weights, as well as the sets Ab initio calculations based on density functional theory (DFT) of parameters (μl, σl ), are adjusted to minimize a cost function, have demonstrated their predictive power in many areas of P s P s condensed matter physics,20 including the vibrational and Δ ¼ c  ðENN  EDFT s Þ2 þ j~ Fi;DFT  ~s Fi;NN j2 Þ, of the differ- s i thermal properties of materials.21,22 Nevertheless, because of ences between the network- and DFT-calculated energies and the high computational load involved, sometimes their use is forces of all the training systems, s. This minimization is per- not practical. This is often the case for the calculation of the formed by stochastic gradient descent and backpropagation interatomic force constants (IFCs) needed to obtain the using tensorflow’s32 GPU implementation. We used 32 phonon dispersion and the thermal conductivity in systems different values of μ and α, and neural network with 2 hidden featuring an interface, which normally need large simulation layers, i.e. 3 layers of coefficients, with 200 and 100 units, c = cells for their description. Lattice properties can also be 0.01, obtaining a RMSE of 0.05 eV Å−1 for forces and 1.3 meV computed using computationally lighter, but still reliable per atom for energies. methods, such as carefully parameterized interatomic Even in the presence of a periodic reconstruction, two dis- potentials.23,24 A common drawback of these methods, crete translation symmetries can be identified along directions however, is their limited transferability, and their use in a parallel to an interface between two crystalline regions. different environment from the one used for the parameteriza- Therefore, elastic scattering of phonons at the interface is tion is in principle questionable. Polytypism is a textbook case subject to the conservation of a two-dimensional parallel in this sense, and indeed it is not guaranteed that a potential momentum q∥. The problem of heat-carrier reflection and parameterized to reproduce the properties of a ZB semi- transmission by an interface can thus be decoupled into an conductor works equally well with the WZ crystal. infinite set of one-dimensional scattering problems, one for To address this concern in binary semiconductors, in this each possible value of q∥ in a two-dimensional parallel article we use a many-body interatomic potential that builds Brillouin zone, BZ∥. The full solution to each of those pro- on the previous work of Vashishta and coworkers,25,26 and blems consists in obtaining the detailed transmission prob- whose parameterization for GaP27 and InP28 reproduces the abilities per unit time among phonon modes on the left- and properties of different crystal phases, including ZB, WZ, and rock salt. As we are mostly interested in the vibrational pro- right-hand sides of the interface, constrained by the conserva- perties, we calculated the phonon dispersions of the ZB tion of energy. For instance, for a mode incident upon the and WZ crystals and compared them with those obtained by interface from the left lead with parallel and perpendicular first-principles DFT calculations. We found that the ZB components of the momentum (q∥, q⊥) and angular frequency and the WZ phase of both GaP and InP are indeed described ω, one first has to search for all the phonon modes in the left to a similar level of accuracy. Acoustic branches, whose and the right sides compatible with both q∥ and ω, and then accuracy is important in the study of interface or dislocation calculate the amplitude of each allowed scattering process. scattering,29 are reasonably well reproduced (see ESI†). There This involves exploring an effective one-dimensional Brillouin is a certain overestimation of the group velocity of the trans- zone BZ⊥(q∥), looking for values of q⊥ and phonon branch verse modes in the case of GaP and InP, but it is similar for indices yielding the same frequency as the incident phonon. both polytypes. The main difference with respect to a typical one-dimensional In contrast, common interatomic force fields for Si, such as problem involving scattering by a barrier is the lack of time- the Tersoff30 or the Stillinger–Weber potentials,31 give a rather reversal symmetry, which only applies to the line crossing the poor description of vibrational properties.29 Therefore, we Γ point. used machine-learning techniques to build a neural-network The detailed procedure described above is computationally (NN) potential for Si that enables phonon dispersions with expensive and has a tendency to numerical instability because ab initio quality for both the cubic diamond (3C) and the hex- of the need to solve the inverse problem of finding wave agonal diamond (lonsdaleite, 2H) polytypes. The energy of the vectors as functions of frequency. Therefore, an alternative 16008 | Nanoscale, 2019, 11, 16007–16016 This journal is © The Royal Society of Chemistry 2019 View Article Online Nanoscale Paper approach is often used, both in actual one-dimensional Green’s functions of the leads calculated by a real-space renor- thermal transport problems and for interfaces: computing the malization-group method, decimation.36,37 With those ingredi- total phonon transmission at each frequency as the trace of a ents we obtain the mode-by-mode energy transmission coeffi- matrix product, which is independent from the basis used for cients between phonons at both sides of the interface34 and, This article is licensed under a Creative Commons Attribution-NonCommercial 3.0 Unported Licence. the wave functions. An example of such approach is the Caroli by summing over channels, the total phonon transmission for (or three-region) formula. The phonon transmission is enough a particular q∥. An average over the two-dimensional parallel to calculate the total thermal conductance in the Landauer– BZ is required to compute the total transmission determining Büttiker formalism and can be used to implement semi- the thermal conductance. The final step is to isolate the effect detailed interface scattering in the context of device simu- of the interface, whose thermal resistance per unit area is lations.33 However, it does not contain enough information to extracted using the standard approach of subtracting the con- Open Access Article. Published on 12 August 2019. Downloaded on 7/21/2021 11:30:42 AM. understand how phonon modes at either side are connected to tribution of the leads: 1  1  one another. Fortunately, a systematic approach to the full cal- culation has emerged in recent years,34 which we implemented Ri ¼ G1 joint system  Gleft þ G1 right : ð1Þ 2 and used for the present work. In this method, the phonon modes in the leads at a given frequency are extracted from the The thermal conductances in eqn (1) are computed using solutions to ordinary eigenvalue problems for the so-called the Landauer formula38,39 as Bloch matrices, each of which is defined as the product of a ð   1 @n0 ðω; T Þ Green’s function times a hopping IFC matrix, i.e., the IFC GðTÞ ¼ E TðEÞ dE ð2Þ 2π ℏ Ω @T block connecting adjacent segments of the leads. The group velocities and transmission are then computed through where n0 is the equilibrium Bose–Einstein distribution, Ω is further matrix operations, and the propagating character of the area of the unit cell parallel to the interface, and ħ is the each mode is only analyzed once, so that inconsistencies reduced Planck constant. between different steps of the calculation are avoided.34 We also perform NEMD calculations on GaP and InP boxes We consider the transport direction to be the cubic [111] containing a crystal-phase interface that separates a ZB and a crystal axis, parallel to the hexagonal [0001] axis, and hence WZ segment of 40 nm each with the LAMMPS code.40 After an the layer stacking order goes from ABCABC in the ZB to initial structural relaxation with a standard conjugate gradient ABABAB in the WZ. This crystallographic axis corresponds to algorithm we connect the ends of the computational cell to the growth orientation of NWs featuring crystal-phase junc- two Nosé–Hoover thermostats set at a hot and a cold tempera- tions and twin boundaries reported in the experiments. We ture, TH and TC. We then evaluate the resulting heat flux from take the z axis as perpendicular to the interface and therefore the work done by the thermostats.41,42 The total simulation parallel to the transport direction. We study bulk junctions as time is 2.5 ns and we use a timestep of 1 fs. The heat flux and an effective model of NWs with realistic diameters, which have temperature are averaged over the last nanosecond, after bulk-like phonon dispersions. Phonon scattering by the side- checking that the system reached the nonequilibrium steady walls of the NW, which to a first approximation is independent state prior to that. The temperature as a function of position, on the interface scattering that we study here, is thus not expli- T (z), is computed as 2Ek/(3kB), where Ek is the average kinetic citly considered. energy per atom within a slice of thickness δz of the compu- The general workflow, in the transport geometry described tational cell. The TBR is calculated as the ratio between the above, is as follows. We obtain a full set of harmonic IFCs for a temperature discontinuity across the interface (ΔT ) and the supercell containing two mirror copies of the interface with heat current through it ( J), as well as with the nonequilibrium the Phonopy code.35 The supercell is long enough that inter- thermodynamics formalism of ref. 43; the two approaches actions between those two copies are negligible, and the yield values in good agreement. We use a fixed ΔT = TH − TC of regions in between far from the interfaces can be taken as 100 K, similarly to other NEMD calculations,42,44,45 but vary TH models for the bulk leads. From this set we extract all the IFC and TC in order to sample different interface temperatures. 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 Crystal-phase interfaces in GaP and enforce the acoustic sum rules on these constants using a InP Lagrange-multiplier approach to eliminate spurious numerical effects and strictly enforce the short-sighted character of the Our first target systems are crystal-phase junctions in GaP and interactions in this picture. The next step is to sample the two- InP. In Fig. 1 we plot the phonon transmission, T(E), as a dimensional BZ of parallel wave vectors using a 61 × 61 regular function of energy. The side panels show the same magnitude grid. For each q∥ we apply the procedure described in ref. 34. for the leads, i.e., for infinite, unperturbed all-ZB and all-WZ The ingredients are the effective complex IFCs computed by systems. From this energy-resolved result it is already possible Fourier-transforming the real-space IFCs along the parallel to draw the conclusion that the degradation of the trans- directions for that value of q∥, the Green’s function of the scat- mission induced by the interface is rather moderate in the tering region obtained by direct inversion, and the surface case of GaP, while it is more sizeable for InP, particularly in This journal is © The Royal Society of Chemistry 2019 Nanoscale, 2019, 11, 16007–16016 | 16009 View Article Online Paper Nanoscale This article is licensed under a Creative Commons Attribution-NonCommercial 3.0 Unported Licence. Open Access Article. Published on 12 August 2019. Downloaded on 7/21/2021 11:30:42 AM. Fig. 1 Transmission through a (a) GaP and (b) InP crystal-phase interface. The harmonic transmission functions of the perfect, infinite leads are shown in the side panels. The ball-and-stick models show the change from ABCABC stacking along the [111] cubic axis in the left lead to ABABAB stacking along the [0001] hexagonal axis in the right lead (blue, yellow, and red spheres represent Ga, P, and In atoms). (c) TBR as a function of inter- face temperature for the two crystal-phase interfaces obtained using the T(E) from panels (a) and (b). Data points show the TBR computed from NEMD, which includes anharmonic effects at the interface. The uncertainty in the estimate of the TBR from NEMD, arising from the fluctuating character of the temperature profile, is 9.6 × 10−11 m2 K W−1 and 5.9 × 10−10 m2 K W−1 for GaP and InP, respectively. the low-frequency part of the spectrum. These considerations are confirmed by the quantitative analysis of the TBR, which we compute from eqn (1) and display in Fig. 1(c). We find an extremely low TBR for the GaP crystal-phase junction, which at room temperature is approximately one order of magnitude smaller than in GaN–AlN46 or Si–Ge heterointerfaces43 (see also Table 1). The InP crystal-phase junction, on the other hand, is considerably more resistive and takes values in the same order of magnitude as those of the GaN–AlN system. 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 T(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 absol- Fig. 2 (Top) Absolute value of the difference between the transmission of the ZB and the WZ leads, ΔTZB/WZ(E), for GaP and InP plotted in ute value of the difference between the transmission of the ZB Fig. 1. (Bottom) Definite integral from E = 0 of the differences depicted and the WZ crystal-phase, TZB(E) and TWZ(E), for both GaP in the top panel. Table 1 TBR at 300 K of the homointerfaces studied in this work, com- and InP. As can be seen from looking directly at the curves or pared to the values of representative III–V and group-IV at their integrals (depicted in the bottom panel of the same heterointerfaces figure), the differences between the ZB and WZ transmissions are mostly concentrated in the low-frequency range in the case Interface R300K i [m2 K W−1] of InP. Hence they are indeed expected to result in a larger GaPphase 4.26 × 10−10 resistance, since this range is especially critical for thermal GaPtwin 3.54 × 10−10 InPphase 3.76 × 10−9 transport because of the higher group velocities and lower InPtwin 9.58 × 10−10 scattering rates of phonons with such energies. This obser- Siphase 2.89 × 10−10 vation agrees well with the calculated TBR of Fig. 1(c). It is also Sitwin 2.39 × 10−10 GaN–AlN46 3.33 × 10−9 related to the significantly larger ratio between the thermal Si–Ge43 3.77 × 10−9 conductivities of the cubic and hexagonal phases in InP with 16010 | Nanoscale, 2019, 11, 16007–16016 This journal is © The Royal Society of Chemistry 2019 View Article Online Nanoscale Paper respect to GaP, as discussed in detail in ref. 8. Nonetheless, it included. A similar agreement has been recently reported for is already clear that this qualitative analysis will not be appli- PbTiO3 domain walls.49 The behavior of the InP crystal-phase cable to the case of twin boundaries where, no matter the com- interface, on the other hand, is quite different: at T > 300 K pound, the transmissions of the perfect leads are identical by we find that the TBR is five times lower than the prediction This article is licensed under a Creative Commons Attribution-NonCommercial 3.0 Unported Licence. construction. of the harmonic approximation and of the same order of Within the NEGF scheme, the interface is described at the the one of GaP. The much larger significance of anharmonicity harmonic level and phonon–phonon scattering processes are in InP can be probably traced to the higher degree of mis- neglected. In a bulk material anharmonicities result in resis- match between the ZB and WZ phases at the harmonic level in tive collisions, i.e. Umklapp scattering, and are the main cause the regions most relevant for thermal transport, as discussed for the finite value of the thermal conductivity in solids. For in previous paragraphs based on the results presented in Open Access Article. Published on 12 August 2019. Downloaded on 7/21/2021 11:30:42 AM. an interface, conversely, anharmonic processes can increase Fig. 2. the conductance by opening new transmission channels, ana- logously to the case of inelastic electron tunneling spectroscopy.47,48 This is easily illustrated considering the lim- Twin boundaries in GaP and InP iting case in which the leads have misaligned phonon gaps: a phonon of energy ħω1 of one of the leads cannot be trans- Although twin boundaries are known to occur in bulk polycrys- mitted across the interface if the other lead has no available talline samples,50–53 their use in the context of NWs, where states at that energy. With anharmonicities in play, on the they can be arranged parallel to each other and with a well other hand, three-phonon process such as ħω1 + ħω2 → ħω3 defined periodicity4,5 gives rise to very special systems, and a and ħω1 → ħω2 + ħω3 become possible at the interface and few isolated reports54,55 hint at the possibility to extend such increase the probability of transmission. Our NEMD simu- arrangements to bulk systems. Here the crystal phase is the lations, which naturally account for anharmonicities of all same on the two sides of the interface, but the ZB structure orders, allow us to quantify the influence of the effect just undergoes a rotation of 60° and forms a twin boundary planar described by comparing the results with the NEGF calculations defect. carried out using the same potentials. It must be kept in We find very low TBRs for both GaP and InP twin interfaces, mind, however, that within NEMD the atoms move classically closely resembling the results obtained for the crystal-phase following Newton’s equations of motion, so in the low-temp- interface of GaP. Specifically, the values of the TBR at T > erature regime where the effects of quantum statistics are rele- 100 K are approximately 5 × 10−10 m2 K W−1 (see Fig. 3). vant these results should only be taken qualitatively. Besides the relevance that they have for applications, the Such a comparison is presented in Fig. 1(c). The good results obtained with the twin interfaces are important agreement between NEGF and NEMD for the GaP crystal-phase because they directly challenge common phenomenological junction indicates that in this case anharmonicities at the models for the calculation of the TBR such as the acoustic mis- interface do not play an important role and the interface match model (AMM) or the diffuse mismatch model (DMM), thermal resistance decreases only marginally when they are either within the Debye approximation56 or in the full-band Fig. 3 Transmission through a (a) GaP and (b) InP twin boundary interface. The harmonic transmission functions of the perfect, infinite leads are shown in the side panels. The ball-and-stick models show the change from ABCABC to CBACBA stacking along the [111] cubic axis in the left and right leads (blue, yellow, and red spheres represent Ga, P, and In atoms). (c) TBR as a function of interface temperature for the two twin boundary interfaces obtained from the T(E) of panels (a) and (b). This journal is © The Royal Society of Chemistry 2019 Nanoscale, 2019, 11, 16007–16016 | 16011 View Article Online Paper Nanoscale formalism.57,58 Both of these models derive the TBR from the interface, an indication of the small, but non-negligible relax- mismatch of the vibrational properties of the materials ation of the atomic structure. Phase field crystal methods have involved. In a twin boundary, however, the phonon dispersions been shown to be able to capture the effects of a finite inter- on the two sides of the interface are the same. Therefore the face thickness and of the related strain in semiconductor This article is licensed under a Creative Commons Attribution-NonCommercial 3.0 Unported Licence. AMM would predict a zero TBR, in disagreement with our systems.61 We argue that our results can be used to inform one results, which yield a small, but finite value. Incidentally, we of such phase field models, which can then be used for com- observe that the twin homojunction is the prototype of a flat putationally lighter simulations. interface, where the AMM should provide more accurate To get additional insight into the microscopic origin of the results. The DMM, on the other hand, would yield a nonzero thermal resistance at a twin interface, we also analyze the TBR, but it would have no way to discriminate between a twin detailed transmission function in twinned InP. There is a clear Open Access Article. Published on 12 August 2019. Downloaded on 7/21/2021 11:30:42 AM. interface and a virtual interface with truly zero resistance, i.e., one-to-one correspondence between phonon modes at both one separating two halves of a fully homogeneous sample con- sides, and indeed the overall phonon transmissions of both sidered as different regions. This is a known computational leads are exactly identical. However, although the 2D Brillouin artifact of the DMM. zones from which q∥ is sampled are equivalent, they are not For a closer look at the scattering properties of the twin identical because the respective mode polarization vectors are boundary, we have analyzed some local properties of the inter- rotated by a 60° angle. This is well illustrated by Fig. 5(a), face itself. In Fig. 4 we plot the average nearest-neighbor where we plot the phonon dispersion of both leads for a distance, the on-site and the nearest-neighbor harmonic IFCs specific low-symmetry q∥ value, (qx, qy) = (cos(π/8), sin(π/8)). moving across the interface. As expected, in the case of the Since q∥ is conserved in the scattering process, this becomes a crystal-phase interface these magnitudes change from a ZB to one-dimensional scattering problem between effective leads a WZ bulk value and the spatial extent of this transition is a with different phonon spectra due to the rotation of the mode good estimate of the thickness of the interface.43,59,60 In the polarizations. The less-than-perfect transmission found for case of the twin boundary the bulk values on either side are this value of q∥, which we present in Fig. 5(b), is a clear reflec- the same, but all these magnitudes feature a bump at the Fig. 5 (a) Phonon dispersion of both leads in an InP twin interface and Fig. 4 (a) Nearest-neighbor distance, (b) on-site IFC and (c) nearest- (b) transmission for (qx, qy) = (cos(π/8), sin(π/8)). (c) Detailed transmission neighbor IFC as a function of the position along the transport direction for the same q∥: the numbers of rows and columns are equal to the for a crystal-phase interface and a twin boundary in InP. The on-site IFC number of modes in each lead at each sampled frequency, annotated in is computed as the Frobenius norm of the Jacobian of the forces of red; the color of each square gives the energy transmission coefficient   @Fm  between each pair of modes according to the color code on the right- atom m with respect to its own displacement,    @rm  , where m is an F   hand side. Therefore, and for instance, at E = 1.84 meV there are no pro- @Fm  pagating modes, at E = 5.53 meV there are two modes on each side with In atom; similarly, the nearest-neighbor IFC is defined as    @rn  , where F rather low transmission, and at E = 13.82 meV there is a single mode n is a P atom and a first neighbor of m such that zm < zn. with near ideal transmission. 16012 | Nanoscale, 2019, 11, 16007–16016 This journal is © The Royal Society of Chemistry 2019 View Article Online Nanoscale Paper tion of this fact. If we look at how modes incident on the inter- wetting angle of the catalytic droplet during growth.62 In the face from the left can be transmitted to the right lead, it case of group-IV semiconductors, on the other hand, the becomes even more apparent that, even when there is only one mechanisms are different and interfaces between cubic (3C) vibrational mode available at each side at a particular energy, and hexagonal (2H) phases have been obtained through shear This article is licensed under a Creative Commons Attribution-NonCommercial 3.0 Unported Licence. the coupling between them is far from ideal. To a certain stress,18,63,64 at NW kinks65,66 or taking advantage of the cata- extent this can be understood in terms of the overlap between lytic properties of magnetite nanocrystals.67 Consequently, the the rotated phonon polarizations. level of control of the positioning of hexagonal segments is These two effects –the structural distortion localized at the considerably lower than what can be achieved in III–V NWs. interface and the rotation of the polarization vectors– are the Nevertheless, isolated and non-periodic cubic–hexagonal inter- physical reasons behind the TBR of the twin boundary. It is faces in Si and Ge NWs are observed68 and thus it is interest- Open Access Article. Published on 12 August 2019. Downloaded on 7/21/2021 11:30:42 AM. useful to compare it with the TBR of a simple stacking defect ing to study phonon transport across these interfaces as well. such as ABCABC|AB|ABCABC, i.e. an all-ZB system with a Similar considerations can be made for twin boundaries that missing layer. In this interface there is no rotation of the polar- have been reported in Si NWs69,70 and that can be positioned ization vectors between the leads, which are identical, so all with increasingly deterministic precision.71,72 the scattering must come only from the local distortion at the The transmission and TBR for a crystal-phase interface interface. The resulting TBR, plotted in Fig. 3 for the case of between 3C- and 2H-Si and for a twin boundary are shown in InP, is slightly larger than the one of the twin boundary, indi- Fig. 6. We obtain a remarkably low thermal resistance associ- cating that in this case a larger local distortion (the interface ated to both the twin boundary and the crystal-phase interface. region is thicker than in the twin boundary) can have, on its The values are only slightly lower than those obtained for GaP own, an effect comparable to the more complex scattering at a and InP in previous sections, which can be somewhat surpris- twin boundary. Indeed, these results are a useful reminder of ing given that the additional contribution to interface mis- the fact that, while it is certainly true that homointerfaces are match coming from the misalignment of polar In/Ga–P bonds flat and have no chemical intermixing, they do nonetheless is obviously absent in Si. It must be taken into account, have an effective thickness, and that distances between atomic however, that the level of ionicity of GaP and InP is relatively planes (and thus interatomic forces) are altered within a finite, low, with those bonds being mostly covalent. The similar if narrow, spatial region. 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) Crystal-phase interfaces and twin have a certain degree of universality. boundaries in Si These results have been obtained with a carefully optimized NN potential that gives phonon dispersions in excellent agree- The formation of WZ segments in III–V NWs is today well ment with DFT. For comparison we also calculate T(E) and understood, and is known to be directly correlated to the the TBR using the IFCs computed with the Tersoff potential.30 Fig. 6 Transmission through a Si (a) phase and (b) twin boundary interface. The harmonic transmission functions of the perfect, infinite leads are shown in the side panels. The ball-and-stick models show the change from ABCABC to CBACBA stacking along the [111] cubic axis in the transition from the left to the right lead. We plot the results obtained with both the optimized NN potential (continuous lines), that yields phonon dispersions with ab initio accuracy, and the Tersoff potential (dashedlines). (c) TBR as a function of interface temperature for the two interfaces obtained from the NN-potential results for T(E) of panels (a) and (b). This journal is © The Royal Society of Chemistry 2019 Nanoscale, 2019, 11, 16007–16016 | 16013 View Article Online Paper Nanoscale As can be seen, neither the T(E) of any of the leads nor that of appear throughout the transmission spectrum, an evidence a system with an interface are well reproduced by the Tersoff that the two interfaces give rise to destructive interference potential, which in particular overestimates the phonon band- between the incoming and the reflected waves. These results width by as much as 10 meV. These results point to the need suggest that if we added more barriers, i.e. twin boundaries, This article is licensed under a Creative Commons Attribution-NonCommercial 3.0 Unported Licence. of a well tested potential that provides phonon dispersions in the antiresonances would result in forbidden bandgaps, just good agreement with DFT benchmarks in order to obtain like it is observed with heterointerfaces. In other words, on a reliable results for the phononic properties. Nevertheless, we qualitative level twinning superlattices can behave like conven- note that, due to a cancellation of errors in the three T(E), the tional superlattices, as also hinted by the NEMD simulations estimate of the TBR computed from the Tersoff potential is of Porter and coworkers.75 reasonably accurate, reinforcing the idea of some relatively uni- However, these considerations are valid only in the purely Open Access Article. Published on 12 August 2019. Downloaded on 7/21/2021 11:30:42 AM. versal geometric effects in play as stated above. coherent and harmonic regime, where it is assumed that phonons retain their coherence while traveling from one inter- face to the next. As a matter of fact, it has been demonstrated76 Multiple interfaces and interference that, as far as thermal transport is concerned, phonon scatter- ing by different periods of a superlattice can be treated as inco- effects herent. This assertion is strongly supported by the accuracy and The discussion in this paper is formulated in terms of the predictive value afforded by calculations based on that hypoth- transmission of a single interface. Multiple interfaces can in esis. Moreover, regarding forbidden regions, the realization of principle create interference effects, and it has been specu- phononic metamaterials for thermal transport applications is lated that they could constitute the basis for phononic meta- hindered by the extremely long wavelength of the relevant materials with tunable vibrational properties19 and thermal phonons. Whether twinning superlattices can mitigate these conductivities.73,74 Although an exhaustive study of crystal- effects is still an open question. Thus, we argue that a more rea- phase and twinning superlattices is beyond the scope of this listic description of heat transport in twinning superlattices work, as a test on the method we have examined the inter- with multiple interfaces calls for a theoretical approach where ference effects in a double-interface setting. To this end, we phonons are scattered elastically at the interfaces, according to have considered phonon transport across two twin boundaries the detailed transmission discussed in this work, and predomi- separated by 27 In–P layers (∼9 nm) and in Fig. 7 we plot the nantly inelastically in between.33 Experimental measurements transmission T(E) of this double barrier for the same q∥ of will also be crucial to shed light on these issues. Fig. 5. In the same figure we also plot the same result for the case of the isolated twin boundary for comparison. As can be seen, the addition of a second barrier does not Conclusions necessarily imply a marked decrease of the transmission. This In summary, we have carried out an atomistic theoretical study is a known behavior in the harmonic regime in the limit of of phonon scattering at crystal-phase interfaces and twin large reflection coefficients: if a phonon is efficiently stopped boundaries, both of which are commonly observed in NWs by an interface, adding a second identical interface will not and constitute the building blocks for crystal-phase and have a marked effect on the overall transmission. More impor- twinning superlattices, respectively. We have computed and tantly, however, we observe as well that several antiresonances analyzed the TBR showing that (i) in twin boundaries it is much lower than in typical (i.e. heterogeneous) semiconductor interfaces; (ii) in crystal-phase interfaces it can be as low as in twin boundaries, but in general it depends on the mismatch between the phonon dispersions of the two polytypes involved and anharmonicity can lower it at higher temperatures; (iii) in the elusive case of the twin boundary, it originates from a com- bined effect of the local structural distortion at the interface and the rotation of the polarization vectors. We have showed that phenomenological models often used to estimate the TBR cannot provide a description of scattering in these subtle interfaces, not even at a qualitative level. Finally, we have discussed the implications for superlattices, showing that nothing prevents homointerfaces to host wave interference effects similarly to those predicted in heterointerfaces. Conflicts of interest Fig. 7 Transmission through single and double InP twin boundaries for a parallel wave vector of (qx, qy) = (cos(π/8), sin(π/8)). There are no conflicts to declare. 16014 | Nanoscale, 2019, 11, 16007–16016 This journal is © The Royal Society of Chemistry 2019 View Article Online Nanoscale Paper Acknowledgements 16 I.-J. Chen, S. Lehmann, M. Nilsson, P. Kivisaari, H. Linke, K. A. Dick and C. Thelander, Nano Lett., 2017, 17, 902– We acknowledge financial support by the Ministerio de 908. Economía, Industria y Competitividad (MINECO) under grants 17 Y. Qian and Q. Yang, Nano Lett., 2017, 17, 7183–7190. This article is licensed under a Creative Commons Attribution-NonCommercial 3.0 Unported Licence. FEDER-MAT2017-90024-P, TEC2015-67462-C2-1-R (MINECO/ 18 C. Fasolato, M. De Luca, D. Djomani, L. Vincent, FEDER) and RTI2018-097876-B-C21 (MICINN/FEDER), the C. Renard, G. Di Iorio, V. Paillard, M. Amato, R. Rurali and Severo Ochoa Centres of Excellence Program under I. Zardo, Nano Lett., 2018, 18, 7075–7084. Grant SEV-2015-0496 and by the Generalitat de Catalunya 19 M. De Luca, C. Fasolato, M. A. Verheijen, Y. Ren, under grants no. 2014 SGR 301 and 2017 SGR M. Y. Swinkels, S. Kölling, E. P. A. M. Bakkers, R. Rurali, 1506. M. L. S. is funded through a Juan de la Cierva X. Cartoixà and I. Zardo, Nano Lett., 2019, 19, 4702–4711. Open Access Article. Published on 12 August 2019. Downloaded on 7/21/2021 11:30:42 AM. fellowship. M. R. M. acknowledges support from FPU2016/ 20 R. M. Martin, Electronic structure: basic theory and practical 02565. We thank J. A. Seijas-Bellido for preliminary tests with methods, Cambridge Univ. Press, Cambridge, 2004. the Vashishta potential. 21 W. Li, J. Carrete, N. A. Katcho and N. Mingo, Comput. Phys. Commun., 2014, 185, 1747–1758. 22 J. Carrete, B. Vermeersch, A. Katre, A. van Roekeghem, References T. Wang, G. K. Madsen and N. Mingo, Comput. Phys. Commun., 2017, 220, 351–362. 1 P. Caroff, J. Bolinsson and J. Johansson, IEEE J. Sel. Top. 23 D. Frenkel and B. Smit, Understanding Molecular Quantum Electron., 2011, 17, 829–846. Simulations: From Algorithms to Applications, Academic 2 P. Caroff, K. A. Dick, J. Johansson, M. E. Messing, Press, Inc., 1996. K. Deppert and L. Samuelson, Nat. Nanotechnol., 2009, 4, 24 A. Rohskopf, H. R. Seyf, K. Gordiz, T. Tadano and A. Henry, 50–55. npj Comput. Mater., 2017, 3, 27. 3 K. A. Dick, C. Thelander, L. Samuelson and P. Caroff, Nano 25 P. Vashishta, R. K. Kalia, J. P. Rino and I. Ebbsjö, Phys. Rev. Lett., 2010, 10, 3494–3499. B: Condens. Matter Mater. Phys., 1990, 41, 12197–12209. 4 R. E. Algra, M. A. Verheijen, M. T. Borgström, L.-F. Feiner, 26 I. Ebbsjö, R. K. Kalia, A. Nakano, J. P. Rino and G. Immink, W. J. P. van Enckevort, E. Vlieg and P. Vashishta, J. Appl. Phys., 2000, 87, 7708–7711. E. P. A. M. Bakkers, Nature, 2008, 456, 369–372. 27 C. I. Ribeiro-Silva, J. P. Rino, L. G. V. Gonçalves and 5 T. Burgess, S. Breuer, P. Caroff, J. Wong-Leung, Q. Gao, A. Picinin, J. Phys.: Condens. Matter, 2011, 23, 055801. H. Hoe Tan and C. Jagadish, ACS Nano, 2013, 7, 8105–8114. 28 P. S. Branicio, J. P. Rino, C. K. Gan and H. Tsuzuki, J. Phys.: 6 A. Togo, L. Chaput and I. Tanaka, Phys. Rev. B: Condens. Condens. Matter, 2009, 21, 095002. Matter Mater. Phys., 2015, 91, 094306. 29 T. Wang, J. Carrete, N. Mingo and G. K. H. Madsen, ACS 7 M. Raya-Moreno, H. Aramberri, J. A. Seijas-Bellido, Appl. Mater. Interfaces, 2019, 11, 8175–8181. X. Cartoixà and R. Rurali, Appl. Phys. Lett., 2017, 111, 032107. 30 J. Tersoff, Phys. Rev. B: Condens. Matter Mater. Phys., 1989, 8 M. Raya-Moreno, R. Rurali and X. Cartoixà, 2019, arXiv: 39, 5566–5568. 1901.03268. 31 F. H. Stillinger and T. A. Weber, Phys. Rev. B: Condens. 9 A. Weathers, A. L. Moore, M. T. Pettes, D. Salta, J. Kim, Matter Mater. Phys., 1985, 31, 5262–5271. K. Dick, L. Samuelson, H. Linke, P. Caroff and L. Shi, MRS 32 M. Abadi, et al., TensorFlow: Large-Scale Machine Learning Proc., 2012, 1404, mrsf11–1404-w04–06. on Heterogeneous Systems, 2015; https://www.tensorflow.org/, 10 D. Majumdar, S. Biswas, T. Ghoshal, J. D. Holmes and Software available from tensorflow.org. A. Singha, ACS Appl. Mater. Interfaces, 2015, 7, 24679– 33 A. van Roekeghem, B. Vermeersch, J. Carrete and 24685. N. Mingo, Phys. Rev. Appl., 2019, 11, 034036. 11 C. Thelander, P. Caroff, S. Plissard, A. W. Dey and 34 Z.-Y. Ong, J. Appl. Phys., 2018, 124, 151101. K. A. Dick, Nano Lett., 2011, 11, 2424–2429. 35 A. Togo and I. Tanaka, Scr. Mater., 2015, 108, 1–5. 12 S. Assali, I. Zardo, S. Plissard, D. Kriegner, M. A. Verheijen, 36 F. Guinea, C. Tejedor, F. Flores and E. Louis, Phys. Rev. B: G. Bauer, A. Meijerink, A. Belabbes, F. Bechstedt, Condens. Matter Mater. Phys., 1983, 28, 4397–4402. J. E. M. Haverkort and E. P. A. M. Bakkers, Nano Lett., 2013, 37 W. Zhang, T. S. Fisher and N. Mingo, Numer. Heat Transfer, 13, 1559–1563. Part B, 2007, 51, 333–349. 13 B. Loitsch, J. Winnerl, G. Grimaldi, J. Wierzbowski, 38 L. G. C. Rego and G. Kirczenow, Phys. Rev. Lett., 1998, 81, D. Rudolph, S. Morkötter, M. Döblinger, G. Abstreiter, 232–235. G. Koblmüller and J. J. Finley, Nano Lett., 2015, 15, 7544– 39 N. Mingo, Phys. Rev. B: Condens. Matter Mater. Phys., 2006, 7551. 74, 125402. 14 M. Amato, T. Kaewmaraya, A. Zobelli, M. Palummo and 40 S. Plimpton, J. Comp. Physiol., 1995, 117, 1–19. R. Rurali, Nano Lett., 2016, 16, 5694–5700. 41 J. Lee, V. Varshney, A. K. Roy, J. B. Ferguson and 15 X. Cartoixà, M. Palummo, H. I. T. Hauge, B. L. Farmer, Nano Lett., 2012, 12, 3491–3496. E. P. A. M. Bakkers and R. Rurali, Nano Lett., 2017, 17, 42 R. Rurali, X. Cartoixà and L. Colombo, Phys. Rev. B: 4753–4758. Condens. Matter Mater. Phys., 2014, 90, 041408. This journal is © The Royal Society of Chemistry 2019 Nanoscale, 2019, 11, 16007–16016 | 16015 View Article Online Paper Nanoscale 43 R. Rurali, L. Colombo, X. Cartoixà, O. Wilhelmsen, 60 R. Dettori, C. Melis, X. Cartoixà, R. Rurali and L. Colombo, T. T. Trinh, D. Bedeaux and S. Kjelstrup, Phys. Chem. Chem. Adv. Phys.: X, 2016, 1, 246–261. Phys., 2016, 18, 13741–13745. 61 K.-A. Wu and P. W. Voorhees, Phys. Rev. B: Condens. Matter 44 P. K. Schelling, S. R. Phillpot and P. Keblinski, Phys. Rev. B: Mater. Phys., 2009, 80, 125408. This article is licensed under a Creative Commons Attribution-NonCommercial 3.0 Unported Licence. Condens. Matter Mater. Phys., 2002, 65, 144306. 62 F. Glas, J.-C. Harmand and G. Patriarche, Phys. Rev. Lett., 45 D. G. Cahill, W. K. Ford, K. E. Goodson, G. D. Mahan, 2007, 99, 146101. A. Majumdar, H. J. Maris, R. Merlin and S. R. Phillpot, 63 L. Vincent, G. Patriarche, G. Hallais, C. Renard, C. Gardès, J. Appl. Phys., 2003, 93, 793–818. D. Troadec and D. Bouchier, Nano Lett., 2014, 14, 4828– 46 C. A. Polanco and L. Lindsay, Phys. Rev. B, 2019, 99, 4836. 075202. 64 L. Vincent, D. Djomani, M. Fakfakh, C. Renard, B. Belier, Open Access Article. Published on 12 August 2019. Downloaded on 7/21/2021 11:30:42 AM. 47 B. C. Stipe, M. A. Rezaei and W. Ho, Science, 1998, 1732–1735. D. Bouchier and G. Patriarche, Nanotechnology, 2018, 29, 48 N. Lorente and M. Persson, Phys. Rev. Lett., 2000, 85, 2997– 125601. 3000. 65 N. Jeon, S. A. Dayeh and L. J. Lauhon, Nano Lett., 2013, 13, 49 J. A. Seijas-Bellido, C. Escorihuela-Salayero, M. Royo, 3947–3952. M. P. Ljungberg, J. C. Wojdeł, J. Íñiguez and R. Rurali, 66 Z. He, J.-L. Maurice, Q. Li and D. Pribat, Nanoscale, 2019, Phys. Rev. B, 2017, 96, 140101. 11, 4846–4853. 50 G. Gottstein, Physical Foundations of Materials Science, 67 S. Biswas, J. Doherty, D. Majumdar, T. Ghoshal, K. Rahme, Springer-Verlag, Berlin-Heidelberg, 2004. M. Conroy, A. Singha, M. A. Morris and J. D. Holmes, 51 C. Kimmer, S. Aubry, A. Skye and P. K. Schelling, Phys. Rev. Chem. Mater., 2015, 27, 3408–3416. B: Condens. Matter Mater. Phys., 2007, 75, 144105. 68 F. J. Lopez, E. R. Hemesath and L. J. Lauhon, Nano Lett., 52 S. Aubry, C. J. Kimmer, A. Skye and P. K. Schelling, Phys. 2009, 9, 2774–2779. Rev. B: Condens. Matter Mater. Phys., 2008, 78, 064112. 69 F. M. Davidson, D. C. Lee, D. D. Fanfair and B. A. Korgel, 53 R. Aghababaei, G. Anciaux and J.-F. Molinari, Appl. Phys. J. Phys. Chem. C, 2007, 111, 2929–2935. Lett., 2014, 105, 194102. 70 J. Arbiol, A. Fontcuberta i Morral, S. Estradé, F. Peiró, 54 H. Hibino and T. Ogino, Mater. Sci. Eng., B, 2001, 87, 214– B. Kalache, P. Roca i Cabarrocas and J. R. Morante, J. Appl. 221. Phys., 2008, 104, 064312. 55 A. Fissel, E. Bugiel, C. Wang and H. Osten, J. Cryst. Growth, 71 N. Shin, M. Chi, J. Y. Howe and M. A. Filler, Nano Lett., 2006, 290, 392–397. 2013, 13, 1928–1933. 56 E. T. Swartz and R. O. Pohl, Rev. Mod. Phys., 1989, 61, 605– 72 N. Shin, M. Chi and M. A. Filler, ACS Nano, 2013, 7, 8206– 668. 8213. 57 P. Reddy, K. Castelino and A. Majumdar, Appl. Phys. Lett., 73 M. Maldovan, Nature, 2015, 14, 209–217. 2005, 87, 211908. 74 M. Maldovan, Nat. Mater., 2015, 14, 667–674. 58 J. Larroque, P. Dollfus and J. Saint-Martin, J. Appl. Phys., 75 A. Porter, C. Tran and F. Sansoz, Phys. Rev. B: Condens. 2018, 123, 025702. Matter Mater. Phys., 2016, 93, 195431. 59 S. Kjelstrup and D. Bedeaux, Non-Equilibrium 76 J. Carrete, B. Vermeersch, L. Thumfart, R. R. Kakodkar, Thermodynamics of Heterogeneous Systems, World Scientific, G. Trevisi, P. Frigeri, L. Seravalli, J. P. Feser, A. Rastelli and Singapore, 2008, vol. 16. N. Mingo, J. Phys. Chem. C, 2018, 122, 4054–4062. 16016 | Nanoscale, 2019, 11, 16007–16016 This journal is © The Royal Society of Chemistry 2019