Remarkable Role of Grain Boundaries in the Thermal Transport Properties of Phosphorene.

In this work, we study the effect of grain boundary (GB) on the thermal transport of phosphorene by using molecular dynamics simulations. By exploring a total of 19 GBs with different GB defect types and densities, we find that there is a relatively high Kapitza thermal boundary resistance at these boundaries. By analyzing the spatial distributions of the heat flux, we find that this high thermal boundary resistance can be attributed to the strong phonon-boundary scattering at the GBs. With the same type of defect, the thermal boundary resistance is found to increase with the increase of the defect density along the GBs, which can be attributed to the nonuniform distribution of stress and lattice distortion. Finally, we investigate the anisotropy in the thermal conductivity of phosphorene with GBs and reveal a strikingly high anisotropy ratio of thermal conductivities, which is found to arise from the different influences of boundaries on the thermal transport along the zigzag and armchair directions. Our results highlight the importance of GBs in the transport behavior of phosphorene and the need to include their effects in the thermal management of phosphorene-based electronic devices.

Introduction

Strong anisotropy in many of the physical properties of phosphorene (a monolayer of black phosphorus) is one of the most intriguing features because of its unique puckered structure. For example, its band structure exhibits a highly anisotropic dispersion in the vicinity of the gap, where the electrons behave like Dirac fermions along one direction, while they behave like classical ones along the other direction.[1,2] In addition, previous experimental angle-resolved transport studies of few-layer phosphorene also showed a significantly high anisotropy in hole mobility[3−6] as well as in-plane optical conductivity at room temperature.[7] Furthermore, a strong anisotropy in mechanical properties was also observed. It was shown that the elastic constants and the ultimate strain along the zigzag direction are much larger than those along the armchair direction.[8] Recently, experimental[9,10] and theoretical[11,12] studies revealed that phosphorene has a remarkable anisotropy in thermal conductivity. More specifically, the thermal conductivity along the zigzag, armchair, and through-plane directions was found to be 83, 28, and 6.5 W m–1 K–1, respectively, at room temperature.[9,13] These highly anisotropic physical properties have also enabled the exploration of the potential applications of phosphorene. For example, phonons and electrons motion preferably along different directions in phosphorene, which indicates that thermal and electronic transport may be decoupled in phosphorene and potentially opens up new routes for improving the efficiency of thermoelectric converting[14−16] and thermal management in phosphorene-based electrical devices.[17,18]

Previous studies on defects in two-dimensional (2D) materials, such as the vacancy and grain boundary (GB), have shown that they can significantly influence their structural and physical properties.[19−23] In particular, it was shown that these defects in phosphorene are energetically stable and of relatively high concentrations because of their lower formation energies.[24,25] As a result, significant efforts have been made to study the effects of these defects on many of its structural and physical properties. These studies have shown that defects, such as single-vacancy, double-vacancy, and Stone–Wales defects, and GBs affect the physical properties of phosphorene very differently. For example, double-vacancy and Stone–Wales defects have less effects on the band gap of phosphorene. However, single vacancies were shown to influence the electronic properties of phosphorene greatly, and the band gap vanishes in phosphorene containing single vacancies.[25,26] It was also shown that a single-vacancy defect behaves like a p-type impurity, which can greatly increase the electrical conductivity along the zigzag direction.[26] As for GBs in phosphorene, a thorough analysis on their possible structures and energetic stabilities was performed based on density functional theory (DFT) calculations.[24] It was found that GBs do not remarkably affect the electronic properties of phosphorene. For example, its band gap is well preserved, and its electron mobility is only moderately reduced.[24] In terms of mechanical properties, the presence of GBs results in a decrease in fracture strain and Young’s modulus, whereas it has less influence on fracture strength.[27]

We note that unlike the extensive studies that have been performed on the effects of GBs on the electronic and mechanical properties of phosphorene, no systematic studies have been reported on the effects of GBs on the thermal transport properties of phosphorene. Hence, it is of great interest to understand the effects of the GB defect type and density on the thermal conductance and the anisotropy in the thermal conductance of phosphorene. In this work, we would like to investigate the effects of the GB defect type and density on the in-plane thermal transport behavior in GB-contained phosphorene by using molecular dynamics (MD) simulation.

