Thermal Conductivity and Tensile Response of Phosphorene

Jun 14, 2017 - *(A.W.-K.L.) Telephone/Fax: +65 6790 5296. E-mail: [email protected]., *(K.Z.) Telephone: +65 6790 5499. ...... Buscema , M.; Groenendi...
2 downloads 7 Views 4MB Size
Download PDF
Article pubs.acs.org/JPCC

Thermal Conductivity and Tensile Response of Phosphorene Nanosheets with Vacancy Defects Bo Liu,†,‡ Lichun Bai,‡ Elena A. Korznikova,§ Sergey V. Dmitriev,§,∥ Adrian Wing-Keung Law,*,†,⊥ and Kun Zhou*,†,‡ †

Environmental Process Modelling Centre, Nanyang Environment and Water Research Institute, Nanyang Technological University, 50 Nanyang Avenue, Singapore 639798, Singapore ‡ School of Mechanical and Aerospace Engineering, Nanyang Technological University, 50 Nanyang Avenue, Singapore 639798, Singapore § Institute for Metals Superplasticity Problems, Russian Academy of Sciences, 39 Khalturin Street, 450001 Ufa, Russia ∥ National Research Tomsk State University, Lenin Avenue, 36, 634050 Tomsk, Russia ⊥ School of Civil and Environmental Engineering, Nanyang Technological University, 50 Nanyang Avenue, Singapore 639798, Singapore ABSTRACT: Phosphorene, a newly discovered member of the two-dimensional nanomaterials family, has attracted intensive interest recently. In this paper, molecular dynamic simulations have been conducted to investigate the thermal conductivity and tensile response of phosphorene nanosheets containing vacancy defects. Three vacancy types, the single-vacancy and two types of divacancies, have been considered. The simulation results show that both the thermal conductivity and mechanical strength are highly anisotropic and are largely reduced by the defects. Interestingly, the effects of the defects are also anisotropic and the anisotropy occurs between the two types of divacancies. Along the armchair direction, one type divacancy leads to better thermal and mechanical properties than the other type, while this relation reverses along the zigzag direction. For the thermal conductivity, the effects of strain and temperature have also been studied. As the temperature increases, the thermal conductivity decreases, particularly along the armchair direction. Moreover, the temperature dependence is weakened when the defects are presented. Under the tensile loading, regardless of the presence of defects, the thermal conductivity along the armchair direction increases with the strain monotonically, while that along the zigzag direction shows a nonmonotonic dependence. Overall, these findings provide helpful insights for the understanding of the thermal and mechanical properties of phosphorene, and have significance for its future applications in novel electronic devices. Besides being applied in field-effect transistors, phosphorene is also found to have a figure of merit higher than 1.0 at 300 K20 and can be further enhanced by applying tensile strain.21 The high ZT value of phosphorene makes it a competitive material to be utilized in thermoelectric devices.22−24 Since the thermoelectric efficiency is highly determined by the thermal conductivity, so far several simulation works have been devoted to studying the thermal conductivity of phosphorene. On the basis of the first-principles calculations, Liu and Chang found that the thermal transport in phosphorene is highly anisotropic.25 They reported that the thermal conductivity along the zigzag direction is 48.9 W m−1 K−1, which is almost two times that of 27.8 W m−1 K−1 along the armchair direction. Zhu et al.26 also reported an anisotropic ratio of 4, and the thermal conductivity along the armchair and zigzag

1. INTRODUCTION Two-dimensional (2D) nanomaterials, characterized by a thickness in the order of subnanometer, have attracted tremendous interests over the past decade due to their extraordinary properties as well as their promising applications in advanced nanotechnologies. The existing library of 2D nanomaterials includes graphene,1 silicene,2−5 hexagonal boron nitride (h-BN),6 molybdenum disulfide (MoS2)7−10 and phosphorene.11−13 Among these 2D nanomaterials, phosphorene, a counterpart of black phosphorus in the form of single layer phosphorus atoms bonded with sp3-hybridization, is one of the most recently discovered 2D materials.11 Unlike graphene and MoS2, phosphorene possesses, at the same time, a large and tunable band gap from 0.3 to 1.0 eV14 and a high carrier mobility15 of up to 10 000 cm2 V−1 s−1. Such unique properties might make phosphorene more promising than graphene in nanoelectronics, such as field-effect transistors16−18 and phototransistors.19 © XXXX American Chemical Society

Received: March 28, 2017 Revised: May 16, 2017

A

DOI: 10.1021/acs.jpcc.7b02933 J. Phys. Chem. C XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry C

Figure 1. Front and top views of the atomic structures of pristine phosphorene and phosphorene with vacancy defects. (a) Pristine phosphorene with the P atoms at the top and bottom atomic planes being depicted by the gray and green colors, respectively. (b) Different types of bonds and angles in phosphorene. (c) One P atom at the top atomic plane is removed, leading to the formation of a single-vacancy (SV). (d) Two neighboring P atoms from the same atomic plane are removed, forming the first type of divacancy (DV1). (e) Two neighboring P atoms from different atomic planes are removed, resulting in the second type of divacancy (DV2). The deleted atoms are circled by the red dash lines.

directions are found to be 24.3 and 83.5 W m−1 K−1, respectively. Because of the anisotropy of the thermal conductivity, the thermoelectric efficiency is also found to be anisotropic according to the work by Leonardo et al.27 The anisotropy of the thermal properties of phosphorene has further been confirmed by molecular dynamics (MD) simulations. For example, Xu et al. demonstrated that the thermal anisotropy is mainly due to the anisotropy of phonon group velocity and lifetime.28 Zhang et al. reported that the thermal conductivity and its anisotropy are insensitive to the number of phosphorene layers.29 Furthermore, they found that tensile strain leads to an anomalous increase in the in-plane thermal conductivity of phosphorene, in particular in the armchair direction. However, this finding is opposite to the results reported by Ong et al., who found that the thermal conductivity along the armchair direction decreases when a tensile strain is applied.30 During the fabrication and application of nanomaterials, defects such as vacancies, line defects, and grain boundaries are usually inevitable and phosphorene is not exceptional.31−35 Previous first-principles studies have shown that the defect formation energy in phosphorene is much lower than that in graphene and silicene,36 and the vacancy defect can contribute to greatly increasing the electrical current along the zigzag direction, particularly for the single vacancy.37 In addition, it has also been demonstrated that the single-vacancy can induce the hole-doping effect in phosphorene and result in local magnetic moments,31,38 which makes phosphorene a promising spintronic material.39 For the divacancy and the Stone−Wales defects, they have little effect on the electronic properties of phosphorene, and thus its semiconducting properties can be well preserved with band gap values similar to that of pristine phosphorene.31 So far, no studies have been reported on the effect of vacancy defects on the thermal transport properties of phosphorene. In addition, to take advantages of the unique properties caused by the vacancy defects in practical application, their effect on the

