**3. Effect of different physical factors**

#### **3.1. Length dependence**

As for the length dependence of the thermal conductivity of silicene, regarding the conver‐ gence with respect to the *q* mesh, it seems that for infinite size (without boundary scattering), the thermal conductivity tends to diverge with denser *q* mesh, while for finite size the thermal conductivity converges. The possible divergence of thermal conductivity for 2D materials has raised a lot of debate recently. For example, in [23] it was claimed that the thermal conductivity of silicene would diverge with the sample size. Similar conclusions have been drawn in some other literature for graphene [57]. However, there have been strong debates on the possible divergence of thermal conductivity of graphene. On the other hand, Fugallo et al. argue that the thermal conductivity of graphene converges when the simulated sampling size goes up to 1 mm [58]. In their work, exact phonon BTE is solved and first-principles calculations are employed to extract harmonic and anharmonic IFCs. Barbarino et al. also reach the same conclusion with approach-to-equilibrium molecular dynamics simulations for graphene sample of 0.1 mm in size [59]. By using a finite *q* mesh, the extremely long wavelength acoustic phonon modes are actually excluded, which are believed to be responsible for the possible divergence of the thermal conductivity [57]. For finite sample size, the boundary scattering imposes a limit on the phonon mean free path (MFP) to avoid divergence. For real applications, a finite sample has to be used, and the wrinkles and defects are generally unavoidable, so the sample cannot have diverged thermal conductivity.

As for monolayer phosphorene, Zhu et al. claimed that there exists a compelling coexistence of size-dependent and size-independent thermal conductivities along zigzag and armchair directions, respectively [53], which is distinctly different from isotropic and divergent thermal conductivities in two-dimensional graphene and silicene. The size-dependent thermal conductivity is due to the lifetime of acoustic phonon modes quickly blowing up approaching the Γ point of the Brillouin zone [57]. The contributions of FA, TA, LA, and optical phonon modes to the total thermal conductivity are plotted to further understand the dependence of *κ*zz on the size of the *q* mesh. It is found that the thermal conductivities contributed by FA, LA, and all optical phonons do not change with the increasing size of the *q* mesh. The divergent *κ*zz is mainly due to the contribution from TA branch.

#### **3.2. Strain**

Despite recent efforts to describe the properties of unstrained materials, in real applications, nanoscale devices usually contain residual strain after fabrication [60]. It is thus important to

investigate possible strain effects on the properties of materials, especially phonon transport properties.

which thus leads to its small contribution to the thermal conductivity. Note that FA contributes to the thermal conductivity along zigzag direction more than that along armchair direction. The reason might lie in the feature of the hinge-like structure that it is more "uneven" along

As for the length dependence of the thermal conductivity of silicene, regarding the conver‐ gence with respect to the *q* mesh, it seems that for infinite size (without boundary scattering), the thermal conductivity tends to diverge with denser *q* mesh, while for finite size the thermal conductivity converges. The possible divergence of thermal conductivity for 2D materials has raised a lot of debate recently. For example, in [23] it was claimed that the thermal conductivity of silicene would diverge with the sample size. Similar conclusions have been drawn in some other literature for graphene [57]. However, there have been strong debates on the possible divergence of thermal conductivity of graphene. On the other hand, Fugallo et al. argue that the thermal conductivity of graphene converges when the simulated sampling size goes up to 1 mm [58]. In their work, exact phonon BTE is solved and first-principles calculations are employed to extract harmonic and anharmonic IFCs. Barbarino et al. also reach the same conclusion with approach-to-equilibrium molecular dynamics simulations for graphene sample of 0.1 mm in size [59]. By using a finite *q* mesh, the extremely long wavelength acoustic phonon modes are actually excluded, which are believed to be responsible for the possible divergence of the thermal conductivity [57]. For finite sample size, the boundary scattering imposes a limit on the phonon mean free path (MFP) to avoid divergence. For real applications, a finite sample has to be used, and the wrinkles and defects are generally unavoidable, so the

As for monolayer phosphorene, Zhu et al. claimed that there exists a compelling coexistence of size-dependent and size-independent thermal conductivities along zigzag and armchair directions, respectively [53], which is distinctly different from isotropic and divergent thermal conductivities in two-dimensional graphene and silicene. The size-dependent thermal conductivity is due to the lifetime of acoustic phonon modes quickly blowing up approaching the Γ point of the Brillouin zone [57]. The contributions of FA, TA, LA, and optical phonon modes to the total thermal conductivity are plotted to further understand the dependence of *κ*zz on the size of the *q* mesh. It is found that the thermal conductivities contributed by FA, LA, and all optical phonons do not change with the increasing size of the *q* mesh. The divergent

