Prediction of novel ground-state structures and analysis of phonon transport in two-dimensional Ge_xS_y compounds

Abstract

We conducted this study to explore the ground-state structures of two-dimensional (2D) variable-composition Ge_xS_y compounds, driven by the polymorphic characteristics of bulk germanium sulfides and the promising thermoelectric performance of 2D GeS (Pmn2₁). To accomplish this, we utilized the highly successful evolutionary-algorithm-based code USPEX in conjunction with VASP for total energy calculations, leading to the discovery of three previously unexplored structures of Ge₂S (P2/c), GeS (P3̄m1), and GeS₂ (P2₁/c). These 2D materials exhibit significantly lower formation energies compared to their reported counterparts. We thoroughly scrutinized the structural stability and subsequently analyzed their electronic structures. Our analysis reveals a nearly direct band gap of 0.12/0.84 eV with the PBE/HSE06 functional for 2D Ge₂S and an indirect band gap for 2D GeS and GeS₂. Their semiconducting nature highlights the crucial importance of lattice thermal conductivity (κ_l), which we determined by solving the Boltzmann transport equation for phonons. Importantly, we predict a room temperature κ_l value of 6.82 W m⁻¹ K⁻¹ for GeS, lower than its 2D orthorhombic counterpart. In the case of GeS₂, we observed an anisotropic κ_l value of 16.95/10.68 W m⁻¹ K⁻¹ along the zigzag/armchair directions at 300 K, with an in-plane anisotropy ratio of 1.59, surpassing that of 2D IV–VI compounds. We delve into detailed discussions regarding the role of lattice anharmonicity, group velocities, phonon lifetimes, and three-phonon weighted phase space in the overall thermal conductivity analysis.

---

1 Introduction

The demand for two-dimensional (2D) materials has significantly increased due to their unique properties and potential for miniaturization in future devices. One such material is graphene, derived from graphite, which exhibits exceptional characteristics such as a high carrier mobility and lattice thermal conductivity (κ_l), making it suitable for efficient thermal transport.^1–3 However, in the context of thermoelectric devices, the efficiency at temperature T is quantified by a dimensionless figure of merit, , where σ, S, κ_e, and κ_l represent the electrical conductivity, Seebeck coefficient, electronic thermal conductivity, and lattice thermal conductivity, respectively. It is worth noting that ZT is inversely proportional to the κ_l. Therefore, the high κ_l and metallic nature of intrinsic graphene pose challenges for its direct use in thermoelectric devices, and various methods like nanostructuring and doping are employed to mitigate these issues.⁴ Consequently, alternative materials are being explored, with 2D transition metal dichalcogenides (TMDs) gaining considerable attention. Nonetheless, most TMDs exhibit relatively high lattice thermal conductivity,⁵ prompting researchers to redirect their focus toward 2D IV–VI materials. These materials have emerged as a prominent research area due to their ability to achieve ultra-low lattice thermal conductivity and higher Seebeck coefficients.⁶

The suitable band gap, earth abundance, and non-toxic characteristics of germanium sulfides and its constituent elements make the orthorhombic (Pnma) phase of GeS the primary focus of researchers. The monolayer GeS (Pmn2₁) can be obtained through mechanical exfoliation from its bulk and has shown great potential for thermoelectric device applications due to its high electron mobility and ultra-low lattice thermal conductivity. Similarly, the 1T structure of 2D GeS₂⁷ and the low-symmetric (P2₁2₁2) phase of two-dimensional (2D) Ge₂S⁸ have also exhibited excellent properties, particularly as thermoelectric and anode materials. However, while the bulk counterparts of these Ge_xS_y compounds exhibit various polymorphs, such as cubic,⁹ orthorhombic,¹⁰ and amorphous¹¹ phases of GeS, and monoclinic,¹² orthorhombic,¹³ and beta¹⁴ phases of GeS₂, it remains unexplored whether their monolayer forms possess unique and undiscovered structures. This intriguing possibility suggests that the monolayers of these compounds may exhibit distinct properties, thereby warranting further investigation and exploration.

In the pursuit of obtaining 2D counterparts of materials, researchers employ various experimental procedures. These methods can be of two categories: top-down and bottom-up approaches. In the first approach, a single layer is exfoliated from the bulk material, which is particularly feasible for layered materials held together by van der Waals forces such as graphene¹ or 2D TMDs. Conversely, in the bottom-up approach, layers grow on a substrate from vapors of precursors through methods like physical or chemical vapor depositions.¹⁵ Bottom-up approaches are employed for synthesizing 2D sheets from non-layered bulk materials due to their high exfoliation energy. They have also enabled the synthesis of 2D structures that lack bulk counterparts. For instance, various 2D boron allotropes have been successfully prepared on substrates such as silver surfaces.¹⁶ The 2D hexagonal layer of GeSe, synthesized through chemical vapor deposition,¹⁷ shares a similar structure with our predicted 2D GeS in this work and also lacks a bulk counterpart. In parallel, density functional theory (DFT) can predict a single layer by cutting it from the layered bulk and allowing for relaxation with sufficient vacuum in the out-of-plane direction. In contrast, the non-layered bulk materials lack preferred planes for cleaving a single layer, leading to many possible 2D structures along various planes. Among these, the one with the lowest formation energy is energetically favorable, called the ground-state 2D structure of the given composition. The ground state structures of materials have particular importance due to their higher possibility of experimental realization. Another strategy involves replacing atoms in already explored monolayer lattices with the elements under investigation and searching for the most energetically favorable configuration, such as in the case of germanene and silicene.¹⁸ However, manually exploring all possible atomic arrangements can be a tedious task. Fortunately, significant advancement in computational tools for predicting ground state structures has revolutionized the field of materials science, allowing the prediction of novel 2D materials like borophenes before their experimental synthesis.^16,19 Among these tools, the Universal Structure Predictor: Evolutionary Xtallography (USPEX)^20,21 stands out with its high success rate among genetic-algorithm-based search codes. USPEX offers a specific 2D search capability, making it particularly suitable for identifying the lowest-energy phases in compounds with variable compositions, such as 2D Cd_xTe_y.²² By harnessing the power of USPEX, one can efficiently search for and identify the energetically stable structures of 2D Ge_xS_y compounds.