mechanical stability and strength of phosphorene nanosheets should also be understood. Hence, this work is devoted to investigating the effect of vacancy defects on the thermal transport and mechanical properties of phosphorene via MD simulations. The in-plane thermal conductivity and tensile response of single-layer phosphorene nanosheets are studied. The single-vacancy and two types of divacancies are considered. The coupling effect of defect, strain and temperature on the thermal conductivity are investigated. The elastic modulus, fracture strain and strength are calculated to characterize the tensile response. Detailed analyses of the phonon modes and the structural deformation process are conducted to understand the underlying mechanisms.

2. MODELING Figure 1 shows the atomic structures of pristine single-layer phosphorene and those containing single-vacancy and divacancy defects. Because of the intrinsic puckered structure of phosphorene, the P atoms locate at two atomic planes with a distance of ∼1.90 Å. Therefore, there are two types of bonds as well as two types of bond angles in phosphorene. The bonds connecting two P atoms from the same atomic plane are denoted as the intralayer bonds in this study, while those connecting the P atoms from different atomic planes are denoted as the interlayer bonds. Similarly, the bond angles formed by two intralayer bonds are denoted as intralayer angles, while those formed by one intralayer bond and one interlayer bond are denoted as interlayer angles. To create a single-vacancy (SV), one P atom is removed from either the top or the bottom atomic plane (Figure 1b). When creating a divacancy (DV), the two neighboring P atoms to be removed can be selected either from the same atomic plane or different planes.39 Hence, two types of divacancies can be created, denoted as DV1 and DV2, and shown in Figure 1, parts c and d, respectively. It is noted that the formation energy of DV2 is slightly higher than that of DV1 according to the B

DOI: 10.1021/acs.jpcc.7b02933 J. Phys. Chem. C XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry C previous first-principles study.39 However, both types of divacancies are considered here. The simulation is conducted by using the large-scale atomic/ molecular massively parallel simulator (LAMMPS) package.40 Periodic boundary conditions are applied along the in-plane xand y-directions to remove the edge effect, with the freeboundary condition along the out-of-plane z-direction. To model the interactions of P−P atoms, the Stillinger−Weber (SW) potential with a recently developed parametrization set for the phosphorene monolayer is adopted.41 This potential fully describes the buckled structure of phosphorene and produces phonon dispersion curve that is in good agreement with first-principles calculations. In addition, it is capable of providing reasonably accurate information for nonlinear effects of the covalent systems which are significant for MD simulations of the thermal conductivity. This potential has been adopted in several previous studies to calculate the thermal and mechanical properties of phosphorene,29,41−44 and can be expressed as Φ=

∑ V2 + ∑ i DV1 > DV2 for the zigzag direction. To examine whether these orders are valid under different temperatures, additional simulations have been conducted for all the vacancy types with the concentration of 0.20 nm−2 at temperatures ranging from 175 to 300 K. The simulation results show that these orders are not changed by the temperature. As will be discussed later, the dependence of the vacancy effect on the chirality may be related to the differences in the numbers of bonds and bond angles of the nanosheets with different types of vacancies. In parts c and d of Figure 3, the thermal conductivity λ values for different types of vacancies with a concentration of 0.20 nm−2 at various temperatures are normalized with their corresponding values at 175 K to show the dependence of λ on the temperature. Here, the temperature is taken as the temperature of the heat sinks, and the temperature difference between the heat sinks and the heat source is kept at 50 K. For comparison, λ for the pristine phosphorene is also presented. As the temperature increases, λ decreases monotonically. Moreover, λ along the armchair direction is more sensitive to temperature. For instance, when the temperature increases from 175 to 300 K, λ for the pristine phosphorene decreases by 43% along the armchair, but only 33% along the zigzag direction. When the vacancy defects are present, the dependence on temperature weakens. For the nanosheets with the DV2 vacancies, the λ along the zigzag direction decreases by only 10% as the temperature increases from 175 to 300 K. Thermal transport in phosphorene is dominated by phonons. As the phonons transport in nanostructures, they experience intensive scatterings caused by each other, the defects, the interfaces and boundaries. According to the basic kinetic theory of particles, the thermal conductivity of a solid can be approximated by λ = CV vl/3, where CV is the constant volume heat capacity, v is the group velocity of the phonons and can be deviated from the phonon dispersion relation, and l is the length of mean free path (MFP) of the phonon and represents the characteristic length for the scatterings of phonons. The heat capacity CV is calculated as

CV = 3NkB

∫0

∞

⎛ hω ⎞2 ⎛ hω ⎞ g (ω)⎜ ⎟ exp⎜ ⎟ ⎝ kBT ⎠ ⎝ kBT ⎠

⎛ ⎛ hω ⎞ ⎞2 ⎟⎟ dω − /⎜⎜exp⎜ 1 ⎟ ⎝ ⎝ kBT ⎠ ⎠

(5)

