Cite this Article
Xu Run-Feng, Han Kui, Li Hai-Peng. Effect of isotope doping on phonon thermal conductivity of silicene nanoribbons: A molecular dynamics study. Chinese Physics B, 2018, 27(2): 026801  
Permissions
Effect of isotope doping on phonon thermal conductivity of silicene nanoribbons: A molecular dynamics study
Xu Run-Feng, Han Kui †, Li Hai-Peng ‡
School of Physical Science and Technology, China University of Mining and Technology (CUMT), Xuzhou 221116, China

 

† Corresponding author. E-mail: han6409@263.net haipli@cumt.edu.cn

Project supported by the National Natural Science Foundation of China (Grant Nos. 11504418 and 11447033), the Natural Science Fund for Colleges and Universities in Jiangsu Province, China (Grant No. 16KJB460022), and the Fundamental Research Funds for the Central Universities of CUMT, China (Grant No. 2015XKMS075).

Abstract

Silicene, a silicon analogue of graphene, has attracted increasing research attention in recent years because of its unique electrical and thermal conductivities. In this study, phonon thermal conductivity and its isotopic doping effect in silicene nanoribbons (SNRs) are investigated by using molecular dynamics simulations. The calculated thermal conductivities are approximately 32 W/mK and 35 W/mK for armchair-edged SNRs and zigzag-edged SNRs, respectively, which show anisotropic behaviors. Isotope doping induces mass disorder in the lattice, which results in increased phonon scattering, thus reducing the thermal conductivity. The phonon thermal conductivity of isotopic doped SNR is dependent on the concentration and arrangement pattern of dopants. A maximum reduction of about 15% is obtained at 50% randomly isotopic doping with30Si. In addition, ordered doping (i.e., isotope superlattice) leads to a much larger reduction in thermal conductivity than random doping for the same doping concentration. Particularly, the periodicity of the doping superlattice structure has a significant influence on the thermal conductivity of SNR. Phonon spectrum analysis is also used to qualitatively explain the mechanism of thermal conductivity change induced by isotopic doping. This study highlights the importance of isotopic doping in tuning the thermal properties of silicene, thus guiding defect engineering of the thermal properties of two-dimensional silicon materials.

PACS: 68.65.-k;44.10.+i;65.80.-g;66.70.-f
Keyword:silicene;phonon thermal conductivity;isotope doping;molecular dynamics simulations
1. Introduction

In 2004, Novoselov et al. found the most suitable thin flake to prepare graphene from graphite by manual selection.[1] Although graphene is a unique layered material offering many potential novel applications, its compatibility with the existing semiconductor industry still faces severe challenges. Silicene, the counterpart of graphene in the silicon realm has recently attracted great attention in both theory[2–5] and experiment[6,7] due to its expected compatibility with current silicon technology. In particular, the experimental realization of silicene[8,9] opens up possibilities for exploring its properties and potential applications. Unlike the flat configuration of graphene, free-standing silicene is not stable in a flat configuration, in which atoms form a slightly buckled monolayer structure because silicon atoms tend to bond in an sp2–sp3 mixing hybridization instead of a pure sp2 hybridization.[10]

Thermoelectric (TE) material is a kind of functional material that can convert heat energy directly into electric energy. The performance of TE material is measured by the figure of merit,

where S, σ, T, and κ are the Seebeck coefficient, electrical conductivity, absolute temperature, and thermal conductivity, respectively. A simple approach to efficiency improvement is to minimize the thermal conductivity or maximize the power factor S2σ. Many previous studies have shown that the thermal conductivity of free-standing silicene (~ 20 W/mK–70 W/mK[11–13]) is considerably smaller than that of graphene, indicating that silicene has a greater advantage in usage of TE materials. Finding effective ways to improve the efficiency has become a key issue in current TE research. Silicene possesses excellent electronic properties and low thermal conductivity, and is thus a potential material for TE applications. The key to improving the TE efficiency of silicene relies on opening a bandgap to enhance the thermopower and suppress the lattice thermal conductivity. Given their buckled atomic structure,[14] silicene nanoribbons (SNRs) have a nonzero energy gap, which can be tuned further by applying external transverse electric fields and doping.[15] A few strategies have also been used to reduce thermal transport in SNRs, including adding substrates,[16,17] surface functionalization,[18] defect,[11,19,20] and strain.[21] These studies provide significant guidance for experimental realization.