This study predicts novel ground-state structures of 2D Ge_xS_y compounds with three different compositions: Ge₂S, GeS, and GeS₂, that exhibit higher stability than their previously reported counterparts. Specifically, we focus on their crystal structures, electronic structure, and phonon transport properties. The 2D GeS phase demonstrates lower lattice thermal conductivity compared to its orthorhombic monolayer counterpart. We delve into the underlying mechanisms responsible for the low κ_l of 2D GeS, examining factors such as anharmonicity, phonon group velocity, their lifetimes, and weighted phase space for three-phonon scattering. Moreover, 2D GeS₂ shows intriguing in-plane anisotropic elastic stiffness and κ_l, while 2D Ge₂S has comparatively higher κ_l. The interplay of anisotropic acoustic and low-frequency optical phonon group velocity and relatively higher lifetimes to their distinguished κ_l in 2D GeS₂ makes this study of fundamental interest.

2 Computational methodology

The search for lowest-energy structures of 2D germanium sulfides involves USPEX^20,21 in conjunction with the Vienna Ab initio Simulation Package (VASP).²³ Our variable-compositional search is constrained to a layer thickness of 0–6 Å and the constituent atoms in a unit cell of 2–12. In this search, a total of 80 generations are considered. The first generation contains 200 randomly created structures, while the subsequent generations have 80 structures each that evolve through natural selection and variation operators such as mutation, crossover, and symmetry operations. The relative energetic stability of different compositions and structures was determined by comparing the formation energy per atom (E_f), which is calculated using the equation:

In this equation, E(Ge_xS_y), E(Ge), and E(S) represent the total energies of the lowest possible 2D configurations of Ge_xS_y, elemental germanium, and sulfur, respectively. The total energy of each structure is calculated using the projector-augmented wave (PAW)²⁴ method, which is well-suited for handling core electrons. To account for the electronic exchange–correlation interaction, the generalized gradient approximation (GGA)²⁵ with the Perdew–Burke–Ernzerhof parameterization (PBE) is employed. To maximize computational efficiency, a strategy is adopted that includes a relatively coarse k-mesh with a resolution of 2π × (0.08–0.04) Å⁻¹, a smaller cutoff energy of 450 eV for plane wave expansion, a vacuum size of 15 Å, and exclusion of the zero damping of Grimme (DFT-D3)²⁶ for the van der Waals corrections in the first 60 generations, followed by its inclusion in the final 20 generations.

To ensure accuracy, the resulting 2D Ge_xS_y structures undergo re-optimization using a larger cutoff energy of 540 eV, a vacuum size of 20 Å, and a denser Γ-center k-mesh with a resolution of 2π × 0.02 Å⁻¹. The convergence criteria for energy is set to 10⁻⁸ eV and for forces to 10⁻³ eV Å⁻¹. Electronic band structures are computed using both the GGA and Heyd–Scuseria–Ernzerhof (HSE06)²⁷ methods. For ab initio molecular dynamics (AIMD) simulations, we employ the NVT ensemble with a step size of 1 fs and a total time of 10 ps at temperatures of 300 K, 600 K, and 900 K. We use the Nosé–Hoover thermostat²⁸ to maintain a constant temperature during simulations. The Phonopy code²⁹ is used to calculate the phonon dispersion based on the finite displacement method with supercells. We use supercells of 6 × 5 × 1 for Ge₂S, 5 × 5 × 1 for GeS, and 3 × 3 × 1 for GeS₂, which are the same-sized supercells as those used for AIMD simulations.

We utilized an iterative self-consistent approach implemented in the ShengBTE package³⁰ to solve the phonon Boltzmann transport equation (PBTE) and predict the lattice thermal conductivity. This method relies on two sets of interatomic force constants (IFCs): harmonic (second-order) and anharmonic (third-order). The harmonic IFCs are obtained using the Phonopy code²⁹ and anharmonic IFCs through the thirdorder.py python script based on the method developed by Lindsay et al.³¹ in combination with VASP. With these sets of IFCs, one can express the αβ component of the lattice thermal conductivity tensor as:³⁰

Here, λ(k,p) represents a phonon mode characterized by the wave vector k and branch number p. Δ^β_λ with the dimensions of velocity is the correction to relaxation time approximation (RTA) prediction. V denotes the unit-cell volume, ω_λ is the phonon frequency, N is the number of q points in the first Brillouin zone (BZ) for a Γ-centered regular grid, T represents the temperature, v^α_λ is the group velocity of the phonon in the α direction, τ_λ is the λ phonon lifetime, and n₀ is the phonon occupation number given by the Bose–Einstein distribution. Here, the inverse of τ_λ, i.e., 1/τ_λ is the total scattering rate, which is determined through Matthiessen's rule given as:

In this equation, 1/τ^3ph_λ represents the three-phonon scattering rate, while 1/τ^iso_λ and 1/τ^b_λ denote the scattering rates associated with phonon–isotope and phonon–boundary interactions, respectively. For a detailed expression of each scattering rate, please refer to ref. 30.

The obtained κ_l is rescaled for 2D Ge_xS_y by multiplying a factor L_z/d_eff, where L_z is the total length of the unit cell along the z direction and d_eff is the effective thickness of the materials obtained from the corresponding bulk structure. The values of d_eff are 7.579, 9.068, and 7.710 in units Å for Ge₂S, GeS, and GeS₂, respectively (see Fig. S1 of the ESI†). The long-range interactions are included by Born effective charge tensors (see Table S1 of the ESI†) computed with the density functional perturbation theory (DFPT).³²

3 Results and discussion

3.1 Structure prediction and stability

The comparative energetic stability among various structures of a material can be represented by the energy convex hull (ECH). The ECH of the 2D Ge–S system (see Fig. 1) is constructed by plotting the formation energy per atom, given by eqn (1), of a given composition ratio in all possible structures along the y-axis and composition ratio S/(Ge + S) along the x-axis. In the literature, several polymorphs of bulk GeS, such as orthorhombic α, cubic, and amorphous phases,^9,11,33 have been reported due to different bonding modes. ECH construction is an efficient and powerful tool for identifying the most stable structure and composition,³⁴ especially for materials with numerous polymorphs. These most stable structures will lie on the convex hull, such as in the Ge–S system, Ge₂S (magenta), GeS (cyan), and GeS₂ (dark-blue), as shown in Fig. 1. These lowest-energy 2D Ge₂S, GeS, and GeS₂ (discussed later) have different structures and lower formation energies than the previously reported ones. For comparison, the formation energies of the already reported 2D Ge₂S,⁸ GeS,³⁵ and GeS₂⁷ structures are represented in Fig. 1 by magenta, cyan, and blue triangles, respectively. The previously known 2D Ge_xS_y structures, i.e., Ge₂S, GeS, and GeS₂, have higher formation energies than our predicted counterpart monolayers by 85, 47, and 84 in units of meV per atom, respectively, and consequently occupy positions above the convex hull. Worth noticing from the ECH analysis is the superior energetic stability of our predicted Ge_xS_y structures, which indicates easier experimental accessibility. Therefore, they deserve detailed stability analysis, and we constrain the discussion to the structures on the convex hull only.