where N denotes the number of atoms, kB is the Boltzmann constant, h is the Planck constant, T is the temperature, and g(ω) is the phonon spectrum power at the frequency ω which can be calculated by performing the fast Fourier transform on the velocity autocorrelation function. The MFP length l of a material is very sensitive to its atomic structure, and is related to the phonon relaxation time τ and group velocity as l = vτ. In the approximation of τ, the effects of various scattering mechanisms which limit l are additive in terms of τ−1 = ∑τi−1, where i numbers all the scattering processes including phonon− phonon scattering and phonon scattering by point defects, dislocations, grain boundaries and so on. To understand the effect of vacancy defect on the thermal conductivity of phosphorene by applying the above particle kinetic theory, the phonon dispersion curve and phonon power spectrum are calculated. During the calculation of the dispersion curve, a supercell containing 5 × 5 unit cells is adopted (Figure 4a) and a vacancy defect is created in the center. By adopting such 5 × 5 supercell, the corresponding defect concentration is around 0.275 nm−2 which is close to the highest concentration considered in this study (Figure 3). It should also be noted that by adopting such supercell, the vacancy defect is treated to be periodically distributed inside the materials. This is different from the random distribution pattern considered for the calculation of the thermal conductivity. However, it is believed that the dispersion relation calculated can still be used to understand the phonon behaviors of the defective phosphorene. In Figure 4b, the low-frequency domain of the phonon dispersion relation for the phosphorene with the SV defect (blue color) is plotted in comparison with that for the pristine phosphorene (red color). The size of the half Brillouin zone is folded to π/5a in the ΓX direction and π/5b in the ΓY direction. With the presence of the SV defect, the acoustic phonon branches are slightly flattened as indicated by the black arrows. This result is consistent with the previous study on the phosphorene phonon crystal.43 Similar results have also been E

DOI: 10.1021/acs.jpcc.7b02933 J. Phys. Chem. C XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry C

and armchair directions of the defective phosphorene nanosheets would be no larger than 10 nm. It was previously reported that for the pristine phosphorene, the lengths of phonon mean free path along the zigzag and armchair direction are 39.5 and 18.9 nm, respectively.29 As a result, the reduction of the length of the phonon mean free path caused by the vacancy defects is larger along the zigzag direction than that along the armchair direction. Hence, based on the relation λ = CVvl/3, the higher sensitivity of the thermal conductivity along the zigzag direction to the vacancy defected can be explained. According to Klemens and Pedraza’s work,50 the phonon relaxation time τp due to phonon-point defect scattering can be expressed as τp−1 = ηCω2g(ω) with C being the defect concentration and η being a parameter determined by the defect nature. Hence, as the phonon power spectrum g(ω) and defect concentration increase, the phonon relaxation time τp decreases and thus the thermal conductivity decreases according to the relation λ = CVvl/3. It is worth mentioning that the change of the heat capacity CV would also lead to the reduction of thermal conductivity. By using eq 5, CV for both the pristine and defective phosphorene have been estimated. The results show that the vacancy defects only lead to less than 5% reduction of CV. Therefore, it can be concluded that the reduction of the phonon relaxation time or reduction of the phonon mean free path caused by the phononvacancy scattering is the dominant reason for the reduction of the thermal conductivity. This also explains the weakening of the temperature dependence of the thermal conductivity for the defective phosphorene. As the temperature increases, more phonons, particularly those with higher frequencies, would be activated. The increased phonon population leads to higher phonon−phonon scattering rate and thus lower thermal conductivity. However, in the presence of the vacancy defects, the existence of the phonon-defect scattering weakens the effect of the phonon−phonon scattering enhancement as the total phonon relaxation time is now determined by both scattering processes according to the relation τ−1 = ∑τi−1. Therefore, the dependence of the thermal conductivity on the temperature weakens. The mechanical strain has been frequently adopted to tune the properties of nanomaterials. Previous MD simulation works have demonstrated that the thermal conductivities of graphene and some other nanostructures including carbon nanotube and silicon nanowires, decrease when the tensile strain is

obtained for the defective phosphorene with the DV1 and DV2 vacancies and thus are not presented here. As the slopes of the acoustic branches give the phonon group velocities, the slight flattening of the acoustic branches demonstrates that the phonon group velocities are not significantly affected by the vacancy defects and thus are not the dominant reason for the reduction of the thermal conductivity. Figure 5 shows the phonon power spectrum of the defective phosphorene with the SV vacancy at the defect concentration

Figure 5. Phonon power spectrum (red line) of the defective phosphorene with the SV vacancy at the defect concentration of 0.20 nm−2 and that for the pristine phosphorene (black line).

of ∼0.20 nm−2 in comparison with that for the pristine phosphorene. It is observed that the phonon spectra overlap well in the frequency region above 5 THz and the characteristic peaks at ∼6 and ∼14 THz are also preserved well. However, the peak heights are reduced by the vacancy defects. Moreover, in the region below 5 THz, the power spectrum for the defective phosphorene is enhanced. This indicates that the phonon modes with higher frequencies in the range of 6−15 THz are scattered by the vacancy defects and turn into phonon modes with lower frequencies. The same phenomenon can be observed for the defective phosphorene nanosheets with higher defect concentrations and for the DV1 and DV2 defects. As the phonons transport through the phosphorene nanosheets, they are frequently scattered by the vacancy defects, thus leading to the reduction of the thermal conductivity. For the smallest defect concentration of 0.04 nm−2 considered, the average distance between two vacancy defects is around 10 nm. Hence, the length of the phonon mean free path for both the zigzag

Figure 6. Thermal conductivity of the phosphorene nanosheets as a function of the tensile strain for the three types of vacancy defects at the defect concentrations of 0.08 and 0.20 nm−2: (a) armchair direction; (b) zigzag direction. The thermal conductivity of the pristine phosphorene is also plotted for comparison. All the data points are normalized by the thermal conductivity at the zero-strain state. F

DOI: 10.1021/acs.jpcc.7b02933 J. Phys. Chem. C XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry C

Figure 7. Dependence of the (a, b) fracture strain, (c, d) fracture strength, and (e, f) elastic modulus on the defect concentration for the three vacancy types. The left panel shows the results for the armchair direction while the right panel for the zigzag direction.

