Lattice thermal conductivity of silicon monolayer in biphenylene network

Since the successful exfoliation of graphene from bulk graphite, the research of two-dimensional (2D) materials has attracted great attention because of their unique mechanical, thermal, photoelectric, piezoelectric, and thermoelectric properties.^1–4 Due to these unique characteristics, 2D materials have a very wide range of applications, such as healthcare, energy storage, supercapacitors, and so on.^5–8 For example, the zero bandgap and absorption properties of graphene allow it to have a wide range of applications in ultrafast photonics and supercomputers.^9,10 Due to its large surface area, electrical conductivity, capacity to immobilize different molecules, and high electron transfer rate, graphene has been widely used to design biosensors.^11,12 After graphene, two-dimensional materials of other homologous elements, such as hexagonal (honeycomb) silicene,^13,14 germanene, and stanene, were theoretically predicted and studied.^15–20 Also, many materials have been prepared and synthesized experimentally. For instance, it has been reported that Li et al. prepared large area graphene films by chemical vapor deposition method.²¹ Liu et al. successfully prepared hexagonal silicene by depositing Si atoms onto the Ag(111) surface.²² Germanene was also successfully prepared on the surface of Al(111) by Wang and Uhrberg.²³ 

Recently, Fan et al. synthesized quasi-one-dimensional carbon nanoribbons in a biphenylene network on the surface of Au(111),²⁴ which consists of octagonal (o), hexagonal (h), and square (s) structures, referred to as ohs structures. This is a very important achievement, as it is a non-hexagonal, metallic, nanoscale carbon monolayer structure predicted in various theoretical studies.^25–28 Salih et al. verified the structural stability of the C ohs monolayer based on molecular dynamics analysis and found its strong metallic properties. The C ohs structures are highly resistant to the penetration of oxygen atoms and molecules.²⁹ Therefore, it can be used for anticorrosive coatings on solid surfaces. Moreover, the C ohs monolayer and its nanoribbons have great potential for future nanoelectronics as well.

Silicon, a homolog of carbon, is also predicted by theory to have an ohs structure.²⁹ However, it should be noted that the ohs structure of silicon is a low-buckled structure, not planar as in C ohs monolayer. It has been shown that it can also form quasi-1D nanotubes and nanoribbons with different edge geometries.²⁹ In specific bilayers and multilayers, strong vertical bonds are formed between the layers, and a network of connected octagonal and square rings can be formed between the layers. Moreover, its vertical heterostructure can form an effective metal junction and Schottky barrier, which can be used in electronic devices such as diodes.^29,30 It is important to note that the various applications of a 2D material in nano-electronics are closely related to its thermal transport properties, which calls for systematic and deep investigation of its phonon transport properties. Although the stability and electronic properties of Si ohs monolayer are reported detailedly, research on thermal transport properties is still lacking.

In this work, the lattice thermal conductivity κ of Si ohs monolayer was investigated by solving the Boltzmann transport equation (BTE) based on ab initio density functional theory (DFT). The comparison with C ohs monolayer, graphene, and hexagonal silicene is presented. The room temperatures κ of Si ohs monolayer along the two crystallographic directions are 2.46 and 3.25 W m⁻¹ K⁻¹, respectively. In addition, we also analyzed the mechanical properties of the material. We carefully examined the percentage contribution, phonon group velocity, and relaxation time of each phonon branch, to understand the underlying mechanism behind the low κ of Si ohs monolayer.

All calculations were performed using the Vienna ab initio simulation package (VASP)^31,32 based on density functional theory (DFT). For the DFT calculations, the projector-augmented-wave (PAW)³³ pseudopotentials and Perdew–Burke–Ernzerhof (PBE)³⁴ exchange-correlation functionals were chosen. The plane-wave energy cutoff was set to 319 eV, which is 30% higher than the recommended value. A vacuum region of 18 Å was added perpendicular to the crystal plane to avoid interactions between adjacent layers. The energy convergence was 10⁻⁸ eV and the maximal residual Hellmann Feynman forces were reduced to 10⁻⁴ eV/Å, while a 4 × 6 × 1 k-mesh was used during the structural relaxation. For structural relaxation, we performed convergence tests for energy cutoff and k-mesh sampling points, as detailed in Fig. S1 of the supplementary material.

The harmonic second-order interatomic force constants (IFCs) were obtained by the finite displacement method as implemented using the PHONOPY code.³⁵ A 5 × 5 × 1 supercell and 1 × 1 × 1 Monkhorst–Pack k-mesh³⁶ were used to ensure convergence. For the selection of supercells in the phonon spectrum calculation, we performed tests, as detailed in Fig. S2 in the supplementary material. The calculation of the thermal conductivity requires not only the harmonic IFCs mentioned above but also the anharmonic third-order IFCs. AlmaBTE code³⁷ was used to calculate the third-order IFCs and thermal transport properties. It is a software package for solving the spatially and temporally dependent Boltzmann transport equations for phonons, using only ab initio calculations as input. After careful testing, the interaction cutoff distance of 0.80 nm was chosen to ensure the convergence of lattice thermal conductivity.