Fig. 1 This figure shows the evolution of the energy convex hull (ECH) during the USPEX search, where the solid bounding curve indicates the final ECH of the 2D Ge–S system. The energetically most stable structures of 2D Ge_xS_y are highlighted with various colored circles, whereas those previously reported ones are represented with filled triangles.

The phonon band structure, which provides critical insights into the dynamic stability of the newly predicted 2D phases of the Ge₂S, GeS, and GeS₂, are determined and shown in Fig. 2. These results reveal the dynamic stability of the 2D Ge₂S, GeS, and GeS₂, as there are no imaginary frequencies in Fig. 2(a)–(c). There are a total of 3n phonon branches, where n is the number of atoms in the primitive unit cell. The first three lowest-frequency phonon branches are called acoustic phonons, which are further classified into transverse acoustic (TA), longitudinal acoustic (LA), and z-axis acoustic (ZA). The ZA mode, also called the flexural mode, represents the out-of-plane atomic vibrations and has a parabolic wave vector dependence contrary to the TA and LA, which have linear dispersion near the Γ-point. In Fig. 2(c), the LA mode along the Γ–Y direction has softer behavior (smaller slope) than along the Γ–X direction. This difference in the phonon dispersion suggests the presence of anisotropic lattice thermal conductivity, which will be discussed later. The remaining 3n − 3 phonon branches are known as optical branches, divided into low and high-frequency bands separated by a gap. This gap is 16/41 cm⁻¹ in Ge₂S/GeS and 88 cm⁻¹ in GeS₂. The common observation from Fig. 2 is the overlapping of the low-frequency optical phonons with acoustic phonons. It is one of the indicators of the high possibility of three-phonon scattering, which will eventually decrease lattice thermal conductivity. The maximum phonon frequency (ω_max) in all predicted structures is below 400 cm⁻¹. Specifically, ω_max is 365 cm⁻¹ for Ge₂S, smaller than its P2₁2₁2 phase (561 cm⁻¹),⁸ and 337 cm⁻¹ for GeS, in the same range as the Pmn2₁ phase (∼330 cm⁻¹).⁶ω_max is 397 cm⁻¹ for 2D GeS₂, comparable to its 1T phase (∼400 cm⁻¹).⁷ Additionally, the dispersion in the high-frequency optical branches is comparatively lower in GeS₂ than in Ge₂S and GeS. The non-analytical correction to phonon spectra (see Fig. S8 of ESI†) exhibits a negligible effect on the dispersion of 2D Ge₂S and GeS₂. However, in 2D GeS, a splitting of high-frequency optical modes near the Γ-point of the BZ is observed. Despite this splitting, it has a negligible impact on the calculated lattice thermal conductivity due to their shorter lifetimes, as discussed later.

Fig. 2 Phonon band structures of 2D (a) Ge₂S, (b) GeS, and (c) GeS₂. The green, blue, red, and black colors represent the longitudinal acoustic (LA), transverse acoustic (TA), flexural acoustic (ZA), and optical mode phonons, respectively.

In addition to energetic and dynamic stability, the elastic stability of the predicted 2D Ge_xS_y provides more insights into their viability. The elastic stability of our predicted 2D monoclinic Ge₂S (P2/c), GeS₂ (P2₁/c), and trigonal GeS (P3̄m1) is confirmed using the corresponding stability criteria³⁶ based on the elastic constants presented in Table 1. Specifically, GeS with a trigonal structure satisfies the criteria C₁₁ > 0, C₆₆ > 0, and C₁₁ > |C₁₂|. The elastic constants of GeS₂ satisfy the criteria for a rectangular lattice, i.e., C₁₁ > 0, C₆₆ > 0, and C₁₁ × C₂₂ > C₁₂². Similarly, due to its oblique structure, Ge₂S satisfies the criterion det|C_ij| > 0, with i, j = 1, 2, 6.

Table 1 Structural and elastic properties of 2D Ge_xS_y compounds. The Young's moduli and Poisson's ratios were calculated from the elastic constants using the corresponding formulas: Y_x = (C₁₁C₂₂ − C₁₂²)/C₂₂, Y_y = (C₁₁C₂₂ − C₁₂²)/C₁₁, ν_x = C₁₂/C₂₂, and ν_y = C₁₂/C₁₁

Structure (space group)	Lattice parameters (a, b), γ (°)	Elastic constants (N m^−1) C11, C22, C12, C66	Young's moduli Yx, Yy (in N m^−1)	Poisson's ratio νx, νy
Ge2S (P2/c)	4.98, 5.56, 45°	47.80, 54.68, 8.07, 27.96 C26 = 10.03, C16 = 5.84	46.6, 53.3	0.147, 0.169
GeS (P3̄m1)	3.63, 3.63, 120°	85.40, 85.40, 28.46, 28.47	75.9, 75.9	0.333, 0.333
GeS2 (P21/c)	6.87, 6.17, 90°	64.58, 27.84, 4.22, 12.56	63.8, 27.6	0.152, 0.065

---

Generally, 2D materials are more vulnerable to high temperatures compared to their bulk. Therefore, analyzing the thermal stability of predicted 2D Ge_xS_y compounds at elevated temperatures is worthwhile. For this purpose, AIMD simulations are performed to investigate the thermodynamic stability of newly predicted monolayers at different temperatures of 300, 600, and 900 K. The simulation results are shown in Fig. 3 at 300 K and Fig. S2 and S3 of the ESI† at temperatures of 600 and 900 K, respectively. These results reveal that 2D GeS and GeS₂ exhibit higher structural stability because no significant changes were observed in the overall morphology or bond lengths at 300, 600, and 900 K. However, the thermally stable 2D Ge₂S at 300 and 600 K becomes unstable when the temperature rises to 900 K, as evident from the broken bonds in the inset of Fig. S3(a) in the ESI.†