applied.51,52 This strain-dependence of thermal conductivity is common in bulk materials, such as silicon53 and diamond.54 Interestingly, it has been reported that for the 2D nanomaterials with intrinsic buckling, including silicene and phosphorene, their thermal conductivity shows an anomalous dependence on the tensile strain and increases with it.29,55,56 Recent MD simulations have shown that the thermal conductivity of pristine phosphorene along the armchair direction can be enhanced by around 15% with a tensile strain of 8%. However, the enhancement along the zigzag direction is much weaker.29 When the vacancy defects are considered to be presented in the phosphorene nanosheets, it is interesting to examine whether the anomalous strain-dependence of their thermal conductivity remains. Figure 6 shows the effect of tensile strain on the thermal conductivity of the pristine and defective

phosphorene nanosheets along the armchair and zigzag directions. The thermal conductivities at various strains are normalized by their corresponding values at the zero-strain state. The defect concentrations of 0.08 and 0.20 nm−2 are considered for each vacancy types. The fracture strains of the defective phosphorene along the armchair and zigzag directions are around 10% and 4%, respectively, as will be shown later. Therefore, for the armchair direction, the strain ranging from 0 to 8% is considered, while for the zigzag direction, the strain ranges from 0 to 3%. It is observed that the thermal conductivity along the armchair direction for the pristine phosphorene nanosheets increases monotonically with the tensile strain. For the thermal conductivity along the zigzag direction, it first slightly increases with the tensile strain and then decreases when the strain G

DOI: 10.1021/acs.jpcc.7b02933 J. Phys. Chem. C XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry C

Figure 8. (a) Evolution of the average bond lengths and (b) evolution of the average bond angles before the fracture of phosphorene nanosheets for the tensile loading along the armchair and zigzag directions. (c) The separate and total potential energy changes per atom caused by the bond length and angle changes. (d) Decomposition of the potential energy change per atom into the intralayer and intralayer components. Here, “IntraBond” stands for the intralayer bond, similar to the others. (e) Snapshot of the fracture of the phosphorene nanosheet along the armchair and (f) along the zigzag direction.

lead to a significant reduction in the phonon scattering rate between the in-plane acoustic and out-of-plane flexural modes, thus increasing their lifetimes and the overall thermal conductivity. Hu et al.55 suggested that the strain induced thermal conductivity enhancement of silicene is related to the phonon group velocity and the unique deformation characteristics of the hexagonal lattice. Therefore, a thorough understanding of the strain-dependence of the thermal conductivity of phosphorene would require the detailed analysis of the interactions among the phonon modes from different dispersion branches, their lifetime and group velocity, as well as the mechanical deformation patterns of the nanosheets. This will be presented in a future study. 3.2. Tensile Response. The results presented above indicate that the presence of the vacancy defects would be efficient in tuning the thermal properties of phosphorene. However, the vacancy defects would also highly degrade the mechanical properties. Hence, to make use of the defective engineering in modifying the thermal properties of phosphorene, its tensile response should also be understood. Figure 7

exceeds 1%. These results are in good agreement with the previous study.29 With the presence of the vacancy defects, such strain-dependence can still be observed and is not significantly affected by the defect concentration. However, instead of increasing linearly with the tensile strain as reported,29 the simulation results show that the thermal conductivity along the armchair direction increases with the strain in a nonlinear manner. It first increases significantly when the strain is smaller than 3% and the increasing trend tends to slow down as the strain further increases, which is similar to the phenomenon observed for silicene.55,56 By analyzing the phonon power spectrum of phosphorene, Zhang et al.29 attributed the enhancement of the thermal conductivity of phosphorene mainly to the stiffening of the phonon modes caused by the tensile strain along the armchair direction. As proposed by Xie et al.,56 for 2D nanomaterials with intrinsic buckling such as silicene, the strain-dependent phonon lifetime is the main reason for the enhancement of the thermal conductivity. This is because, under the tensile strain, the buckling nanosheets tend to be more planar. This would H

DOI: 10.1021/acs.jpcc.7b02933 J. Phys. Chem. C XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry C Table 1. Number of Intralayer and Interlayer Bonds/Angles at the Defect Concentration of 0.20 nm−2 pristine_armchair SV_armchair DV1_armchair DV2_armchair pristine_zigzag SV_zigzag DV1_zigzag DV2_zigzag

intrabond

interbond

intraangle

interangle

no. of atoms

1600 1580 1570 1560 1584 1564 1554 1544

800 790 780 790 792 782 772 782

1600 1570 1560 1540 1584 1554 1544 1524

3200 3140 3100 3120 3168 3108 3068 3088

1600 1590 1580 1580 1584 1574 1564 1564

the intralayer angle increases but the interlayer angle is almost unchanged. To further identify the dominant factor that determines the tensile response of the phosphorene, the separate and total potential energy changes caused by the changes of the bond length and bond angle are plotted in Figure 8c. It is found that the total potential energy change for the tensile loading along the zigzag direction is much faster than that along the armchair direction, which explains the much higher elastic modulus for the zigzag direction. For the armchair direction, the potential energy change caused by the bond elongation is close to that caused by the bond angle expansion when the strain is small. When the strain is larger than 5%, the energy change caused by the former becomes lower than that by the latter and the difference increases with the strain. This indicates the deformation of the bond angles is the main response to the tensile loading along the armchair direction. For the zigzag direction, the bond elongation and bond angle expansion contribute almost equally to the overall potential energy change. To identify which bond or angle type contributes more to the mechanical strength, the total potential energy changes are further decomposed into the four intralayer and interlayer bond/angle components as shown in Figure 8d. It is clearly seen that for the zigzag direction, the interlayer bonds and angles almost do not contribute to the potential energy change, while the intralayer bonds and angles dominate the potential energy change with the former playing a more significant role. For the tensile loading along the armchair direction, when the tensile strain is smaller 3%, the intralayer and interlayer bonds and angles contribute almost equally to the potential energy change. As the strain increases to be larger than 3%, the energy change caused by the interlayer angles becomes much larger than the other three components. Overall, the fracture strain and strength along the armchair direction are mainly determined by the interlayer angles, while those along the zigzag direction are determined by both the intralayer bonds and angles. This observation at the temperature of 300 K is in good agreement with the recent study on the fracture patterns of phosphorene at 1 K.57 Such observation can further be evidenced by the snapshots of the fracture of the phosphorene nanosheets along the armchair and zigzag directions, as shown in Figure 8, parts e and f. For the tensile loading along the armchair direction, the fracture is realized via the breaking of the interlayer bonds and the fracture surface is fairly smooth. However, the fracture for the tensile loading along the zigzag direction, the fracture surface is rough due to the breaking of the intralayer bonds. On the basis of the understanding of the deformation pattern, the effect of the vacancy defects on the mechanical properties can be explained. In Table 1, the numbers of