Computational Methods

It is known that GBs in phosphorene are composed of different defected rings (e.g., octagon, pentagon, heptagon, and hexagon rings). Based on the defect type, arrangement, and density, GBs can make various misorientation angles θ between two phosphorene domains. Among the various potential GBs, 19 possible GB structures have been studied using DFT calculations. As shown in Figure S1, these GBs can be categorized into four types.[24] The type I GBs, denoted as a-(5|7) GBs, are composed of heptagon and pentagon rings (5|7 pairs), which orientate to the left direction and are asymmetric about the interface, as shown in N1–N5 (from right to left) of Figure S1. Figure shows an a-(5|7) GB with a misorientation of 20.16°. As shown in N1–N5 of Figure S1, the misorientation angles of a-(5|7) GBs change from 4.39, 9.55, 14.06, 20.16 to 37.42° as the linear defect density increases. As shown in Table S1 and Figure , a smaller dH indicates a higher defect density along the GB (y) direction. The type II GBs are denoted as s-(5|7) GBs, which include 5|7 pairs and are symmetric about the interface, as shown in N6–N10 of Figure S1. The type III GBs are denoted as (5|8|7) GBs, which comprise octagon, heptagon, and pentagon pairs (5|8|7) and are symmetric or asymmetric about the interface, as shown in N11–N15 of Figure S1. In addition to the above GB structures, all other GBs are categorized into type IV GBs, which include 5|8|8|7 pairs, 5|6|7 pairs, or 5|8|7–5|7 pairs, as shown in N16–N19 of Figure S1. The structural parameters of these 19 GBs are given in Table S1. The details of the structural parameters and electronic properties of these GBs can be found in ref (24).

Figure 1

Structure of an a-(5|7) GB with a misorientation angle of θ = 20.16°.

Here, we adopted MD simulation to study the thermal transport across the GBs of phosphorene using the LAMMPS code.[28] An obvious advantage of MD simulation is that it does not depend upon any thermodynamic limit assumption; thus, it is suitable for modeling nanostructures with the real geometric shape. It has been widely applied to investigate the thermal transport properties of different 2D nanostructures, such as graphene,[29] MoS2,[30,31] and phosphorene.[32]

In the present MD simulations, the P–P covalent interactions are described by the Stillinger–Weber (SW) interatomic potential,[33] which is well-established to reproduce the characteristics of phosphorene predicted by first-principles calculations, including the crystal structure, phonon dispersion, and cohesive energy.[34] During the simulations, the atomic motion is calculated from the integration of Newton’s equations with the velocity Verlet algorithm, and the time step is set as 0.25 fs. First, the system was relaxed under the isothermal–isobaric (NPT) ensemble and then equilibrated under the canonical NVT (constant atom number, volume, and temperature) ensemble with a temperature of 300 K for 200 ps using a Nosé–Hoover thermostat.[35] Next, the system was further relaxed under the NVE (constant atom number and volume, and no thermostat) ensemble for additional 200 ps. In this period, we monitored the temperature and total energy (potential energy and kinetic energy) of the system. It was found that the temperature kept constant with small fluctuations around 300 K, and the total energy was conserved, which indicated that the system has reached the equilibrium state. It should be noted that during MD simulations, all the GB structures exhibit robust stability without bond breaking and reconstruction, which indicates that the SW potential is suitable for studying the properties of different GB structures in this work.

To study the thermal properties of the defected phosphorene, the total heat flux J going through the system was calculated by[36]where V is the volume of the simulated system, ε and v denote the energy and velocity of atom i, respectively, is the distance between atom i and atom j, and F and F are the multibody force vectors between atoms i, j, and k. In this work, according to the interlayer distance in bulk phosphorus, the thickness of the phosphorene sheet is set as 5.29 Å.[37]

Results and Discussion

Kapitza Thermal Boundary Resistance