The optimized Si ohs monolayer is shown in Fig. 1, which is consisting of octagons, hexagons, and squares. The primitive cell of Si ohs structure contains six silicon atoms, with the lattice constants of a = 7.13 Å, b = 5.73 Å, and the total buckling height of 1.04 Å, which are in good agreement with the previous work.²⁹ Different from the planar geometry of C biphenylene, the 2D Si ohs structure has a complex buckling configuration, with its atoms located on four different lateral planes. Thus, the structure has three different buckling distances, as shown in Fig. 1 and listed in Table I. This structure also has four different bond lengths, as listed in Table I. In addition, there are several different bond angles in this structure, which are 151.312° and 138.843° in the octagon; 126.081° and 97.685° in the hexagon, and all 90° in the square.

Figure 2 shows the calculated phonon dispersion along Γ-X-L-Y-Γ for Si ohs monolayer. As there are six silicon atoms in a primitive cell, 18 phonon branches can be found in the phonon spectra, including three acoustic phonon branches and 15 optical phonon branches. The phonon spectra have no negative frequencies, indicating that the structure of the material is dynamically stable (see supplementary material Fig. S3 for a more detailed zoomed-in partial view). Moreover, the longitudinal acoustic branch (LA) and the transverse acoustic branch (TA) of the Si ohs monolayer are linear near the Γ point, while the out-of-plane acoustic branch (ZA) deviates from the linearity close to the Γ point due to sufficiently weak interfacial interactions. It shows quadratic dependency, which is a very common phenomenon in 2D materials.^38–40 It can be seen that the LA branch couples with the low-frequency optical branches in the range of 2.21–2.90 THz. Note that this acoustic-optical coupling generally suppresses thermal transport.

To further verify the mechanical stability of the Si ohs monolayer, we calculated the elastic constants of the material. There are four independent elastic constants: C₁₁ = 39.944 N/m, C₂₂ = 38.372 N/m, C₁₂ = 16.132 N/m, and C₆₆ = 0.017 N/m. They satisfy the Born–Huang criterion⁴¹ C₁₁C₂₂ – C₁₂² > 0 and C₆₆ > 0, indicating the structure is mechanically stable.

The mechanical properties of a material such as elastic modulus and Poisson’s ratio, are key factors that must be considered in its application. In this study, the anisotropic Young’s modulus, shear modulus, and Poisson’s ratio of Si ohs monolayer were calculated as presented in Table II. Unlike lots of anisotropic material structures, mechanical properties, such as Young’s modulus and Poisson’s ratio of the material, are essentially isotropic. Compared to Young’s modulus of graphene (∼340 N/m)^42,43 and hexagonal silicene (∼61 N/m),^42,44 Si ohs monolayer has a lower stiffness, is prone to deformation, and is more flexible. Note that large elastic properties usually imply strong harmonic properties, such as high phonon group velocity and Debye temperature. Poisson’s ratio is related to the volume change during uniaxial deformation and provides into the brittle/ductile behavior of materials. The value is usually from −1 to 0.5,⁴² while the higher the value of Poisson’s ratio, the more plastic the material is. The Poisson’s ratios of about 0.4 indicate that the Si ohs monolayer has good ductility.

To obtain the κ for a 2D material, an effective thickness must be defined. Here, we define the buckling height of Si ohs monolayer plus the van der Waals diameter of the silicon atom as its effective thickness, which is 5.236 Å. Figure 3(a) shows that the κ varies with the cutoff distance at 300 K. In addition, it can be found that the κ starts to converge at 0.55 nm. Here, we have calculated the κ of the material as a function of temperature using two methods, the relaxation time approximation (RTA) and the full solution of the Boltzmann transport equation (BTE). The calculated κ of Si ohs monolayer in the full solution (the RTA) is 2.46 (2.50) W m⁻¹ K⁻¹ along the a-axis direction and 3.25 (2.91) W m⁻¹ K⁻¹ along the b-axis direction, at 300 K. The κ along b-axis is about 31.68% (16.55%) higher than the one along a-axis in the full solution (RTA). The result of the full solution is very close to the one of RTA, suggesting that the normal processes affect thermal transport little. It should be noted that our discussion is based on the results of the full solution, as the full solution is more reasonable theoretically mainly. Note that such low κ would be particularly beneficial for thermoelectrics. The κ was found to have a T⁻¹ dependence, which confirms that the Umklapp processes play a dominant role in phonon thermal transport in the considered temperature range. In addition, due to its anisotropic structure, the κ of the Si ohs monolayer has a significant anisotropy.