shows the effects of the defect on the elastic modulus, fracture strain and strength of phosphorene. Here, the elastic modulus is obtained from the Hooke’s law based on the stress−strain curve for the strain below 2%. Parts a, c, and e of Figure 7 on the left panel show the results for the armchair direction, while the rest on the right panel show the results for the zigzag direction. For the pristine phosphorene, i.e., with zero defect concentration, the fracture strain, strength and elastic modulus along the armchair direction are 0.120, 2.16, and 20.84 GPa, respectively. These values for the zigzag direction are 0.056, 4.14, and 98.14 GPa, respectively. These results are in good consistency with previous studies.42,57 It is noted that the fracture strain for the armchair direction is much higher than that for the zigzag direction, but the fracture strength and elastic modulus are much lower, indicating that the mechanical properties of phosphorene are also highly anisotropic, in similarity to its thermal properties. As the defect concentration increases, regardless of the defect types, the fracture strain, strength and elastic modulus decrease monotonically. Interestingly, their sensitivities to the increase of defect concentration also depend on the chirality. For the fracture strain and strength, their values along the armchair direction follow the order of SV > DV2 > DV1. However, for the zigzag direction, this order changes to SV > DV1 > DV2. This phenomenon coincide with that observed for the dependence of the thermal conductivity on the defect types (Figure 3, parts a and b). Such coincidence implies that the vacancy defects may intrinsically affect the thermal conductivity and tensile response in the same manner. For the elastic modulus, its values for the three vacancy types are close to each other when the defect concentration is low. As the concentration increases to be larger than 0.12 nm−2, the elastic modulus for the armchair direction follows the order of SV > DV1 ≈ DV2, while for that along the zigzag direction, this order changes to SV ≈ DV1 > DV2. A common feature that shares among the fracture strain, strength and the elastic modulus is that their differences (if any) for different vacancy types increase as the defect concentration increases. To understand the underlying mechanism of the above dependence of the mechanical response, the deformation pattern of the pristine phosphorene is first analyzed. Toward this purpose, the average bond length and angle evolutions before the fracture of the nanosheets are plotted in Figure 8, parts a and b, respectively. It can be observed that for the tensile loading along the armchair direction, both the intralayer and interlayer bond lengths increase with the strain, the intralayer angle decreases slightly, and the interlayer angle increases apparently. For the loading along the zigzag direction, only the intralayer bond length increases significantly and the interlayer bond length is almost unaffected. At the same time, I

DOI: 10.1021/acs.jpcc.7b02933 J. Phys. Chem. C XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry C

along the zigzag direction is much higher than that along the armchair direction, and both of them decrease monotonically with the defect concentration. Moreover, the thermal conductivity along the zigzag direction is more sensitive to the defects and presents higher reductions with the increase of defect concentration. By applying the basic kinetic theory of particles, and analyzing the phonon dispersion curve and power spectrum, it is found that the presence of defects does not significantly reduce the phonon group velocity and the heat capacity, and the reduction of the thermal conductivity is mainly caused by the phonon-defect scattering process. The strain effect on the thermal conductivity of the defective phosphorene nanosheets is studied. When tensile strain is applied along the armchair direction, the thermal conductivity along this direction shows an abnormal increase, regardless of the presence of the vacancy defects and their type. However, when the tensile strain is applied along the zigzag direction, the thermal conductivity along this direction first increases slightly and then decreases, showing a nonmonotonic dependence on the strain. The effect of temperature on the thermal conductivity is also studied. As the temperature increases, the thermal conductivity decreases monotonically. The thermal conductivity along the zigzag direction is found to be more sensitive to the temperature. Furthermore, with the presence of the vacancy defects, the dependence of the thermal conductivity on the temperature is weakened. One interesting phenomenon observed is that besides the anisotropy of the thermal conductivity of the pristine phosphorene, the effect of the vacancy defects on the thermal conductivity is also anisotropic. With the same defect concentration, the thermal conductivity along the armchair direction for different types of vacancies follows the order of SV > DV2 > DV1. However, this order changes to SV > DV1 > DV2 for the zigzag direction. Similar anisotropy is also found for the effect of the vacancy defects on the fracture strain and strength. The presence of defects largely degrades the mechanical strength of phosphorene. With the same defect concentration, the fracture strength and strain along the armchair direction also follow the order of SV > DV2 > DV1, while those for the zigzag direction follow the order of SV > DV1 > DV2. For the elastic modulus, its value along the armchair direction is highest for the SV vacancy, and for the DV1 and DV2 vacancies, the values are close. In other words, the elastic modulus for the armchair direction for different types of vacancies follows the order of SV > DV1 ≈ DV2. For the zigzag direction, the order shifts to SV ≈ DV1 > DV2. Detailed analysis of the bond lengths and angles evolution during the tensile loading process indicates that the mechanical strength along the armchair direction is mainly determined by the interlayer angles, while that along the zigzag direction is determined by both the intralayer bonds and angles. The specific atomic structures of different types of vacancies lead to different numbers of the bonds and angles, thus leading to the anisotropic dependence of the mechanical strength on the defect.