GBs in phosphorene are in the form of one-dimensional interfaces, which cause the change in the heat conductivity in phosphorene. The influence of GBs on thermal transport can be quantified by the Kapitza thermal boundary resistance, which is defined aswhere ΔT and J are the temperature jump and heat flux at the GB, respectively.

Here, we use the nonequilibrium MD (NEMD) method to obtain the Kapitza thermal boundary resistance. The schematic of the atomistic model used in the NEMD simulation is shown in Figure . For the system containing two N5 GBs with opposite orientations, LL = LR = 361.4 Å, LM = 464.6 Å, and the width is 108.6 Å. To eliminate the edge effects on the thermal transport and mimic an infinite width of the system, the periodic boundary condition is applied along the transverse y-direction.

Figure 2

Schematic of GB-contained phosphorene with hot and cold reservoirs applied at the two ends.

To establish a temperature gradient along the longitudinal x-direction and make the heat energy transfer across the GBs, we placed two ends of phosphorene into cold and hot Nosé–Hoover reservoirs, as shown in Figure . The temperatures of the cold and hot reservoirs are TC = 290 K and TH = 310 K, respectively. The simulation was then carried out for 1.0 ns, which is long enough for the system to reach the nonequilibrium steady state, where a constant local heat flux and temperature distribution are well established along the longitudinal direction. After the system reached the nonequilibrium steady state, a time-averaging of the heat flux J and temperature distribution was carried out for 30 ns. The temperature gradient can be obtained according to the slope of the temperature distribution. To avoid the thermal boundary resistance effects, we calculated the slope by fitting the middle region of the temperature profile with a linear function. The averaged results of ΔT and J were then used to calculate the Kapitza thermal boundary resistance R.

A typical temperature distribution along the heat transfer direction in phosphorene containing two N5 GBs at 300 K is shown in Figure a. There is an obvious temperature jump at the GBs, indicating the existence of thermal boundary resistance. The values of temperature jumps ΔT at a-(5|7) GBs are given in Figure b. Based on eqs and 2, the thermal boundary resistances R of a-(5|7) GBs are obtained and shown in Figure c. Because of the heavy computational cost for MD simulations, here, only a-(5|7) GBs and typical N5, N10, N15, and N19, which have the highest defect density along GBs, are calculated for the thermal boundary resistances. It is found that with the increase of the 5|7 defect density along GBs, R increases from 2.39 × 10–10 to 9.93 × 10–10 m2 K W–1, which is 1 order higher than those (from 0.22 × 10–10 to 0.67 × 10–10 m2 K W–1) of graphene with a similar 5|7 defect density along GBs.[38] With the same 5|7 defect density as in the N5 GBs, the R of the N5 GBs is about 4.5 times lower than that of the MoS2–graphene interface (RMoS = 45 × 10–10 m2 K W–1)[39] and that of the phosphorene–graphene interface (RP–G = 42 × 10–10 m2 K W–1).[40] The thermal boundary resistances at the N10, N15, and N19 GBs have also been calculated and are shown in Figure d. It is found that the R of phosphorene GBs can reach as high as 13.81 × 10–10 m2 K W–1. Compared with silicene with the same symmetric 5|7 GBs, the thermal boundary resistance of the N10 GBs is 2.8 times that of its silicone counterpart, which is about 5 × 10–10 m2 K W–1.[41] Furthermore, it is found that the R of GBs in phosphorene is of the same order of magnitude as that of GBs in polycrystalline MoS2, which is in the range of 2.83 × 10–10 to 15.6 × 10–10 m2 K W–1.[42] Our results indicate that GBs play an important role in the thermal transport and thus should be taken into account for the thermal management in 2D material-based devices.

Figure 3

(a) Temperature profile along the heat flux direction in phosphorene with N5 GBs; (b) temperature jump at the a-(5|7) GBs; (c) thermal resistance at a-(5|7) GBs; and (d) thermal resistance of N5, N10, N15, and N19 GBs.