Fig. 3 The total energy per unit cell and temperature profiles of the 2D materials (a) Ge₂S, (b) GeS, and (c) GeS₂ were obtained through AIMD simulations at 300 K. The simulated structures were visualized from the top and side in the insets after a simulation time of 10 ps.

Under real experimental conditions, the 2D layer can interact with various molecules, including nitrogen and oxygen, which are predominant components of the air and can influence their stability. To address this concern, we assessed the stability of the predicted 2D Ge_xS_y against N₂ and O₂ molecules at 0 K. We conclude this section with the remarks that our predicted three compositions of 2D Ge_xS_y are energetically more favorable than previously reported ones, stable dynamically and elastically, and can withstand higher temperatures.

3.2 Structural analyses

The properties of a material, either bulk or two-dimensional, are entangled with its crystal structures.³⁷ Therefore, before delving into the electronic structure and phonon transport properties, it is crucial to investigate the crystal structures of the low-energy lying 2D Ge_xS_y. In Fig. 4(a)–(c), we show the 2D structures of the stable compositions Ge₂S, GeS, and GeS₂, respectively. The top and side views of each structure in these figures will facilitate the readers understanding of the discussion in this section. The dashed lines represent the primitive unit cell in each case. Regarding symmetry, Ge₂S has a monoclinic unit cell with the space group P2/c, the GeS₂ unit cell has the space group P2₁/c, while GeS adopts a trigonal structure with the space group P3̄m1. Due to the low symmetry of 2D Ge₂S, the views of various orientations are shown in Fig. 4(a). It composes two sublayers; the red-dashed and the solid-black rhombic lattices represent the bottom and top sublayers, respectively. The top sublayer can be obtained from the bottom by mirror symmetry along the z-axis plus a glide of 0.5b along the b⃑ lattice vector, where b is the lattice constant. For clarity, we show the top view of only the bottom sublayer on the right side of Fig. 4(a). It consists of irregular (unequal bonds) parallel zigzag chains of Ge atoms in the ab-plane connected through S atoms. Each Ge atom makes four bonds, one with the S atom and three with other Ge atoms, of which two Ge belong to the same sublayer and one Ge to the neighboring sublayer. Hence, the two sublayers are attached by the Ge–Ge bond of length 2.52 Å. The Ge–Ge bond lengths in the zigzag chain are 2.47 Å and 2.51 Å. Other lattice parameters and mechanical properties, such as the elastic constants and Poisson's ratios, are tabulated in Table 1.

Fig. 4 The top and side views of 2D (a) Ge₂S, (b) GeS, and (c) GeS₂ and the primitive unit cell are shown with dashed lines in each top view. All bond lengths and thicknesses of the layers are shown in units of angstroms (Å). The solid circle (red) in (a) represents the bond between two sublayers. The red-dashed lines in (c) represent zigzag and armchair directions.

Fig. 4(b) shows the top and side views of 2D GeS. It consists of two oppositely buckled honeycomb sublayers with AB stacking. In this structure, the in-plane coordinates of the S atoms are (1/3a⃑, 1/3b⃑) and of the Ge atoms are (0,0) and (1/3a⃑, 2/3b⃑), where a⃑ and b⃑ are the lattice vectors. For GeS, a = b = 3.63 Å is the optimized lattice constant. Within each sublayer, Ge/S atoms are bonded with three S/Ge atoms through a bond length of 2.44 Å. The two sublayers are then weakly connected by a Ge–Ge bond of length 2.93 Å. Fig. 4(c) shows that the primitive unit cell of 2D GeS₂ constitutes corner-sharing distorted (unequal edges) GeS₄ tetrahedra. The GeS₄-tetrahedron contains the Ge atom in its center, bonded with four S atoms at the corners. All four Ge–S bonds in the GeS₄ tetrahedron have different lengths, i.e., 2.23, 2.24, 2.26, and 2.28 in Å, which reduces the unit-cell symmetry from orthorhombic Pbcm (β-2D silica³⁸) to monoclinic P2₁/c. The side views in Fig. 4(c) exhibit the zigzag morphology along the a-axis and armchair along the b-axis. The GeS₄ tetrahedra centers align in the ac-plane and buckle in the bc-plane, which leads to a hinge-like structure along the b-axis. This arrangement in 2D GeS₂ results in a rectangular lattice with lattice parameters a = 6.87 Å and b = 6.17 Å. The structural anisotropy of the lattice is reflected in its high elastic anisotropy, where the Young's modulus along the a-axis is 63.8 N m⁻¹, which is more than twice that along the b-axis (27.6 N m⁻¹). Like phosphorene,^39,40 we expect that GeS₂ will have high anisotropic lattice thermal conductivity because it is often directly proportional to the elastic modulus.⁴¹

3.3 Electronic structure

The electronic band structure and the total density of states (DOS) with projected density of states (PDOS) for the three compositions of 2D Ge_xS_y are presented in Fig. 5. A first look shows the semiconducting nature of all these compositions. 2D Ge₂S in Fig. 5(a) depicts a nearly direct and narrow band gap of 0.12/0.84 eV using the PBE/HSE06 functional. The valence band maximum (VBM) and conduction band minimum (CBM) occur along the Γ–C–V–S path. The complete k-path used for plotting this band structure is given in Fig. S4(c) of the ESI.† Interestingly, the k-points corresponding to the band edges, denoted by (C,V), do not occur along the high symmetry k-path of the first Brillouin zone (BZ). We identified this behavior of 2D Ge₂S from Fig. S4(a) (ESI†), where the band edges of the high-symmetry k-path [Fig. S4(b), ESI†] do not match the total DOS. This discrepancy should be attributed to the band structure being limited to only high-symmetry k-points while the total DOS covers the whole first BZ. The k-path is revised, as shown in Fig. S4(c) of the ESI,† by including accurate band edges [C (0.333 0.384 0), V (0.37 0.425 0)] determined from the DOS with a denser k-mesh of reciprocal space resolution of 2π × 0.004 Å⁻¹. Due to the high dispersion of the band structure at the CBM, there is a lower density of states than in the VBM. Furthermore, one can analyze the atomic orbital contributions to the band structure by the PDOS, right panel in Fig. 5(a). This implies that the VBM of 2D Ge₂S is contributed by the p-orbitals of both Ge and S atoms, whereas the CBM mainly originates from the s-orbital of the Ge atom. The 2D GeS and GeS₂ exhibit band structures in Fig. 5(b) and (c) with indirect band gaps. Specifically, the VBM of 2D GeS is located on both sides of the Γ-point, i.e., along Γ–X and Γ–K. The CBM is located at the X point, where the band structure is highly asymmetric, as shown in Fig. 5(b). The asymmetric conduction band around the CBM is a footprint of the electrical anisotropy. The VBM of 2D GeS includes the p and s orbitals of S and Ge atoms, respectively. Its CBM is composed of Ge-s, Ge-p, and S-p orbitals. 2D GeS₂ in Fig. 5(c) has VBM at the Γ point and CBM at the Y point of the first BZ. The valence and conduction bands are flat near the band edges and have a higher density of states than the 2D Ge₂S and GeS₂. The p-orbital of the S atom mainly contributes to the VBM and CBM, and there is also a comparable contribution from the Ge-s orbital to the CBM. The fundamental band gaps, calculated with PBE functionals, are 0.72/1.85 eV for 2D GeS/GeS₂. The dispersions of band structures with the HSE06 functional remain similar to PBE, but the fundamental band gaps open to 1.05 eV for 2D GeS and 2.92 eV for 2D GeS₂. Thus, we have narrow, medium, and wide band gap 2D semiconductors: Ge₂S, GeS, and GeS₂, respectively.