Defects or impurities are ubiquitous in natural materials. As working with defect-free or impurity-free materials is almost impossible, understanding how defects and impurities change the electronic and thermal properties of systems is essential.[22] Yang et al. theoretically investigated the thermal conductivity of isotope-doped silicon nanowires and found that thermal conductivity can be reduced exponentially by isotopic defects at room temperature.[23] Yu et al. reported the enhancement of TE properties by modulating doping in silicon germanium alloy nanocomposites.[24] Sevincli et al. introduced14C isotopes to suppress phonon transport in graphene.[25] Ni et al. studied the tunable band gap and doping type in silicene by surface adsorption.[26] Recently, the reductions of thermal conductivity of silicene through germanium doping and isotope doping were also reported.[27, 28]

Isotope impurities provide a powerful technique for studying the phonon-related properties of nanomaterials. Technologically speaking, isotope doping easily introduces phonon-defect scattering in a silicene system without damaging the electronic quality, thus attracting the intensive interest of researchers. However, the influences of isotope doping, including doping concentration, doping pattern, and doping type on the thermal conductivity of SNRs, remain clear. In light of the extreme sensitivity of electron transport in nanostructures to isotope doping,[29] a parallel, simultaneous investigation of phonon transport is clearly needed. In this work, we introduce the isotope to tailor the phonon thermal conductivity of SNRs with (i) random atomic distribution of isotopes and (ii) superlattice-structured isotope substitution. Phonon thermal conductivity is calculated using reverse non-equilibrium molecular dynamics (RNEMD) simulation. The effects of size and isotope doping on the thermal conductivity are discussed. Phonon density of states (PDOS) analysis is also carried out to explain the mechanism. Our study can guide defect engineering of the TE properties of SNR materials.

Fig. 1. (color online) (a) Schematic diagram of the simulation model for the RNEMD. A small amount of heat is added into the hot region and removed from the two cold regions to create the heat fluxes. Periodic boundary conditions are applied in the x and y directions. (b) Typical temperature profile in an SNR along the x direction at 300 K.
2. Simulation methods

The RNEMD method[30] implemented in an open-source MD simulator LAMMPS[31] was utilized to impose a temperature gradient and obtain the values of thermal conductivity of SNRs. The key idea of the RNEMD was to apply a heat flux to the SNRs and determine the temperature gradient induced by this flux.[32] A time step of 0.5 fs was taken for integration of the atomic equations of motion. Tersoff potential[33] was used to describe the atomic interactions in the simulation systems. Periodic boundary conditions were used in the x and y directions, and free boundary was used in the z direction (Fig. 1(a)). Heat fluxes were imposed in the system by continuously adding heat to the hot region and removing heat from the cold region. The temperature gradient was then generated along the x direction. The initial configuration was relaxed at room temperature, which was achieved by constant volume and constant temperature (NPT) MD simulation for 4 × 106 time steps (2 ns). After the relaxation, the bond lengths of SNRs produced by Tersoff potentials each stayed at around 2.30 Å, a value close to that reported from ab initio simulation.[34] The equilibrium state was reached by the constant volume and constant temperature (NVT) simulation under the Nose–Hoover thermostat with a coupling constant of 0.05 ps.[35] Next, a heat flux J was imposed in these two regions by exchanging the kinetic energies between the hottest atoms in the hot region and the coldest atoms in the cold region in a microcanonical ensemble. Under the microcanonical ensemble, we ran 1.2 × 107 time steps (6 ns) to allow the system to reach steady state. At each time step, a small amount of heat Δε was added into a hot region and removed from the cold region of the SNR. When the steady state was achieved, the heat flux could be calculated from

where A is the cross-section of heat transfer, defined by width times thickness of the SNR. The SNR thickness was assumed to be 0.42 nm[36] and was the time step. The factor of 2 accounted for the fact that the heat currents propagated in two directions away from the hot region.[37] To calculate the temperature gradient, the sample was divided into N slabs in the x direction. The temperature in each slab was calculated based on the velocities of atoms inside that slab from
where Ni is the number of atoms in slab i; mj and vj are the mass and velocity of atom j in the slab, respectively; kB is the Boltzmann constant. In the microcanonical ensemble, we can obtain the time-averaged temperature profile of the system (Fig. 1(b)), where the linear fitting region is used to calculate the temperature gradient δT/δx. Based on Fourier’s law of heat conduction, the thermal conductivity κ is calculated from the following equation:
To reduce the fluctuation in the simulated thermal conductivities, results were averaged over 10 realizations. The phonon density of state (PDOS) was calculated to gain an insight into the mechanism of the phonon thermal conductivity. The phonon spectrum function Pω was calculated from the Fourier transform of the velocity autocorrelation function,
where vj(0) is the average velocity vector of a particle at initial time, vj(t) is its velocity at time t, and ω is the vibrational wavenumber.