To clarify the mechanism of GBs in the thermal transport of phosphorene, we computed the spatial distribution of the atomic heat flux for each P atom under the nonequilibrium steady state, and the result is plotted in Figure . Here, the heat flux was calculated by using eq and averaged over 30 ns in the steady state. The vector arrows show the magnitude and direction of the heat flux, which vividly reflect the evolution of the heat flux path and phonon scattering around the GB regions. When phonons travel across the GBs, they are scattered at the boundary, and this changes their propagating directions, which results in the increase of thermal resistance along the x-direction. Furthermore, it is found that with the increase of the defect density along the GBs from N1 to N5, the phonon scattering at the GB region is enhanced. Because N5 has the highest defect density among the a-(5|7) GBs studied here, the most extensive phonon scatterings occur at this boundary, resulting in the highest boundary thermal resistance among the a-(5|7) GBs. Figure shows the phonon density of states of two sides of the N5 GB. It is clear that there exist mismatches between the two spectra of left and right sides of the GB, indicating the strong phonon scattering at the GB.

Figure 4

Spatial distributions of the atomic heat flux denoted by vector arrows on each P atom of phosphorene with (a) N1, (b) N3, and (c) N5 GBs.

Figure 5

Phonon density of states at two sides of the N5 GB.

It should also be noted that in this simulation, the structure size of the N5 GB is LL = LR = 361.4 Å, LM = 464.6 Å, and the width is 108.6 Å. To confirm that the simulated length is longer enough to include the contribution of long-wavelength phonons to the boundary thermal conductance, two more N5 GB structures with a longer length have been simulated. The structure sizes and results are given in Table S2 in the Supporting Information. It is found that the difference of the boundary thermal resistance between these three structures is less than 5%, which indicates that the structure size used in our simulations is suitable for calculating the boundary thermal resistance.

Thermal Conductivity Anisotropy in Phosphorene with GBs

Strong anisotropy in the thermal property is one of the most intriguing features of pristine phosphorene. Taking advantage of this property, new strategies for enhancing the thermoelectric efficiency have been proposed.[14−16] Here, we investigate the anisotropy ratio of thermal conductivities for all the 19 GBs as shown in Figure S1 with the equilibrium MD (EMD) method. Structural parameters of these 19 phosphorene GBs are shown in Table S1. It should be noted that Figure S1 and Table S1 give the unit cell structure of each phosphorene GB. To avoid the artificial size effect in EMD calculations, based on these unit cells, large supercell structures (∼120 Å × 120 Å) are used for the EMD simulations, and the periodic boundary condition is applied in the in-plane (x and y) directions.

The thermal conductivities of phosphorene (κ and κ in x and y directions) are calculated with the Green–Kubo formula, which was derived from linear response theory and the fluctuation–dissipation theorem. In the Green–Kubo formula, the thermal conductivity is related to the heat flux autocorrelation function[36] aswhere V denotes the system volume, kB is the Boltzmann constant, and ⟨Ji(0)Ji(t)⟩is the average of the heat flux autocorrelation function. In calculating the thermal conductivity, the periodic boundary condition is applied in the in-plane (x and y) directions. In this work, the total simulation time is 25 ns, the correlation length is 50,000 time steps, and the sample interval is 20.

The calculation results shown in Figure reveal that the influence of GBs is much stronger on κ than κ, which is attributed to the combined effects of the defect density, misorientation angle, and GB density. For the type I GBs at room temperature, with the increase of the defect density (from 0.43 to 1.84 nm–1) and misorientation angle (from 4.39 to 37.42°), κ decreases from 37.48 to 14.54 W m–1 K–1. This is because the increased defect density and misorientation angle will enhance the phonon scattering when phonons transfer across the GBs. In addition, the zigzag direction is a preferred direction for heat transfer compared to the armchair direction in phosphorene. Thus, with the increase of the misorientation angle, the contribution of heat transfer along the zigzag direction to κ decreases.

Figure 6

Thermal conductivities in the x-direction, κ (a), in the y-direction, κ (b), and the anisotropy ratio (c) in phosphorene with different types of GBs. The experimental values of the thermal conductivity of pristine phosphorene in x and y directions from ref (9) are denoted by red dashed lines.