Fig. 5 Electronic band structures using PBE (blue)/HSE06 (magenta) functionals and PBE total DOS with the corresponding projected density of states (PDOS) of Ge and S atoms are presented for 2D materials: (a) Ge₂S, (b) GeS, and (c) GeS₂. The zero of the energy axis is referenced to a horizontal red-dashed line that aligns with the valence band maximum (VBM). The band edges, i.e., VBM and CBM (conduction band minimum), are highlighted with solid green circles. The fundamental band gaps (in eV) are mentioned for each case.

3.4 Phonon transport in 2D Ge_xS_y compounds

Lattice thermal conductivity (κ_l) is a key characteristic of phonon transport and is obtained by iteratively solving the linearized phonon Boltzmann transport equation (PBTE). The solution of PBTE relies on second and third-order interatomic force constants (IFCs) as inputs. The influence of four-phonon scattering is disregarded due to its minimal impact on phonon transport in materials like 2D Ge_xS_y, which lack optical-acoustic gaps in phonon dispersions⁴² (see Fig. 2). The third-order IFCs capture the phonon–phonon scattering processes and play a crucial role in predicting reliable thermal conductivity. A cutoff radius determines the region of interatomic interactions for third-order IFC calculations. Selecting an appropriate cutoff radius ensures the appropriate inclusion of anharmonic interactions and improves the accuracy of κ_l. Furthermore, solving the PBTE involves discretizing the reciprocal space into a q-point grid hence achieving an accurate κ_l requires a sufficiently dense q-point grid. Considering inadequate cutoff radii and q-point grids can lead to orders of variations in the reported κ_l values. For instance, in the case of phosphorene, different κ_l values in the literature can be attributed to variations in the chosen cutoff radii and q-point grids.^39,40 Therefore, to avoid ambiguity in κ_l, we present its convergence against the cutoff radius and q-point grid in Fig. S5 of the ESI.† These figures depict an enormous variation of κ_l in the lower cutoff radius, which does not include sufficient nearest neighbors for interactions, and a coarse q-point grid. We achieve sufficient convergent κ_l by considering cutoff radii of 7.4 Å, 10.8 Å, and 6.3 Å for 2D Ge₂S, GeS, and GeS₂, respectively. These include up to the 20th nearest neighbors for Ge₂S and 22nd for GeS and GeS₂. Our q-point grid tests show that κ_l well converges at 60 × 60 × 1, 100 × 100 × 1, and 45 × 45 × 1 for Ge₂S, GeS, and GeS₂, respectively.

Fig. 6(a)–(c) depict the temperature-dependent lattice thermal conductivities of the 2D Ge₂S within the stable temperature range of 200–600 K, and the 2D GeS and GeS₂ compounds over the temperature range of 200–900 K. Specifically, GeS has a room temperature κ_l of 6.82 W m⁻¹ K⁻¹, which is lower than monolayer MoSe₂, WSe₂, and phosphorene.^43,44 Importantly, it is also smaller than its meta-stable orthorhombic phase (10.5/7.8 W m⁻¹ K⁻¹ along the zigzag/armchair directions) but higher than other IV–VI 2D compounds like GeSe (6.7/5.2 W m⁻¹ K⁻¹), SnS (4.7/4.4 W m⁻¹ K⁻¹), and SnSe (2.6/2.4 W m⁻¹ K⁻¹) along the respective directions.⁶ On the other hand, 2D Ge₂S exhibits a room temperature κ_l of ∼29 W m⁻¹ K⁻¹, which is the highest among the predicted 2D Ge_xS_y compounds and closer to that of phosphorene in the zigzag direction.⁴⁴ The novel stable phase of 2D GeS₂ has 16.95/10.68 W m⁻¹ K⁻¹ along the zigzag/armchair directions at 300 K, larger than the reported value for its meta-stable 1T phase.⁷ The anisotropy of κ_l can be quantified by the ratio κ_aa/κ_bb, where κ_aa and κ_bb represent the κ_l values along a- and b-axes, respectively. In the case of GeS₂, we designate the a-axis as the zigzag direction and the b-axis as the armchair direction (see Fig. 4(c)). For 2D GeS, this ratio is exactly one, indicating isotropic phonon transport. Similarly, for 2D Ge₂S, the ratio is nearly one, implying a predominantly isotropic nature of κ_l. Interestingly, the in-plane anisotropy ratio κ_aa/κ_bb is 1.59 for the predicted stable 2D GeS₂, which is higher than other orthorhombic 2D IV–VI compounds such as GeS (1.35), GeSe (1.29), SnS (1.07), and SnSe (1.08), indicating giant in-plane anisotropic lattice thermal conductivity. Although phosphorene exhibits a larger in-plane anisotropy ratio of 2.21⁴⁴ compared to 2D GeS₂, it is prone to stability issues under ambient temperatures.⁴⁵ These intrinsic anisotropic materials are rare and of particular importance in nanodevice applications. Moreover, Fig. 6(a)–(c) show a consistent decrease in κ_l as the temperature increases, displaying a strong adherence to the 1/T relationship, as indicated by the solid fitting line in each plot. This trend is commonly observed in intrinsic 2D materials and is attributed to the dominant occurrence of three-phonon scattering via the umklapp process, which depends on temperature directly. The umklapp scattering process opposes heat flow by scattering phonons in the opposite direction, impeding their effective transport. In contrast, the phonon–phonon scattering by the normal process could not provide resistance to phonon transport.⁴⁶