different types of bonds and angles are summarized for different vacancy types at the defect concentration of 0.20 nm−2. For the armchair direction, the interlayer angle is the dominant factor for the fracture strain and strength. The number of interlayer angles for the SV type of vacancy is the highest while that for the DV1 type is the lowest. Therefore, the fracture strain and strength follow the order of SV > DV2 > DV1 (Figure 7, parts a and c). Since each vacancy leads to the same reduction of the number of the interlayer angles (6, 10, 8 for the SV, DV1, DV2 vacancies, respectively), the fracture strain/strength differences for different vacancy types increase as the defect concentration increases. For the loading along the zigzag direction, both the intralayer bonds and angles play significant roles in determining the fracture strain and strength. The number of intralayer bonds/angles for the SV type is the highest but that for the DV2 type instead of the DV1 type is the lowest. As a result, the fracture strain and strength follow the order of SV > DV1 > DV2 for the zigzag direction (Figures 7, parts b and d). For the elastic modulus, it is determined based on the stress−strain curve when the strain is lower than 2%. In this strain region, the four components, namely, the intralayer/ interlayer bond/angle, contribute equally to the potential energy change and thus play equal roles in determining the elastic modulus for the tensile loading along the armchair direction. The total number of bonds and angles for the SV type is the highest, while those for the DV1 and DV2 types are equal. Hence, the elastic modulus follows the order of SV > DV1 ≈ DV2 along the armchair direction. However, for the zigzag direction, the intralayer components, namely, the intralayer bonds and angles, play more important roles in determining the potential energy change and thus the elastic modulus, even in the strain region of DV2 along the zigzag direction. Considering that the thermal conductivity of the defective phosphorene follows the same orders as the fracture strain/ strength, the dependence of the thermal conductivity on the defect types may also be explained in a similar way. When the heat transport is along the armchair direction, the interlayer bonds are parallel with the heat flow direction. More interlayer bonds mean more heat transport channels. Therefore, the DV1 vacancies lead to the lowest thermal conductivity along the armchair direction since they cause the largest reduction in the interlayer bonds. Similarly, when the heat transport is along the zigzag direction, the interlayer bonds are perpendicular to the heat flow direction. The intralayer bonds are the main channel for thermal transport. Hence, the DV2 vacancies lead to the lowest thermal conductivity along the zigzag direction since they result in the least number of intralayer bonds. However, this explanation is phenomenological, and an understanding of the effect of the vacancy types requires further analysis of the roles of the intralayer and intralayer bonds in determining the thermal transport behaviors of phosphorene.

■

4. CONCLUSIONS The thermal conductivity and tensile response of the phosphorene nanosheets with vacancy defects are investigated systematically via MD simulations. Three types of vacancies, namely, the SV, DV1, and DV2 vacancies are considered. It is found that the thermal conductivity of the phosphorene nanosheets is highly anisotropic. The thermal conductivity

AUTHOR INFORMATION

Corresponding Authors

*(A.W.-K.L.) Telephone/Fax: +65 6790 5296. E-mail: [email protected]. *(K.Z.) Telephone: +65 6790 5499. Fax: +65 6792 4062. Email: [email protected]. J

DOI: 10.1021/acs.jpcc.7b02933 J. Phys. Chem. C XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry C ORCID

(16) Li, L. K.; Yu, Y. J.; Ye, G. J.; Ge, Q. Q.; Ou, X. D.; Wu, H.; Feng, D. L.; Chen, X. H.; Zhang, Y. B. Black Phosphorus Field-Effect Transistors. Nat. Nanotechnol. 2014, 9, 372−7. (17) Koenig, S. P.; Doganov, R. A.; Schmidt, H.; Castro Neto, A. H.; Ozyilmaz, B. Electric Field Effect in Ultrathin Black Phosphorus. Appl. Phys. Lett. 2014, 104, 103106. (18) Buscema, M.; Groenendijk, D. J.; Blanter, S. I.; Steele, G. A.; van der Zant, H. S. J.; Castellanos-Gomez, A. Fast and Broadband Photoresponse of Few-Layer Black Phosphorus Field-Effect Transistors. Nano Lett. 2014, 14, 3347−52. (19) Deng, Y. X.; Luo, Z.; Conrad, N. J.; Liu, H.; Gong, Y. J.; Najmaei, S.; Ajayan, P. M.; Lou, J.; Xu, X. F.; Ye, P. D. Black Phosphorus-Monolayer MoS2 van der Waals Heterojunction p-n Diode. ACS Nano 2014, 8, 8292−9. (20) Fei, R. X.; Faghaninia, A.; Soklaski, R.; Yan, J. A.; Lo, C.; Yang, L. Enhanced Thermoelectric Efficiency via Orthogonal Electrical and Thermal Conductances in Phosphorene. Nano Lett. 2014, 14, 6393−9. (21) Lv, H. Y.; Lu, W. J.; Shao, D. F.; Sun, Y. P. Enhanced Thermoelectric Performance of Phosphorene by Strain-induced Band Convergence. Phys. Rev. B: Condens. Matter Mater. Phys. 2014, 90, 085433. (22) Zhang, J.; Liu, H. J.; Cheng, L.; Wei, J.; Liang, J. H.; Fan, D. D.; Shi, J.; Tang, X. F.; Zhang, Q. J. Phosphorene Nanoribbon as a Promising Candidate for Thermoelectric Applications. Sci. Rep. 2015, 4, 6452. (23) Sevik, C.; Sevincli, H. Promising Thermoelectric Properties of Phosphorenes. Nanotechnology 2016, 27, 355705. (24) Zhou, H. B.; Cai, Y. Q.; Zhang, G.; Zhang, Y. W. Thermoelectric Properties of Phosphorene at the Nanoscale. J. Mater. Res. 2016, 31, 3179−86. (25) Liu, T. H.; Chang, C. C. Anisotropic Thermal Transport in Phosphorene: Effects of Crystal Orientation. Nanoscale 2015, 7, 10648−54. (26) Zhu, L. Y.; Zhang, G.; Li, B. W. Coexistence of Size-Dependent and Size-Independent Thermal Conductivities in Phosphorene. Phys. Rev. B: Condens. Matter Mater. Phys. 2014, 90, 214302. (27) Medrano Sandonas, L.; Teich, D.; Gutierrez, R.; Lorenz, T.; Pecchia, A.; Seifert, G.; Cuniberti, G. Anisotropic Thermoelectric Response in Two-Dimensional Puckered Structures. J. Phys. Chem. C 2016, 120, 18841−9. (28) Xu, W.; Zhu, L. Y.; Cai, Y. Q.; Zhang, G.; Li, B. W. Direction Dependent Thermal Conductivity of Monolayer Phosphorene: Parameterization of Stillinger-Weber Potential and Molecular Dynamics Study. J. Appl. Phys. 2015, 117, 214308. (29) Zhang, Y. Y.; Pei, Q. X.; Jiang, J. W.; Wei, N.; Zhang, Y. W. Thermal Conductivities of Single- and Multi-Layer Phosphorene: a Molecular Dynamics Study. Nanoscale 2016, 8, 483−91. (30) Ong, Z. Y.; Cai, Y. Q.; Zhang, G.; Zhang, Y. W. Strong Thermal Transport Anisotropy and Strain Modulation in Single-Layer Phosphorene. J. Phys. Chem. C 2014, 118, 25272−7. (31) Hu, W.; Yang, J. L. Defects in Phosphorene. J. Phys. Chem. C 2015, 119, 20474−80. (32) Guo, Y. Y.; Qiao, C.; Wang, A. H.; Zhang, J. P.; Wang, S. Y.; Su, W. S.; Jia, Y. The Fracture Behaviors of Monolayer Phosphorene With Grain Boundaries under Tension: a Molecular Dynamics Study. Phys. Chem. Chem. Phys. 2016, 18, 20562−70. (33) Guo, Y.; Zhou, S.; Zhang, J. F.; Bai, Y. Z.; Zhao, J. J. Atomic Structures and Electronic Properties of Phosphorene Grain Boundaries. 2D. 2D Mater. 2016, 3, 025008. (34) Jang, W.; Kang, K.; Soon, A. Remarkably Low-Energy oneDimensional Fault Line Defects in Single-Layered Phosphorene. Nanoscale 2015, 7, 19073−9. (35) Cai, Y. Q.; Ke, Q. Q.; Zhang, G.; Yakobson, B. I.; Zhang, Y. W. Highly Itinerant Atomic Vacancies in Phosphorene. J. Am. Chem. Soc. 2016, 138, 10199−206. (36) Zhang, R. Q.; Wu, X. J.; Yang, J. L. Blockage of Ultrafast and Directional Diffusion of Li Atoms on Phosphorene with Intrinsic Defects. Nanoscale 2016, 8, 4001−6.