Despite recent efforts to describe the properties of unstrained materials, in real applications, nanoscale devices usually contain residual strain after fabrication [60]. It is thus important to

armchair direction because of the up and down of the sublayers chains.

206 Two-dimensional Materials - Synthesis, Characterization and Potential Applications

**3. Effect of different physical factors**

sample cannot have diverged thermal conductivity.

*κ*zz is mainly due to the contribution from TA branch.

**3.2. Strain**

**3.1. Length dependence**

As for graphene, the application of tensile strain will give three noticeable changes that can affect *κ*: (1) lowered TO and LO phonon branches; (2) linearized zone-center ZA phonon branch; (3) weakened anharmonic IFCs [11]. Full solution of the BTE gives *κ* relatively strain independent up to 1% isotropic strain, which is a balanced result of the following aspects. With the increase of tensile strain, the ZA branch near the zone center becomes linear with increasing slope, leading to the increase of near zone center ZA phonon velocities and the decrease of density of ZA phonons. The character of the acoustic phonon-phonon scattering also changes due to these changes of the phonon dispersion. In addition, the tensile strain weakens the anharmonic IFCs, which acts to increase *κ*, but this is balanced by the changes in the phonon dispersion that tend to decrease *κ*. Furthermore, the lowered TO and LO branches play only a minor role in reducing *κ* through increased *acoustic* + *acoustic* ↔ *optic* three-phonon scatter‐ ing.

It was found that a mechanical tensile strain less than 5% could tune the electronic structure of silicene [61] and larger tensile strain (7.5%) could induce a semimetal-metal transition [62]. On the other hand, using first-principles it has been demonstrated that the silicene structure remains buckled even when 12.5% tensile strain is applied [62, 63]. In comparison to the structural and electronic properties, the strain effect on the lattice thermal conductivity of silicene is also investigated. Pei et al. [64] and Hu et al. [22] investigated the effect of uniaxial strain on the thermal conductivity based on the classical NEMD method. Pei et al. studied tensile strain up to 12% and concluded that the thermal conductivity first increases slightly (around 10% increment) and then decreases with an increased amount of tensile strain. Hu et al. found that the thermal conductivity of silicene sheet and silicene nanoribbon experiences monotonic increase by a factor of two with tensile strain up to 18%. The modified embeddedatom method (MEAM) [65] and original Tersoff potential [66] were used in their simulations, respectively. However, both potentials are developed for bulk silicon, thus directly applying those potentials to the new 2D silicene structure is questionable. For example, the Tersoff potential cannot even reproduce the buckled structure of silicene and the MEAM potential seems to overestimate the buckling distance. It is well known that the interatomic potential directly determines the quality of classical molecular dynamics simulation. Therefore, in order to precisely predict the strain effect on the lattice thermal conductivity of silicene and identify the underlying mechanism, it is necessary to calculate the lattice thermal conductivity of silicene under different strains using a more accurate method.

The strain-dependent thermal conductivity of monolayer silicene was recently studied based on single mode RTA and iterative solution of the BTE by Xie et al. [28], where the harmonic and anharmonic force constants are determined using first-principles calculations. Both methods yield a similar trend in the change of thermal conductivity with respect to tensile strain. It is shown that within 10% tensile strain, the thermal conductivity of silicene first increases dramatically and then decreases slightly. The maximum thermal conductivity was found when 4% tensile strain was applied, and the value was about 7.5 times that of the unstrained case. Such a dramatic change is quite unusual for solid materials, and could be used as a thermal switch together with thermal diodes to build thermal circuits. This trend is mainly due to the strain-dependent phonon lifetime, which is related to the variations of both harmonic and anharmonic force constants under strain. FA phonon lifetimes increase signifi‐ cantly under tensile strain because the structure becomes more planar, which leads to a large increase of their contribution to overall thermal conductivity, but is not the major reason for the significant change of overall thermal conductivity within 4% strain. The significant enhancement of thermal conductivity from 0 to 4% strain is mainly due to the reduced scattering of TA and LA phonons with FA phonons. The result suggests that other 2D materials with intrinsic buckling may have similar strain dependence of thermal conductivity, which is left for further investigation. Almost at the same time, Kuang et al. used a similar method to study the phonon transport properties of silicene with strain applied, and obtained similar results [29]. They also found a strong size dependence of *κ* for silicene with tensile strain, i.e., divergent *κ* with increasing system size. However, based on their calculations, the intrinsic room temperature *κ* for unstrained silicene converges with system size to 19.34 W/mK at 178 nm. They analyzed that the convergence behavior of thermal conductivity may be significantly affected by the out-of-plane acoustic phonon branch. The divergence of thermal conductivity with respect to system size is resulted from the linear behavior of FA at low frequencies, which is very different from the familiar quadratic nature for the corresponding branch in unstrained graphene. The origin of the size effect stems from nonzero contributions of FA modes at the long wavelength limit. Although physically this still demands further careful investigation, technologically speaking, using a larger supercell for the calculations of harmonic IFCs and the phonon dispersion should effectively suppress the divergence of thermal conductivity with respect to system size.