Fig. 6 The lattice thermal conductivity as a function of temperature and the corresponding 1/T fitting for 2D materials (a) Ge₂S, (b) GeS, and (c) GeS₂. The cumulative κ_l is plotted against the phonon mean free path (MFP) and fitted using the expression for (d) Ge₂S, (e) GeS, and (f) GeS₂, where the value (in nm) of the fitting parameter Λ₀ along the a/b-axis is mentioned in each case.

Nanostructuring is a common approach to control the κ_l of a material. Generally, decreasing the grain size causes a reduction in the mean free path (MFP) Λ of phonons, resulting in a lower κ_l while keeping the electronic conductivity unchanged. This behavior offers significant potential for enhancing the system's thermoelectric figure of merit. The cumulative κ_l as a function of MFP in Fig. 6(d)–(f) provides insights into the effect of grain size on the κ_l of 2D Ge_xS_y compounds at 300 K. To investigate the MFP range that dominantly contributes to κ_l, we fit the data with a logistic function , where κ_max is the ultimate lattice thermal conductivity, and Λ₀ is the fitting parameter called the representative MFP at which cumulative κ_l reaches 50% of its maximum value. The higher values of Λ₀, i.e., 569/371 nm for Ge₂S and 125/170 nm for GeS₂ compared to other IV–VI monolayers,⁴⁷ indicate that nanostructuring will be highly effective in reducing their κ_l. In contrast, the low phonon MPF Λ₀ of 16.5 nm for GeS, comparable to orthorhombic IV–VI monolayers, implies that it will be hard to reduce its κ_l with nanostructuring and other methods like alloying might be a convenient way. For a given nanostructure size of 50 nm, the κ_l is reduced to 5.70/7.45 W m⁻¹ K⁻¹ for Ge₂S and 5.54/2.75 W m⁻¹ K⁻¹ for GeS₂, which becomes competitive with κ_l (5.83 W m⁻¹ K⁻¹) of GeS. The evident anisotropy in the Λ₀ for Ge₂S and GeS₂ along the a/b-axes implies that reducing the grain size of the system will have different impacts on the κ_l along the a-axis compared to the b-axis. Similar anisotropic trends in Λ₀ for phosphorene allotropes are reported in previous studies.⁴⁸ This anisotropy offers intriguing possibilities for tailoring the κ_l by providing additional flexibility for design purposes. The analysis of the cumulative κ_l at various temperatures within the range of 300–600 K for 2D Ge₂S and 300–800 K (see Fig. S6 of the ESI†) for 2D GeS and GeS₂ demonstrates a uniform decrease in Λ₀ as the temperature rises.

Understanding the mechanisms involved in heat conduction through lattice vibrations in 2D Ge_xS_y compounds requires comprehensive analysis. We calculate the phonon group velocities using the formula , where λ represents the mode index, ω is the frequency, and q is the wave vector of the phonon, and the results are shown in Fig. 7. Fig. 7 shows the lower group velocities, i.e., below 6 km s⁻¹, in all 2D Ge_xS_y compared to other 2D materials like graphene, silicene,⁴⁹ MoS₂,⁵⁰ and phosphorene,⁴⁴ and the overall lower κ_l can be attributed to it. However, Fig. 7(b) and (e) depict comparable group velocities in the frequency range 200–300 cm⁻¹ of the optical phonon to its acoustic branches and is higher than the optical phonon of Ge₂S and GeS₂. This implies that the lower group velocity alone does not guarantee the lower κ_l, and further investigation is required, particularly in conjunction with phonon lifetimes. The phonon lifetimes, which represent the time a phonon can exist before scattering with other lattice vibrations, for 2D Ge_xS_y compounds, are shown in Fig. 8(a)–(c). Specifically, the short phonon lifetimes of less than 10 ps in 2D GeS indicate frequent scattering events. These scattering events impede the heat transfer through the lattice, resulting in the reduced κ_l of GeS compared to Ge₂S and GeS₂.

Fig. 7 Figures (a)–(c) shows the phonon group velocity (v_λ) along Γ–X for Ge₂S, GeS, and GeS₂, while (d) and (f) represent v_λ along Γ–Y for Ge₂S and GeS₂, and (e) along Γ–K for GeS, respectively.

Fig. 8 The phonon lifetimes of 2D Ge₂S, GeS, and GeS₂ are shown in (a), (b), and (c), respectively. The insets display the enlarged view for visual comparison in the frequency range of 50–150 cm⁻¹.

The anisotropic κ_l observed in 2D GeS₂ can be attributed to the anisotropic group velocities of both acoustic and low-frequency optical phonons in the frequency 50–150 cm⁻¹ range. Fig. 7(b) and (e) provide evidence of larger group velocities along the Γ–X reciprocal lattice direction compared to the Γ–Y direction in GeS₂ for acoustic and low-frequency optical phonons. This anisotropy in group velocities directly translates into anisotropic κ_l along the respective crystallographic directions. The frequency range of 50–150 cm⁻¹, mainly populated by optical phonons in GeS₂ and Ge₂S compositions, plays a crucial role due to the larger phonon lifetimes observed within this range. While heat transfer typically dominates in acoustic phonons, the acoustic-optical (low-frequency) phonon coupling effect allows for exceptions to this trend, as demonstrated in a recent study of nitride perovskite LaWN₃, where 60% of the heat transfer occured through optical phonons.⁵¹ In the case of GeS₂ and Ge₂S, the phonon lifetimes within the mentioned frequency range are approximately 20 ps, which is four times longer than the lifetimes observed in GeS (∼5 ps), evident from the insets showing enlarged views in Fig. 8(a)–(c). This highlights the significant role played by low-frequency optical phonons in heat transport. The frequency-resolved cumulative κ_l analysis (see Fig. S7(c) of the ESI†) further confirms the substantial contributions of low-frequency optical phonons to the overall κ_l. Moreover, the thermal anisotropy in GeS₂ aligns with its elastic properties. The zigzag direction exhibits a larger Young's modulus than the armchair direction, resulting in a larger κ_l along the zigzag direction.