For all GBs with different types, with the increase of the misorientation angle, the thermal conductivity along the y-direction, κ, increases. However, with the increase of the misorientation angle, the defect density along the y-direction also increases, which leads to stronger phonon scatterings at the GBs, and suppresses the heat transfer along the y-direction. These two factors interplay with each other, thus resulting in the slight change of κ among different GBs.

It is well-known that the stress distribution in the region of GBs plays an important role in the phonon transmission and scattering.[29,40] The defects may disrupt regular atomic/bond structures, which induces localized stress around the defects and causes additional phonon scattering. Figure shows the stress fields along the GBs. It is seen clearly that in N1 GBs, the stress fields σ along the x-direction are strongly localized at the core area of defects and decay rapidly with the distance away from the defects. However, in N3 GBs, the stress fields of the defects are largely overlapped, which causes all the atoms/bonds along the GBs to have a nonuniform distribution of compressive and tensile stresses. The large-area distorted lattice structures along the GBs create strong phonon scattering, resulting in a lower thermal conductivity in the x-direction in N3 GBs than in N1 GBs. However, for the stress along the y-direction, in both N1 and N3 GBs, all the stress fields are only localized at the central region of defects. These phenomena are in accord with the values of the thermal conductivity of N1 and N3 GBs, as shown in Figure .

Figure 7

Stress distributions in N1 and N3 GBs. (a) N1 GB structure after relaxation; (b) stress σ in the x-direction of the N1 GB; (c) stress σ in the y-direction of the N1 GB; (d) N3 GB structure after relaxation; (e) stress σ in the x-direction of the N3 GB; and (f) stress σ in the y-direction of the N3 GB.

Because of the different influence of GBs on κ and κ, the dramatic variation of the anisotropy ratio of thermal conductivities, α = κ/κ, has been calculated and is shown in Figure c. With the increase of the defect density/misorientation angle, α decreases from 5.5 to 1.41, from 5.8 to 2.2, and from 3.6 to 1.8 in phosphorene with type I, II, and III GBs, respectively.

One interesting feature observed here is the strikingly high anisotropy ratio α in GB-contained phosphorene. The thermal conductivities of pristine phosphorene are 83 W m–1 K–1 (zigzag) and 28 W m–1 K–1 (armchair) from frequency-dependent time-domain thermoreflectance experimental measurements,[9] 110.7–152.7 W m–1 K–1 (zigzag) and 33.0–63.6 W m–1 K–1 (armchair) from MD calculations,[11,43] and 83.5–110 W m–1 K–1 (zigzag) and 24.3–36 W m–1 K–1 (armchair) from first-principles calculations.[44−47] However, including the GB, the N1 model in type I GBs and the N6 model in type II GBs have a value of α = 5.5 and 5.8, respectively, which is remarkably higher than that in pristine phosphorene (α < 3.5).[9,11,43−49] These results indicate that phosphorene, even with GBs, is able to further enhance the thermal conduction anisotropy and thus provide a unique platform for exploring the anisotropic heat transport in 2D materials.

Conclusions

In the present work, we have investigated the GB effect on the thermal transport in phosphorene by using MD simulations. A total of 19 different types of GBs have been explored. A remarkably high Kapitza thermal resistance at the boundary has been observed, which can be attributed to the strong phonon-boundary scattering at the GBs. It is found that the defect density along the GBs has a significant effect on the thermal boundary resistance, and thus, the defect density can be used to modulate the thermal boundary resistance of phosphorene. Furthermore, the dependence of the thermal conductivity anisotropy on GBs in phosphorene has also been investigated. A strikingly high anisotropy ratio of thermal conductivities is observed, and this strong anisotropy is found to arise from the different influences of GBs on the thermal transport along the zigzag and armchair directions of phosphorene. Our results highlight that the GB effect on the thermal property of phosphorene is an important factor and thus should be taken into account in the applications of phosphorene-based nanoelectronic devices because effective thermal management of localized hot spots in nanoelectronic devices remains challenging.