As for phosphorene, the strain effect on the thermal conductance was studied by Ong et al. using the first-principles-based non-equilibrium Green's function (NEGF) method, which yields the thermal transport behavior in the ballistic limit [67]. They find that the thermal conductance anisotropy with the orientation can be tuned by applying strain. In particular, the thermal conductance of phosphorene in zigzag direction is found to be enhanced when the strain is applied, but decreases when the strain is applied in the armchair direction; whereas the thermal conductance in armchair direction always decreases regardless of the strain direction.

Besides, Zhang et al. performed NEMD simulations to study the strain effect on the thermal conductivity of phosphorene [52]. The results show a clear trend that the thermal conductivity increases with the tensile strain and decreases with the compressive strain, which can be explained from the buckling deformation and is consistent with that of graphene. Moreover, the thermal conductivity along the zigzag direction increases slightly when the tensile strain is 0.01, and thereafter reaches a plateau until the strain level of 0.04, which may be attributed to the tension-induced elongation of the phosphorene sample.

#### **3.3. Thickness dependence**

For the thermal conductivity of few-layer films, there are two aspects affecting the phonon transport: (1) intrinsic properties of few layers, that is, crystal anharmonicity; (2) extrinsic effects such as phonon-boundary or defect scattering. The optothermal Raman study found that *κ* of suspended uncapped few-layer graphene decreases with increasing number of layers, approaching the bulk graphite limit, which was explained by considering the intrinsic quasi-2D crystal properties described by the phonon Umklapp scattering. As the thickness increases, the phonon dispersion changes and more phase-space states become available for phonon scattering, resulting in the decrease of thermal conductivity. The small thickness (< 4) also means that phonons do not have a transverse component of group velocity, leading to weak phonon scattering from the top and bottom boundaries, especially if constant *n* is maintained over the whole area of the flake. The boundary scattering will increase if n > 4, because the transverse component of the group velocity is not equal to zero, and it is harder to maintain constant *n* through the whole area of an FLG flake, resulting in *κ* below the graphite limit. The graphite value recovers for thicker films.

The layer effect on the thermal conductivity of phosphorene was studied by Zhang et al. based on NEMD simulations [52]. It was found that the thermal conductivities along the two in-plane directions are insensitive to the number of layers, which is in sharp contrast to that of graphene, as it was shown that the thermal conductivity in multi-layer graphene decreases with increas‐ ing layer number. The underlying physical origin of the layer-independent thermal conduc‐ tivity of multi-layer phosphorene was explained that unlike graphene with only one-atom thickness, the atoms in single-layer phosphorene are arranged in two sub-layers and formed a puckered geometry, which hinders the out-of-plane (flexural) phonon mode and thus diminishes the layer effect in multi-layer phosphorene. It is worth pointing out that the Lennard-Jones (LJ) potential is used for the van der Waals (vdW) interactions across different layers in multi-layer phosphorene in their work. Based on the studies of Qiao et al. and Hu et al. [31, 56], as confirmed with the real-space wave functions, the coupling between layers in few-layer black phosphorus is mediated by interactions that are much stronger than van der Waals (vdW) as in graphene or TMDCs. Considering the significant effect of interlayer interactions on the thermal conductivity, the use of vdW to quantify the interactions in the MD is questionable. Further works based on accurate first-principles calculations are needed to address the layer-dependent thermal conductivity of few-layers black phosphorus (phosphor‐ ene) in detail, which is ongoing in our group currently.

### **3.4. Effect of substrate**