Importantly, within the dominant heat-carrying frequency range (see Fig. S7 of the ESI†), the phonon lifetimes of 2D Ge₂S and GeS₂ are significantly larger than those of 2D GeS. However, this increase does not correspond to a proportional rise in κ_l. It highlights that phonon–phonon scattering cannot solely be determined by phonon lifetimes. Therefore, we investigate the lattice anharmonicity and the weighted phase space for three-phonon scattering. In Fig. 9, we display the Grüeisen parameter , where a₀ represents the equilibrium lattice constant, serving to quantify the lattice anharmonicity of materials. Additionally, the weighted phase space W of three-phonon scattering (as defined by Wu Li⁵²) for 2D Ge_xS_y is presented in the lower panel of Fig. 9. Both these parameters, i.e., γ and W, exhibit an inverse relationship with phonon lifetimes, providing additional insights. The large magnitude of γ and phase space W in the range of acoustic frequencies, which serve as the primary heat carriers, is evident from Fig. 9(a) and (d) and aligns with the lower κ_l of 2D GeS. On the other hand, a comparison of γ, W between 2D Ge₂S and 2D GeS₂ reveals that while the |γ|_max of Ge₂S (∼5) is greater than that of 2D GeS₂ (∼3), the weighted phase space for Ge₂S decreases significantly compared to GeS₂ as the phonon frequency increases within the frequency range of 150 cm⁻¹. Consequently, the weaker anharmonicity of 2D GeS₂ is compensated by the larger three-phonon weighted phase space, resulting in a smaller κ_l compared to 2D Ge₂S. The large W can be attributed to the higher phonon population in the lower-optical branches, as shown in Fig. 2(c), due to the greater number of atoms in the primitive unit cell of GeS₂, increasing the probability of phonon–phonon scattering.⁵³ Furthermore, the higher phase space [as observed in Fig. 9(e)] of GeS₂ in the frequency range of 150 cm⁻¹ supports our argument regarding its considerable anisotropic κ_l, attributed to both acoustic and low-frequency optical phonons within the 50–150 cm⁻¹ frequency range.

Fig. 9 Phonon-mode-resolved Grüneisen parameter and the room temperature weighted phase space (W) for (a and d) 2D GeS, (b and e) 2D GeS₂, and (c and f) 2D Ge₂S. The gray dashed lines (in c and f) aid in visually comparing the W within the 50–150 cm⁻¹ frequency range.

4 Conclusion

In this manuscript, we present the discovery of novel energetically stable 2D configurations for Ge₂S, GeS, and GeS₂ structures, achieved through the utilization of the USPEX in conjunction with VASP. Their dynamic, elastic, and thermal stabilities are confirmed using state-of-the-art methods. Subsequently, we characterize their electronic structures, revealing the semiconducting nature with narrow, medium, and wide band gaps for Ge₂S, GeS, and GeS₂, respectively. Then, we focus on their lattice thermal conductivity (κ_l), which provides insights into their heat transport properties. All structures exhibit overall low κ_l values, with 2D GeS having the lowest value of 6.82 W m⁻¹ K⁻¹. We observe a significant decrease in κ_l of 2D Ge₂S and GeS₂ with decreasing grain size, suggesting that nanostructuring could further suppress their κ_l. Through comprehensive analysis of the underlying mechanisms, we uncover the factors contributing to the observed low κ_l in 2D GeS and the in-plane anisotropy in 2D GeS₂. In 2D GeS, the combination of shorter phonon lifetimes, strong anharmonicity, higher weighted phase space, and smaller group velocities account for its low κ_l. The κ_l anisotropy of 2D GeS₂ is attributed to the anisotropic group velocities of both acoustic and low-frequency optical phonons. The low-frequency optical phonons in 2D Ge₂S and GeS₂ exhibited lifetimes four times longer than that of 2D GeS, resulting in increased κ_l. These findings contribute to the fundamental understanding of heat transport in 2D Ge_xS_y compounds, which have potential in thermal management and related fields.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

This work was supported by National Research Foundation of Korea (NRF) grants funded by the Ministry of Science and ICT (MSIT), government of South Korea (Grant No. 2021R1F1A1048555).

References