The Debye temperature (Θ_D) is an important fundamental parameter related to harmonic properties closely, which can be calculated from the highest acoustic phonon frequency ω_m with the expression Θ_D = hω_m/k_B, where h is the Planck’s constant and k_B is the Boltzmann’s constant. The calculated Θ_D for Si ohs monolayer is 138.31 K, which is a quite low relative to other 2D materials, such as graphene (∼2300 K),⁴⁵ hexagonal silicene (∼798 K),⁴⁶ and C ohs monolayer (∼1800 K).⁴⁷ It shows a positive correlation with elastic properties. Si ohs monolayer has a lower Young’s modulus than graphene, C ohs monolayer, and hexagonal silicene, hence the lowest Debye temperature. It is also in agreement with the commonly accepted criteria that a low Θ_D usually leads to a low κ.

As presented in Table III, we present the information on the lattice constants and lattice thermal conductivities of several 2D materials. It can be found that the κ of Si ohs monolayer is one order of magnitude lower than that of graphene and also significantly lower than that of hexagonal silicene and C ohs monolayer. Moreover, for the 2D material of the same element, the κ decreases with the increase of its lattice constant as the decrease of interatomic forces.

The percentage contribution to the total κ of each phonon mode is calculated. The percentage contribution of ZA, TA, LA, and optical branches along the a-axis (b-axis) is 2% (3%), 31.2% (13.7%), 16.6% (23.5%), and 52% (62.5%), respectively. In general, the thermal transport of materials is mainly determined by acoustic phonons. However, the acoustic branches of the Si ohs monolayer contribute less than 50% to the total κ, while the optical branches have the major contribution. A similar phenomenon has been reported in previous studies on some other materials, such as SnSe,⁵¹ AlN,⁵² KCuSe,⁵³ and TiS₃.⁵⁴ The ZA mode has a negligible contribution among the three acoustic branches in the Si ohs monolayer. In contrast, in graphene, the contribution of the ZA branch is dominant, which is known to be about 80%.⁵⁵ Due to the planar structure of graphene, the mirror symmetry allows only an even number of ZA phonons for the phonon–phonon scattering processes.^56,57 This makes the relaxation time of the ZA branch much longer than those of the TA and the LA modes, leading to the high contribution of the ZA mode. However, the low buckling (not planar) structure in the Si ohs monolayer has no mirror symmetry, so the contribution of the ZA mode is much lower than that of graphene.

The lattice thermal conductivity of the material can be expressed as⁵⁸: $καα=∑q⃗pCVvαq⃗,p2τq⃗,p$⁠, where C_V is the volumetric specific heat capacity, $vα(q⃗,p)$ is the α(= x, y, z) component of the phonon group velocity of phonon mode with wave vector $q⃗$ and polarization p, and τ is the relaxation time. Based on the equation, it is found both the phonon group velocity (harmonic property) and relaxation time τ (anharmonic property) have tremendous roles in affecting κ. Thus, we can further understand the thermal transport properties by phonon group velocity and τ.

The phonon group velocity is calculated as listed in Table IV and shown in Fig. 4. Although weighted average value about phonons was used in some works^59–61 to analyze the thermal transport properties, the simple arithmetic average value is adopted here as it also works well and shows the availability in many previous studies.^62–65 It can be seen that the maximum value of phonon group velocity on the b-axis is higher than that on the a-axis for all branches. Among the average values of phonon group velocity, only the one of TA modes along the b-axis is lower than that along the a-axis. It indicates that the higher κ along the b-axis than the a-axis. On the whole, these average phonon group velocities don’t differ much for all phonon branches. Moreover, compared to the high phonon group velocity of graphene (up to 23 × 10³ m/s),⁶⁶ hexagonal silicene (up to 9.2 × 10³ m/s)⁶⁷ and C ohs monolayer (up to 19 × 10³ m/s),⁴⁹ the phonon group velocity of Si ohs monolayer is much lower. As mentioned above, Si ohs monolayer has the lowest Young’s modulus and Θ_D among these monolayers, usually indicating the lowest phonon group velocity and κ.