Kun Zhou: 0000-0002-2152-8774 Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS

■

REFERENCES

The authors acknowledge the financial support from the Ministry of Education, Singapore (Academic Research Fund TIER 1 - RG128/14) and the Nanyang Environment and Water Research Institute (Core Fund), Nanyang Technological University, Singapore. S.V.D. acknowledges financial support from the Russian Science Foundation, Grant No. 14-13-00982. E.A.K. thanks the financial support from the Russian Science Foundation, Grant No. 16-12-10175.

(1) Novoselov, K. S.; Geim, A. K.; Morozov, S. V.; Jiang, D.; Zhang, Y.; Dubonos, S. V.; Grigorieva, I. V.; Firsov, A. A. Electric Field Effect in Atomically Thin Carbon Films. Science 2004, 306, 666−669. (2) Aufray, B.; Kara, A.; Vizzini, S.; Oughaddou, H.; Leandri, C.; Ealet, B.; Le Lay, G. Graphene-like Silicon Nanoribbons on Ag(110): a Possible Formation of Silicene. Appl. Phys. Lett. 2010, 96, 183102. (3) De Padova, P.; Kubo, O.; Olivieri, B.; Quaresima, C.; Nakayama, T.; Aono, M.; Le Lay, G. Multilayer Silicene Nanoribbons. Nano Lett. 2012, 12, 5500−5503. (4) Liu, B.; Baimova, J. A.; Reddy, C. D.; Dmitriev, S. V.; Law, W. K.; Feng, X. Q.; Zhou, K. Interface Thermal Conductance and Rectification in Hybrid Graphene/Silicene Monolayer. Carbon 2014, 79, 236−44. (5) Liu, B.; Baimova, J. A.; Reddy, C. D.; Law, A. W.-K.; Dmitriev, S. V.; Wu, H.; Zhou, K. Interfacial Thermal Conductance of a Silicene/ Graphene Bilayer Heterostructure and the Effect of Hydrogenation. ACS Appl. Mater. Interfaces 2014, 6, 18180−8. (6) Song, L.; Ci, L. J.; Lu, H.; Sorokin, P. B.; Jin, C. H.; Ni, J.; Kvashnin, A. G.; Kvashnin, D. G.; Lou, J.; Yakobson, B. I.; et al. Large Scale Growth and Characterization of Atomic Hexagonal Boron Nitride Layers. Nano Lett. 2010, 10, 3209−15. (7) Radisavljevic, B.; Radenovic, A.; Brivio, J.; Giacometti, V.; Kis, A. Single-layer MoS2 Transistors. Nat. Nanotechnol. 2011, 6, 147−50. (8) Liu, K.-K.; Zhang, W.; Lee, Y.-H.; Lin, Y.-C.; Chang, M.-T.; Su, C.; Chang, C.-S.; Li, H.; Shi, Y.; Zhang, H.; et al. Growth of Large-Area and Highly Crystalline MoS2 Thin Layers on Insulating Substrates. Nano Lett. 2012, 12, 1538−44. (9) Liu, B.; Meng, F. M.; Reddy, C. D.; Baimova, J. A.; Srikanth, N.; Dmitriev, S. V.; Zhou, K. Thermal Transport in a Graphene-MoS2 Bilayer Heterostructure: a Molecular Dynamics Study. RSC Adv. 2015, 5, 29193−200. (10) Jiang, J.-W.; Park, H. S.; Rabczuk, T. Molecular Dynamics Simulations of Single-Layer Molybdenum Disulphide (MoS2): Stillinger-Weber Parametrization, Mechanical Properties, and Thermal Conductivity. J. Appl. Phys. 2013, 114, 064307. (11) Reich, E. S. Phosphorene Excites Materials Scientists. Nature 2014, 506, 19. (12) Liu, H.; Neal, A. T.; Zhu, Z.; Luo, Z.; Xu, X. F.; Tomanek, D.; Ye, P. D. Phosphorene: An Unexplored 2D Semiconductor with a High Hole Mobility. ACS Nano 2014, 8, 4033−41. (13) Kistanov, A. A.; Cai, Y. Q.; Zhou, K.; Dmitriev, S. V.; Zhang, Y. W. Large Electronic Anisotropy and Enhanced Chemical Activity of Highly Rippled Phosphorene. J. Phys. Chem. C 2016, 120, 6876−84. (14) Das, S.; Zhang, W.; Demarteau, M.; Hoffmann, A.; Dubey, M.; Roelofs, A. Tunable Transport Gap in Phosphorene. Nano Lett. 2014, 14, 5733−9. (15) Qiao, J. S.; Kong, X. H.; Hu, Z. X.; Yang, F.; Ji, W. Highmobility Transport Anisotropy and Linear Dichroism in Few-Layer Black Phosphorus. Nat. Commun. 2014, 5, 4475. K