1. K. S. Novoselov, A. K. Geim, S. V. Morozov, D. Jiang, Y. Zhang, S. V. Dubonos, I. V. Grigorieva and A. A. Firsov, Science, 2004, 306, 666–669.
2. A. A. Balandin, S. Ghosh, W. Bao, I. Calizo, D. Teweldebrhan, F. Miao and C. N. Lau, Nano Lett., 2008, 8, 902–907.
3. J. Li, M. Chen, A. Samad, H. Dong, A. Ray, J. Zhang, X. Jiang, U. Schwingenschlögl, J. Domke, C. Chen, Y. Han, T. Fritz, R. S. Ruoff, B. Tian and X. Zhang, Nat. Mater., 2022, 21, 740–747.
4. Y. Anno, Y. Imakita, K. Takei, S. Akita and T. Arie, 2D Mater., 2017, 4, 025019.
5. M. Zulfiqar, Y. Zhao, G. Li, Z. Li and J. Ni, Sci. Rep., 2019, 9, 4571.
6. A. Shafique and Y.-H. Shin, Sci. Rep., 2017, 7, 506.
7. X. Wang, W. Feng, C. Shen, Z. Sun, H. Qi, M. Yang, Y. Liu, Y. Wu and X. Wu, Front. Mater., 2021, 8, 709757.
8. T. Jiang, Y.-J. Zhu, X.-J. Ye and C.-S. Liu, J. Appl. Phys., 2022, 132, 075001.
9. R. E. Abutbul, E. Segev, U. Argaman, G. Makov and Y. Golan, Adv. Mater., 2018, 30, 1706285.
10. W. H. Zachariasen, Phys. Rev., 1932, 40, 917–922.
11. M. Shimada and F. Dachille, Inorg. Chem., 1977, 16, 2094–2097.
12. G. Dittmar and H. Schäfer, Acta Crystallogr., Sect. B: Struct. Crystallogr. Cryst. Chem., 1975, 31, 2060–2064.
13. W. H. Zachariasen, J. Chem. Phys., 2004, 4, 618–619.
14. Z. V. Popović, M. Holtz, K. Reimann and K. Syassen, Phys. Status Solidi B, 1996, 198, 533–537.
15. Y. Li, G. Kuang, Z. Jiao, L. Yao and R. Duan, Mater. Res. Express, 2022, 9, 122001.
16. A. J. Mannix, X.-F. Zhou, B. Kiraly, J. D. Wood, D. Alducin, B. D. Myers, X. Liu, B. L. Fisher, U. Santiago, J. R. Guest, M. J. Yacaman, A. Ponce, A. R. Oganov, M. C. Hersam and N. P. Guisinger, Science, 2015, 350, 1513–1516.
17. S. Lee, J.-E. Jung, H.-G. Kim, Y. Lee, J. M. Park, J. Jang, S. Yoon, A. Ghosh, M. Kim, J. Kim, W. Na, J. Kim, H. J. Choi, H. Cheong and K. Kim, Nano Lett., 2021, 21, 4305–4313.
18. S. Cahangirov, M. Topsakal, E. Aktürk, H. Sahin and S. Ciraci, Phys. Rev. Lett., 2009, 102, 236804.
19. X. Yu, L. Li, X.-W. Xu and C.-C. Tang, J. Phys. Chem. C, 2012, 116, 20075–20079.
20. A. R. Oganov and C. W. Glass, J. Chem. Phys., 2006, 124, 244704.
21. A. R. Oganov, A. O. Lyakhov and M. Valle, Acc. Chem. Res., 2011, 44, 227–237.
22. A. Ali and Y.-H. Shin, Phys. Chem. Chem. Phys., 2022, 24, 29772–29780.
23. G. Kresse and J. Furthmüller, Comput. Mater. Sci., 1996, 6, 15–50.
24. G. Kresse and D. Joubert, Phys. Rev. B: Condens. Matter Mater. Phys., 1999, 59, 1758.
25. J. P. Perdew, K. Burke and M. Ernzerhof, Phys. Rev. Lett., 1996, 77, 3865–3868.
26. S. Grimme, J. Antony, S. Ehrlich and H. Krieg, J. Chem. Phys., 2010, 132, 154104.
27. J. Heyd, G. E. Scuseria and M. Ernzerhof, J. Chem. Phys., 2003, 118, 8207–8215.
28. G. J. Martyna, M. L. Klein and M. Tuckerman, J. Chem. Phys., 1998, 97, 2635.
29. A. Togo and I. Tanaka, Scr. Mater., 2015, 108, 1–5.
30. W. Li, J. Carrete, N. A. Katcho and N. Mingo, Comput. Phys. Commun., 2014, 185, 1747–1758.
31. W. Li, L. Lindsay, D. A. Broido, D. A. Stewart and N. Mingo, Phys. Rev. B: Condens. Matter Mater. Phys., 2012, 86, 174307.
32. X. Gonze and C. Lee, Phys. Rev. B: Condens. Matter Mater. Phys., 1997, 55, 10355–10368.
33. B. D. Malone and E. Kaxiras, Phys. Rev. B: Condens. Matter Mater. Phys., 2013, 87, 245312.
34. B. Dong, Z. Wang, N. T. Hung, A. R. Oganov, T. Yang, R. Saito and Z. Zhang, Phys. Rev. Mater., 2019, 3, 013405.
35. A. K. Singh and R. G. Hennig, Appl. Phys. Lett., 2014, 105, 042103.
36. M. Maździarz, 2D Mater., 2019, 6, 048001.
37. X. Li, L. Tao, Z. Chen, H. Fang, X. Li, X. Wang, J.-B. Xu and H. Zhu, Appl. Phys. Rev., 2017, 4, 021306.
38. Z. Gao, X. Dong, N. Li and J. Ren, Nano Lett., 2017, 17, 772–777.
39. L. Zhu, G. Zhang and B. Li, Phys. Rev. B: Condens. Matter Mater. Phys., 2014, 90, 214302.
40. G. Qin, X. Zhang, S.-Y. Yue, Z. Qin, H. Wang, Y. Han and M. Hu, Phys. Rev. B, 2016, 94, 165445.
41. D. R. Clarke, Surf. Coat. Technol., 2003, 163–164, 67–74.
42. T. Feng, L. Lindsay and X. Ruan, Phys. Rev. B, 2017, 96, 161201.
43. S. Kumar and U. Schwingenschlögl, Chem. Mater., 2015, 27, 1278–1284.
44. G. Qin, Q.-B. Yan, Z. Qin, S.-Y. Yue, M. Hu and G. Su, Phys. Chem. Chem. Phys., 2015, 17, 4854–4858.
45. J. D. Wood, S. A. Wells, D. Jariwala, K.-S. Chen, E. Cho, V. K. Sangwan, X. Liu, L. J. Lauhon, T. J. Marks and M. C. Hersam, Nano Lett., 2014, 14, 6964–6970.
46. C. Kittel, Introduction to Solid State Physics, Wiley, 2004.
47. G. Qin, Z. Qin, W.-Z. Fang, L.-C. Zhang, S.-Y. Yue, Q.-B. Yan, M. Hu and G. Su, Nanoscale, 2016, 8, 11306–11319.
48. J. Zhang, H. J. Liu, L. Cheng, J. Wei, J. H. Liang, D. D. Fan, P. H. Jiang and J. Shi, Sci. Rep., 2017, 7, 4623.
49. B. Peng, H. Zhang, H. Shao, Y. Xu, G. Ni, R. Zhang and H. Zhu, Phys. Rev. B, 2016, 94, 245420.
50. B. Peng, H. Zhang, H. Shao, Y. Xu, X. Zhang and H. Zhu, Ann. Phys., 2016, 528, 504–511.
51. Q. Ren, Y. Li, Y. Lun, G. Tang and J. Hong, Phys. Rev. B, 2023, 107, 125206.
52. W. Li and N. Mingo, Phys. Rev. B: Condens. Matter Mater. Phys., 2015, 91, 144304.
53. C. Liu, C. Wu, T. Song, Y. Zhao, J. Yang, P. Lu, G. Zhang and Y. Chen, ACS Appl. Energy Mater., 2022, 5, 15356–15364.

---

Footnote

† Electronic supplementary information (ESI) available. See DOI: https://doi.org/10.1039/d3cp04568d

---