3. Results and discussion
3.1. Thermal conductivity of pristine SNRs

Figure 2 shows the thermal conductivities of SNRs with different lengths and edges. The thermal conductivity of SNR increases with length L increasing when the length of SNR is less than 80 nm. When the SNR length is greater than 80 nm, the thermal conductivity of SNR converges to a stable value. The converged thermal conductivity values are estimated to be 32 W/mK and 35 W/mK for the studied armchair-edged SNRs and zigzag-edged SNRs, respectively, similar to previously reported results.[12,18] The thermal conductivity of SNR shows anisotropy, and the thermal conductivity of armchair-edged SNRs is 10% smaller than that of zigzag-edged SNRs as shown in Fig. 2. Lattice thermal conductivity depends on how far phonons travel between scattering events, that is, how long their mean free path (MFP) is.[38] The MFPs of phonons along the armchair-edged and zigzag-edged directions are estimated to be 1.35L and 1.15L respectively.[28] Thus longer MFP phonons restricted on a nanoscale can be remarkably scattered by the boundary, resulting in a lower thermal conductivity. Our study demonstrates that the lattice thermal conductivity of SNR shows size and edge chirality dependencies; a similar phenomenon is also found in graphene nanoribbons.[39,40]

Fig. 2. (color online) Length-dependent thermal conductivities of SNRs with zigzag and armchair edges. The width of the studied SNRs is about 3 nm.
3.2. Effect of random isotopic doping on thermal conductivity

To reduce the thermal transport,30Si isotope impurity is introduced to increase phonon scattering in SNRs. Isotopic doping can be done either in a random manner or in an ordered manner. Figure 3 displays the calculated thermal conductivity of randomly doped armchair-edged SNRs as a function of30Si doping concentration (ρ). The thermal conductivity of the studied SNR decreases initially to a minimum value and then increases as the doping concentration changes from 0% to 100%, and the minimum value (about 15.5 W/mK) occurs at a doping concentration of around 50%. Particularly, the calculated thermal conductivity of pure30Si SNR (ρ = 100%) is found to be approximately 3.8% smaller than that of pure28Si SNR (ρ = 0%, κ = 18.3 W/mK) because heavier isotope atoms reduce the lattice vibration frequency, thus resulting in lower thermal conductivity. This U-shaped change of thermal conductivity is also found in isotope-doped silicon nanowires[29] and graphene.[41] The isotopic-doping-induced reduction of thermal conductivity is due to a point defect induced by a mass difference in a crystal lattice, which causes phonon scattering; this finding is controlled by the mass and concentration of the individual isotopes that contribute to the disorder.[42]

In addition, we theoretically study κ reduction ratio induced by isotope doping based on a mean-field model[43] on the assumption that the heat transfer rate is proportional to frequency. The mean-field approximation suggests that the heat transfer reduction ξ can be given as

where Mα and Mβ are different point mass, κϕ is the thermal conductivity of isotope substituted nanomaterial, κ0 is the thermal conductivity of the pure nanostructure, ϕ is the molar fraction of the isotope β that replaces atom α in the nanomaterial, and ε is a dimensionless numerical parameter to prevent divergence for the limiting cases. In a completely harmonic mean-field approximation, parameters ε, Mα, and Mβ are determined using the following two boundary conditions: ξ(1)=(mβmα)−1/2, where mα and mβ are the atomic mass of isotopes.

For the studied SNR with buckled structure, anharmonicity considerably contributes to the thermal conductivity. Srilok et al. optimized the above model for anharmonic case.[36] They used necessary corrections to mean field parameters, where κ100, κ50, and κ0, are computed from MD simulations at 100%, 50%, and 0% doping respectively. For the studied SNR, parameters calculated from the above boundary conditions are[43]

The thermal conductivity reduction ratios due to different doping concentrations are calculated by using the mean field model [Eq. (6)]. The results are compared with the MD simulation results and are shown in Fig. 3. The results from the optimized mean field model are in good agreement with the MD calculations.