The relaxation time τ at room temperature of all phonon branches is studied, as shown in Fig. 5. Most TA and LA phonons have a long τ in the range of 0.1–10 ps, while ZA phonons have very short τ in the range of 10⁻³ to 1 ps. The τ of most optical phonons is in the range of 0.1–1 ps, a bit shorter than TA and LA phonons but significantly longer than ZA phonons. Specifically, the average relaxation times for the ZA, TA, LA, and optical modes are 0.06, 1.46, 1.06, and 0.76 ps, respectively. Note that the phonon group velocities are close to each other as listed in Table IV. Thus, τ mainly determines the different percentage contributions for each phonon branch. ZA phonons have quite short τ, leading to a little contribution to κ, while the long τ of TA and LA branches result in a significant contribution. For the optical branches, although the τ are slightly shorter than LA and TA branches, the number of phonon modes is much larger, finally leading to the major contribution. Furthermore, the τ of Si ohs is lower than many other 2D materials, such as graphene (up to 10⁶ ps order of magnitude)⁶⁸ and C ohs monolayer (up to 10⁵ ps order of magnitude).⁶⁸ It is also responsible for the low κ of Si ohs monolayer. The τ of hexagonal silicene (up to 10² ps order of magnitude)⁶⁹ and Si ohs are in the same order of magnitude. However, smaller phonon group velocity mainly leads to the lower κ in Si ohs monolayer than hexagonal silicene.

In practical applications, κ of material is suppressed as any material has a finite size, especially at the nanoscale. Here, the size effect of thermal transport is estimated by the cumulative κ, where the contribution of only phonons with mean free path (MFP) below a threshold is considered. The normalized cumulative κ along both directions with respect to the phonon MFP at 300 K are exhibited in Fig. 6. It can be found that the MFP of phonons with significant contribution is in the range of 1–10² nm along the two directions. In fact, the normalized cumulative κ along the b-axis rises slightly more quickly than the one along the a-axis. The phonons with MFP below 10² nm contribute 84% and 94% of the total κ along the a- and b-axes. When the MFP is smaller than 10 nm, the cumulative κ along both two directions declines sharply.

In summary, based on first-principles calculations, we have investigated the κ of Si ohs monolayer by solving the phonon BTE with AlmaBTE. The phonon dispersion and mechanical properties are studied as well. It is found that the stiffness of Si ohs monolayer is much smaller than that of graphene, C ohs, and hexagonal silicene. It indicates weak harmonic properties, such as lower Θ_D and phonon group velocity. At 300 K, the κ of Si ohs monolayer is 2.46 W m⁻¹ K⁻¹ in the a-axis and 3.25 W m⁻¹ K⁻¹ in the b-axis, showing a significant anisotropy. It is found that the anisotropy of phonon group velocity is responsible for the anisotropy of the κ. The different percentage contributions mainly originate from the τ of phonon branches, as their phonon group velocities are close to each other. The contribution of ZA phonons is negligible as the very short τ. The LA and TA modes contribute to the κ in a non-negligible fraction, as the τ of them are the longest. Although the τ of optical modes is a little shorter than that of LA and TA modes, the contribution of optical modes still exceeds 50% as the number of optical modes is much more than acoustic modes. Compared to graphene, C ohs monolayer, and hexagonal silicene, the weakest harmonic property (lowest Debye temperature and phonon group velocity) and strongest anharmonic property (shortest relaxation time) are the reasons for the lowest κ of Si ohs monolayer. Finally, the size effect is studied. The κ will decline sharply when the size of the sample is smaller than 10 nm. Our findings reveal the underlying physical mechanisms of thermal transport properties for Si ohs monolayer and pave the way for future applications in thermoelectric and thermal management.

The supplementary material contains convergence tests for energy cutoff and k-mesh taking points at structural relaxation, and supercell size of phonon spectrum calculation, and gives a zoomed-in view of the phonon spectrum in the paper.

This work was supported by the National Natural Science Foundation of China (Grant No. 11974100 and 61874160) and the Program for Innovative Research Team (in Science and Technology) at the University of Henan Province (Grant No. 22IRTSTHN012).

The authors have no conflicts to disclose.

Aiqing Guo: Data curation (lead); Formal analysis (lead); Funding acquisition (lead); Investigation (equal); Methodology (equal); Project administration (lead); Supervision (equal); Validation (equal); Writing – original draft (lead); Writing – review & editing (equal). Fengli Cao: Data curation (lead); Formal analysis (lead); Investigation (equal); Methodology (equal); Supervision (equal); Validation (equal); Writing – original draft (lead); Writing – review & editing (equal). Weiwei Ju: Funding acquisition (equal); Methodology (equal); Project administration (equal); Supervision (equal); Writing – review & editing (equal). Zhaowu Wang: Funding acquisition (equal); Project administration (equal); Supervision (equal). Hui Wang: Funding acquisition (equal); Project administration (equal); Supervision (equal). Guo-Ling Li: Funding acquisition (equal); Project administration (equal); Supervision (equal). Gang Liu: Data curation (equal); Formal analysis (equal); Funding acquisition (lead); Investigation (equal); Methodology (lead); Project administration (lead); Supervision (equal); Validation (equal); Writing – review & editing (lead).

The data that support the findings of this study are available from the corresponding author upon reasonable request.