DOI: 10.1021/acs.jpcc.7b02933 J. Phys. Chem. C XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry C (37) Farooq, M. U.; Hashmi, A.; Hong, J. S. Anisotropic Bias Dependent Transport Property of Defective Phosphorene Layer. Sci. Rep. 2015, 5, 12482. (38) Srivastava, P.; Hembram, K.; Mizuseki, H.; Lee, K. R.; Han, S. S.; Kim, S. Tuning the Electronic and Magnetic Properties of Phosphorene by Vacancies and Adatoms. J. Phys. Chem. C 2015, 119, 6530−8. (39) Babar, R.; Kabir, M. Transition Metal and Vacancy Defect Complexes in Phosphorene: A Spintronic Perspective. J. Phys. Chem. C 2016, 120, 14991−5000. (40) Plimpton, S. Fast Parallel Algorithms for Short-Range Molecular-Dynamics. J. Comput. Phys. 1995, 117, 1−19. (41) Jiang, J. W. Parametrization of Stillinger-Weber Potential Based on Valence Force Field Model: Application to Single-Layer MoS2 and Black Phosphorus. Nanotechnology 2015, 26, 315706. (42) Sha, Z. D.; Pei, Q. X.; Zhang, Y. Y.; Zhang, Y. W. Atomic Vacancies Significantly Degrade the Mechanical Properties of Phosphorene. Nanotechnology 2016, 27, 315704. (43) Xu, W.; Zhang, G. Remarkable Reduction of Thermal Conductivity in Phosphorene Phononic Crystal. J. Phys.: Condens. Matter 2016, 28, 175401. (44) Hong, Y.; Zhang, J. C.; Huang, X. P.; Zeng, X. C. Thermal Conductivity of a Two-Dimensional Phosphorene Sheet: a Comparative Study with Graphene. Nanoscale 2015, 7, 18716−24. (45) Hao, F.; Chen, X. First-principles Study of the Defected Phosphorene Under Tensile Strain. J. Appl. Phys. 2016, 120, 165104. (46) Liu, B.; Reddy, C. D.; Jiang, J.; Baimova, J. A.; Dmitriev, S. V.; Nazarov, A. A.; Zhou, K. Morphology and in-plane Thermal Conductivity of Hybrid Graphene Sheets. Appl. Phys. Lett. 2012, 101, 211909. (47) Liu, B.; Reddy, C. D.; Jiang, J. W.; Zhu, H. W.; Baimova, J. A.; Dmitriev, S. V.; Zhou, K. Thermal conductivity of silicene nanosheets and the effect of isotopic doping. J. Phys. D: Appl. Phys. 2014, 47, 165301. (48) Sha, Z. D.; Pei, Q. X.; Ding, Z. W.; Jiang, J. W.; Zhang, Y. W. Mechanical Properties and Fracture Behavior of Single-Layer Phosphorene at Finite Temperatures. J. Phys. D: Appl. Phys. 2015, 48, 395303. (49) Yang, Z. Y.; Zhao, J. H.; Wei, N. Temperature-Dependent Mechanical Properties of Monolayer Black Phosphorus by Molecular Dynamics Simulations. Appl. Phys. Lett. 2015, 107, 023107. (50) Klemens, P. G.; Pedraza, D. F. Thermal-Conductivity of Graphite in the Basal-Plane. Carbon 1994, 32, 735−41. (51) Guo, Z. X.; Zhang, D.; Gong, X. G. Thermal conductivity of graphene nanoribbons. Appl. Phys. Lett. 2009, 95, 163103. (52) Li, X. B.; Maute, K.; Dunn, M. L.; Yang, R. G. Strain Effects on the Thermal Conductivity of Nanostructures. Phys. Rev. B: Condens. Matter Mater. Phys. 2010, 81, 245318. (53) Parrish, K. D.; Jain, A.; Larkin, J. M.; Saidi, W. A.; McGaughey, A. J. H. Origins of Thermal Conductivity Changes in Strained Crystals. Phys. Rev. B: Condens. Matter Mater. Phys. 2014, 90, 235201. (54) Broido, D. A.; Lindsay, L.; Ward, A. Thermal Conductivity of Diamond under Extreme Pressure: A First-Principles Study. Phys. Rev. B: Condens. Matter Mater. Phys. 2012, 86, 115203. (55) Hu, M.; Zhang, X. L.; Poulikakos, D. Anomalous Thermal Response of Silicene to Uniaxial Stretching. Phys. Rev. B: Condens. Matter Mater. Phys. 2013, 87, 195417. (56) Xie, H.; Ouyang, T.; Germaneau, E.; Qin, G. Z.; Hu, M.; Bao, H. Large Tunability of Lattice Thermal Conductivity of Monolayer Silicene via Mechanical Strain. Phys. Rev. B: Condens. Matter Mater. Phys. 2016, 93, 075404. (57) Liu, N.; Hong, J. W.; Pidaparti, R.; Wang, X. Q. Fracture Patterns and the Energy Release Rate of Phosphorene. Nanoscale 2016, 8, 5728−36.

L

DOI: 10.1021/acs.jpcc.7b02933 J. Phys. Chem. C XXXX, XXX, XXX−XXX