Fig. 3. (color online) Plot of the thermal conductivity reduction ratio ξ against doping concentration ρ for30Si randomly doped armchair-edged SNR, calculated with MD simulations (blue solid circle) and mean field model (red circle). In this study, thermal conductivities are normalized by the thermal conductivity of pure28Si SNR. The armchair-edged SNR is 30 nm in length and 3 nm in width.

The PDOS of randomly doped armchair-edged SNR is calculated to understand the mechanism of thermal conductivity reduction induced by isotopic doping, and the results are shown in Fig. 4. We examine the phonon spectra for ρ = 0%, 20%, 40%, 60%, and 100%. With increasing ρ from 0% to 100%, the frequencies of the phonon modes show redshifts. When ρ increases from 0% to 40% or decreases from 100% to 60%, the mass disorder increases, thus lowering the intensities of the PDOS peaks. Compared with a pure system (with no doping), the doping isotope also softens the G-band of the phonons. Hence, isotope substitution induces mass disorder in the lattice, resulting in increased phonon scattering at defective sites because of the differences in the characteristic frequency among phonons. Those impurities cause phonon modes to be localized, thus reducing the thermal conductivity.

Fig. 4. (color online) PDOS of randomly doped armchair-edged SNR with different values of doping concentration ρ of30Si.
3.3. Effect of ordered isotopic doping on thermal conductivity

Next, the thermal conductivity of SNR doped by30Si isotopes is simulated in an ordered manner (Fig. 5(a)). The SNR with a length of 30 nm and a width of 30 nm is used as an example. We investigate the influence of superlattice arrangement on thermal conductivity. Figure 5(b) shows the thermal conductivity of superlattice-structured SNR as a function of period length LP for doping concentration ρ = 40%. When LP increases gradually from 3 nm to 10 nm, the κ value of the studied superlattice-structured SNR decreases from 31 W/mK (LP = 3 nm) to 27 W/mK (LP = 7.5 nm) and to 25 W/mK (LP = 10 nm). This outcome can be explained qualitatively from phonon spectrum variations induced by superlattice arrangement. Figure 6 shows the phonon spectra of pure SNR and ordered doped SNRs with LP = 3 nm and with LP = 10 nm. Peaks of phonon modes corresponding to G-band are suppressed gradually when going from pure SNR to the ordered doped SNRs with LP = 3 nm and the ordered doped SNRs with LP = 10 nm, thus lowering the thermal conductivity.

Fig. 5. (color online) (a) Schematic diagram of the superlattice-structured SNR. (b) Thermal conductivity of superlattice-structured SNR as a function of period length LP for doping concentration ρ = 40%. The heat flows perpendicularly to the doping strip region.
Fig. 6. (color online) PDOS spectra for superlattice-structured SNRs with LP = 3 nm (red line) and LP = 10 nm (blue line) compared with that of pristine SNR (black line).
Fig. 7. (color online) Thermal conductivities of random isotope substitution and superlattice isotope substitution in SNR for the same concentration ρ = 40% compared with the case of pure SNRs (ρ = 0% and ρ = 100%).

The thermal transports in randomly doped SNR and superlattice-doped SNR are compared at the same concentration (ρ = 40%) of30Si. The results are shown in Fig. 7. Compared with the pure SNR (with no doping), the random doping and the superlattice doping can reduce the thermal conductivity. Interestingly, however, the ordered doping (i.e., isotope superlattice) leads to a much larger reduction in the thermal conductivity than the random doping for the same doping concentration.

4 Conclusions

In the present study, we investigated the phonon thermal conductivities of SNRs by using the RNEMD method. Phonon thermal conductivity of SNR shows size and edge chirality dependencies. We mainly focused on the effects of isotopic doping on the phonon thermal conductivity of SNR doped in a random manner and in an ordered manner. The isotopic-doping-induced reduction of thermal conductivity, which is due to a point defect induced by a mass difference in a crystal lattice, causes phonon scattering. Our simulation results indicate that doping concentration, doping pattern, and doping arrangement have important effects on the phonon thermal conductivity of the SNR studied. Detailed analyses of the phonon spectra demonstrate that defect-induced suppression of the phonon modes is responsible for the reduction of thermal conductivity. This study guides defect engineering of the thermal transport properties in silicene materials and has crucial implications for silicene-based thermoelectric applications.