as a thermal switch together with thermal diodes to build thermal circuits. This trend is mainly due to the strain-dependent phonon lifetime, which is related to the variations of both harmonic and anharmonic force constants under strain. FA phonon lifetimes increase signifi‐ cantly under tensile strain because the structure becomes more planar, which leads to a large increase of their contribution to overall thermal conductivity, but is not the major reason for the significant change of overall thermal conductivity within 4% strain. The significant enhancement of thermal conductivity from 0 to 4% strain is mainly due to the reduced scattering of TA and LA phonons with FA phonons. The result suggests that other 2D materials with intrinsic buckling may have similar strain dependence of thermal conductivity, which is left for further investigation. Almost at the same time, Kuang et al. used a similar method to study the phonon transport properties of silicene with strain applied, and obtained similar results [29]. They also found a strong size dependence of *κ* for silicene with tensile strain, i.e., divergent *κ* with increasing system size. However, based on their calculations, the intrinsic room temperature *κ* for unstrained silicene converges with system size to 19.34 W/mK at 178 nm. They analyzed that the convergence behavior of thermal conductivity may be significantly affected by the out-of-plane acoustic phonon branch. The divergence of thermal conductivity with respect to system size is resulted from the linear behavior of FA at low frequencies, which is very different from the familiar quadratic nature for the corresponding branch in unstrained graphene. The origin of the size effect stems from nonzero contributions of FA modes at the long wavelength limit. Although physically this still demands further careful investigation, technologically speaking, using a larger supercell for the calculations of harmonic IFCs and the phonon dispersion should effectively suppress the divergence of

208 Two-dimensional Materials - Synthesis, Characterization and Potential Applications

As for phosphorene, the strain effect on the thermal conductance was studied by Ong et al. using the first-principles-based non-equilibrium Green's function (NEGF) method, which yields the thermal transport behavior in the ballistic limit [67]. They find that the thermal conductance anisotropy with the orientation can be tuned by applying strain. In particular, the thermal conductance of phosphorene in zigzag direction is found to be enhanced when the strain is applied, but decreases when the strain is applied in the armchair direction; whereas the thermal conductance in armchair direction always decreases regardless of the strain

Besides, Zhang et al. performed NEMD simulations to study the strain effect on the thermal conductivity of phosphorene [52]. The results show a clear trend that the thermal conductivity increases with the tensile strain and decreases with the compressive strain, which can be explained from the buckling deformation and is consistent with that of graphene. Moreover, the thermal conductivity along the zigzag direction increases slightly when the tensile strain is 0.01, and thereafter reaches a plateau until the strain level of 0.04, which may be attributed

For the thermal conductivity of few-layer films, there are two aspects affecting the phonon transport: (1) intrinsic properties of few layers, that is, crystal anharmonicity; (2) extrinsic

thermal conductivity with respect to system size.

to the tension-induced elongation of the phosphorene sample.

direction.

**3.3. Thickness dependence**

Supported silicene is, in fact, a common form in realistic applications. For example, it can be transferred onto an insulating substrate and gated electrically. The effects of different sub‐ strates on the thermal transport of silicene were studied by Zhang et al. recently based on NEMD simulations using the optimized SW potential for silicene [68]. They observe that the thermal conductivity of silicene can be bilaterally changed with different surface crystal plane of the substrate. The phenomenon found here is fundamentally different from our general understanding of monolayer graphene supported on substrate, where the substrate always has a negative effect on the in-plane thermal transport. The discrepancy between monolayer graphene and silicene can be explained in terms of different effects induced by the substrate on the phonon transport, specifically, the competition between the out-of-plane flexural modes and the in-plane modes. This mechanism is further linked to the different atomic structure, i.e., for graphene, it is planar (no buckling distance), while silicene has a buckling distance of about 0.42 Å. By performing phonon polarization and spectral energy density (SED) analysis, the authors further revealed the underlying physics of the novel phenomenon in terms of the different impacts on the dominant phonons in the thermal transport of silicene induced by the substrate. These results indicate that by choosing different substrates, the thermal conductivity of 2D silicene can be largely tuned, which paves the way for manipulating the thermal transport properties of silicene for future emerging applications.

Very recently, the thermal conductivity of silicene supported in an amorphous silicon dioxide (SiO2) for temperature ranging from 300 to 900 K was studied by Wang et al. from MD simulations [69]. They found that the thermal conductivity of silicene has a substantial reduction with increasing temperature, and putting silicene on amorphous SiO2 leads to 78% reduction in the overall thermal conductivity of silicene at room temperature. They further compared model-level phonon properties, such as phonon relaxation times and phonon mean free paths (MFPs) of freestanding and supported silicene at 300 K. It is found that the phonon relaxation time in the case of supported silicene is reduced from 1–13 ps to 1 ps, and corre‐ sponding MFPs decrease from 10–120 nm to 0–20 nm. The thermal conductivities of free‐ standing and supported silicene are mainly (more than 85%) contributed by the longitudinal and transverse acoustic phonons, while the out-of-plane acoustic phonons have a negligible contribution of less than 3%. These results are in line with those found previously [68].