Reference
[1] Novoselov K S Geim A K Morozov S V Jiang D Zhang Y Dubonos S V Grigorieva I V Firsov A A 2004 Science 306 666
[2] Liu C C Feng W X Yao Y G 2011 Phys. Rev. Lett. 107 076802
[3] Shirzadi B Yarmohammadi M 2017 Chin. Phys. 26 017203
[4] Chowdhury S Jana D 2016 Rep. Prog. Phys. 79 126501
[5] Ni Z Y Liu Q H Tang K C Zheng J X Zhou J Qin R Gao Z X Yu D P Lu J 2012 Nano Lett. 12 113
[6] Sadeddine S Enriquez H Bendounan A Kumar Das P Vobornik I Kara A Mayne A J Sirotti F Dujardin G Oughaddou H 2017 Sci. Rep. 7 44400
[7] Meng L Wang Y L Zhang L Z Du S X Gao H J 2015 Chin. Phys. 24 086803
[8] Aufray B Kara A Vizzini S Oughaddou H Leandri C Ealet B Le Lay G 2010 Appl. Phys. Lett. 96 183102
[9] Wu K H 2015 Chin. Phys. 24 086802
[10] Pulci O Gori P Marsili M Garbuio V Del Sole R Bechstedt F 2012 Europhys. Lett. 98 37004
[11] Li H P Zhang R Q 2012 Europhys. Lett. 99 36001
[12] Hu M Zhang X L Poulikakos D 2013 Phys. Rev. 87 195417
[13] Ng T Y Yeo J J Liu Z S 2013 Int. J. Mech. Mater. Des. 9 105
[14] Lew Yan Voon L C 2015 Chin. Phys. 24 087309
[15] Drummond N D Zolyomi V Fal’ko V I 2012 Phys. Rev. 85 075423
[16] Wang Z Y Feng T L Ruan X L 2015 J. Appl. Phys. 117 084317
[17] Zhang X L Bao H Hu M 2015 Nanoscale 7 6014
[18] Liu Z Y Wu X F Luo T F 2017 2D Mater. 4 025002
[19] Zhao W Guo Z X Zhang Y Ding J W Zheng X J 2016 Solid State Commun. 227 1
[20] Wirth L J Osborn T H Farajian A A 2016 Appl. Phys. Lett. 109 173102
[21] Kuang Y D Lindsay L Shi S Q Zheng G P 2016 Nanoscale 8 3760
[22] Araujo P T Terrones M Dresselhaus M S 2012 Materials Today 15 98
[23] Yang N Zhang G Li B 2008 Nano Lett. 8 276
[24] Yu B Zebarjadi M Wang H Lukas K Wang H Z Wang D Z Opeil C Dresselhaus M Chen G Ren Z F 2010 Nano Lett. 12 2077
[25] Sevincli H Sevik C Cagin T Cuniberti G 2013 Sci. Rep. 3 1228
[26] Ni Z Zhong H Jiang X Quhe R Luo G Wang Y Y Ye M Yang J Shi J Lu J 2014 Nanoscale 6 7609
[27] Guo Y Zhou S Bai Y Zhao J 2015 J. Supercond. Nov. Magn. 29 717
[28] Liu B Reddy C D Jiang J W Zhu H W Baimova J A Dmitriev S V Zhou K 2014 J. Phys. D: Appl. Phys. 47 165301
[29] Zhang G Zhang Y W 2013 Phys. Status Solidi 7 754
[30] Muller-Plathe F 1997 J. Chem. Phys. 106 6082
[31] Plimpton S 1995 J. Comput. Phys. 117 1
[32] Cao B Y Li Y W 2010 J. Chem. Phys. 133 024106
[33] Tersoff J 1989 Phys. Rev. 39 5566
[34] Lebegue S Eriksson O 2009 Phys. Rev. 79 115409
[35] Nose S 1984 Mol. Phys. 52 255
[36] Srinivasan S Ray U Balasubramanian G 2016 Chem. Phys. Lett. 650 88
[37] Pei Q X Sha Z D Zhang Y W 2011 Carbon 49 4752
[38] Regner K T Sellan D P Su Z H Amon C H McGaughey A J H Malen J A 2013 Nat. Commun. 4 1640
[39] Ye Z Q Cao B Y Yao W J Feng T L Ruan X L 2015 Carbon 93 915
[40] Aksamija Z Knezevic I 2011 Appl. Phys. Lett. 98 141919
[41] Hu J Schiffli S Vallabhaneni A Ruan X Chen Y P 2010 Appl. Phys. Lett. 97 133107
[42] Pei Q X Zhang Y W Sha Z D Shenoy V B 2013 J. Appl. Phys. 114 033526
[43] Balasubramanian G Puri I K Bohm M C Leroy F 2011 Nanoscale 3 3714