Lattice vibrations and electronic properties of GaSe nanosheets from first principles

TL;DR

This study uses first-principles calculations to explore the electronic and vibrational properties of GaSe nanosheets of varying thickness, revealing their semiconducting nature, stability, and vibrational spectra characteristics.

Contribution

It provides a detailed first-principles analysis of the electronic band structure and lattice dynamics of GaSe nanosheets across different layer counts, highlighting their stability and vibrational features.

Findings

GaSe nanosheets are semiconductors with an indirect band gap.

All systems studied are dynamically stable.

Distinct Raman and IR spectral features are identified.

Abstract

Electronic properties and lattice dynamics of bulk $\varepsilon$-GaSe and mono-, bi- and tri-tetralayer GaSe are investigated by means of density functional and density functional perturbation theory. The few-tetralayers systems are semiconductors with an indirect nature of the fundamental band gap and a Mexican-hat-shape is observed at the top of the valence band. The phonon branches analysis reveals the dynamical stability for all systems considered together with the LO-TO splitting breakdown in two-dimensional systems. In-plane (E) and out-of-plane (A) zone-center lattice vibrations dominate the Raman and IR spectra.

Figures

Tables

Table 1. Table 1: Parameters of atomic geometry (in Å) such as lattice constants a 𝑎 a , c 𝑐 c and atomic distances. Electronic band gaps (in eV) are also listed.

xc-functional	$a$	$c$	Se-Se	Ga-Ga	Ga-Se	Se-Se	Indirect	Direct	Fundamental
			(intra)	(intra)	(inter)	(inter)	band gap	band gap	gap (in HSE)
bulk
LDA-D2	3.651	15.336	4.727	2.370	2.415	3.618		0.97
LDA	3.72	15.56	4.73	2.41	2.44	3.73		0.97
GGA-D2	3.662	15.587	4.812	2.432	2.470	3.826		0.92	1.75
Exp. Kuhn et al. (1975); Li and Li (2015)	3.735	15.887	4.766	2.383	2.485	–		2.12
trilayer
LDA-D2	3.704		4.749	2.409	2.437	3.804	1.30	1.34
LDA	3.71		4.74	2.41	2.44	3.82	1.30	1.34
GGA-D2	3.753		4.834	2.446	2.472	3.835	1.31	1.35	2.25
bilayer
LDA-D2	3.704		4.742	2.405	2.437	3.812	1.54	1.63
LDA	3.71		4.73	2.41	2.44	3.75	1.54	1.63
GGA-D2	3.748		4.818	2.428	2.472	3.845	1.60	1.67	2.58
monolayer
LDA	3.73		4.72	2.41			2.16	2.30
GGA-D2	3.745		4.821	2.429			2.20	2.32	3.24

Table 2. Table 2: Calculated optical modes of GaSe nanosheets from 1TL to 3TL and bulk ε 𝜀 \varepsilon -GaSe at zone center characterized by their corresponding irreducible representation. Frequencies are in cm -1 . Raman active modes are labeled with ”R”, IR active modes with ”I”, and experimentally not observed modes with ”N.O.”. E 𝐸 E modes are twofold degenerate, while A 𝐴 A modes are not. Experimental values are extracted from Ref.  Lei  et al. ( 2013 ) and a  Ref. Rahaman2018

Theory

Exp.

1 TL

$\omega$

2 TL

$\omega$

3 TL

$\omega$

Bulk

$\omega$

1-2TL

few TL

Bulk

$E$

22.4

I+R

$E^{{}^{\prime \prime }}$

24.2

R

$E^{{}^{\prime }}$

13.6

I+R

$E^{{}^{\prime }}$

25.5

I+R

$A_1$

26.3

I+R

$A_1^{{}^{\prime }}$

18.4

R

$A_2^{{}^{\prime \prime }}$

26.8

I

$A_2^{{}^{\prime \prime }}$

31.1

I

$E^{{}^{\prime \prime }}$

55.2

R

$E$

57.0

I+R

$E^{{}^{\prime \prime }}$

58.1

R

$E^{{}^{\prime \prime }}$

54.7

R

59

59

59

$E$

58.2

I+R

$E^{{}^{\prime }}$

59.2

I+R

$E^{{}^{\prime \prime }}$

55.7

R

$E^{{}^{\prime \prime }}$

59.5

R

$A_1^{{}^{\prime }}$

128.8

R

$A_1$

128.5

I+R

$A_1^{{}^{\prime }}$

128.3

R

$A_1^{{}^{\prime }}$

127.2

R

132

132

132

$A_1$

131.8

I+R

$A_2^{{}^{\prime \prime }}$

130.1

I

$A_1^{{}^{\prime }}$

129.8

R

$A_1^{{}^{\prime }}$

132.3

R

$E^{{}^{\prime \prime }}$

207.7

R

$E$

208.0

I+R

$E^{{}^{\prime \prime }}$

207.8

R

$E^{{}^{\prime \prime }}$

205.9

R

N.O.

N.O.

208

$E$

209.0

I+R

$E^{{}^{\prime }}$

209.2

I+R

$E^{{}^{\prime \prime }}$

207.1

R

$E^{{}^{\prime \prime }}$

209.2

R

$E^{{}^{\prime }}$

213.5

I+R

$E$

214.0

I+R

$E^{{}^{\prime }}$

213.3

I+R

$E^{{}^{\prime }}$

211.5

I+R

211a

$E$

215.1

I+R

$E^{{}^{\prime \prime }}$

215.4

R

$E^{{}^{\prime }}$

212.6

I+R

$E^{{}^{\prime }}$

215.4

I+R

$A_2^{{}^{\prime \prime }}$

245.2

I

$A_1$

243.9

I+R

$A_2^{{}^{\prime \prime }}$

242.7

I

$A_2^{{}^{\prime \prime }}$

242.9

I

$A_1$

244.5

I+R

$A_2^{{}^{\prime \prime }}$

244.3

I

$A_2^{{}^{\prime \prime }}$

243.6

I

$A_1^{{}^{\prime }}$

244.4

R

$A_1^{{}^{\prime }}$

309.7

R

$A_1$

308.9

I+R

$A_1^{{}^{\prime }}$

307.6

R

$A_1^{{}^{\prime }}$

306.3

R

303.4

303.4

305.2

$A_1$

309.3

I+R

$A_1^{{}^{\prime }}$

309.0

R

$A_1^{{}^{\prime }}$

306.8

R

306a

$A_2^{{}^{\prime \prime }}$

309.0

I

Full text

††preprint: APS/123-QED

I Introduction

The discovery of graphene Geim.Novoselov:2007:NM with its unique physical and chemical properties has triggered the research on other two-dimensional (2D) inorganic materials, like transition metal dichalcogenides and transition metal oxides Miro et al. (2014); Niu and Li (2015). Several three-dimensional (3D) van der Waals (vdW) solids have been shown to be exfoliable down to atomic monolayers (MLs), giving rise to fabrication of several 2D crystals Koski and Cui (2013). Stable single layer or few layers of 2D compounds can be produced via chemical exfoliation or mechanical cleavage starting from the corresponding III-VI layered crystal counterparts Ashton et al. (2017); Bandurin.Tyurnina.ea:2016:NN ; Hu et al. (2012). In parallel, the interesting research field of heterostructures fabrication by stacking different 2D crystals on top of each other has flourished Geim.Grigorieva:2013:N ; Cannuccia2012 ; Prete2017 ; Prete2018 .

Among the transition-metal monochalcogenides MX with M = Ga, In and X = S, Se, Te, GaSe is a prototypical example. It consists of covalently bonded Ga and Se atoms, vertically stacked in four-atom Se-Ga-Ga-Se layers (here on named tetralayers, TL) held together by vdW forces. The single tetralayer thickness is 0.93 nm Hu et al. (2012). There are four possible stacking arrangements for bulk GaSe, designated as $\beta$, $\varepsilon$, $\gamma$, and $\delta$ Kuhn et al. (1975). The $\beta$, $\varepsilon$ polytypes have two-tetralayer hexagonal stacking sequences (2H), i.e., eight atoms, in the unit cell, and they crystallize respectively into $P_{3}/mmc$ and $P\overline{6}m2$ space groups. The $\delta$ polytype is 4H, while the $\gamma$ polytype crystallizes in a rhombic or trigonal lattice (3R). With a band gap of about 2 eV Le Toullec et al. (1975) and a strong optical oscillator strength, bulk GaSe polytypes have potential optoelectronic applications in the visible spectral range. Electronic confinement effects in GaSe nanosheets may shift the optical absorption edge toward higher energies and increase the oscillator strength Lei et al. (2013).

The $\varepsilon$-GaSe is among the polytypes the most common and the most studied one. First-principles calculations Rybkovskiy et al. (2011); Ma et al. (2013) predict the valence band maximum (VBM) to be situated at $\Gamma$ point, and the conduction band minimum (CBM) close to point M of the hexagonal Brillouin zone (BZ) whereas the direct gap is only slightly different from the indirect one. Similar investigations have been performed for GaSe layered structures with respect to the slab thickness Rybkovskiy et al. (2011); Ma et al. (2013); Zólyomi et al. (2013); Zhuang and Hennig (2013); Li and Li (2015) with the inclusion in some cases of the quasi-particle effects within Hedin’s GW approximation Rybkovskiy et al. (2011); Zhuang and Hennig (2013). The increase of the fundamental gap with a decreasing thickness of the slab has been demonstrated. It is accompanied by the formation of a Mexican-hat or camel-back dispersion of the uppermost valence band near the $\Gamma$ point, making the few-tetralayer systems indirect semiconductors.

In contrast to the electronic structure, the phonons of GaSe are less investigated. Since Raman spectroscopy is employed to determine the thickness (and thus, the number of unit layers) of ultrathin flakes in a non-destructive, and unambiguous manner, Raman studies are available now for bulk and a-few-tetralayer systems Lei et al. (2013); Ibragimov1989 ; Late et al. (2012); Zhou et al. (2014); Rodriguez et al. (2014); Rahaman2018 . Theoretical investigations of phonon modes have been only performed for a tetralayer of GaSe Zólyomi et al. (2013); Pandey et al. (2017). For that reason, in this paper we investigate the vibrational properties of $\varepsilon$-GaSe and layered GaSe by varying the slab thickness from bulk, 3 TL down to 1 TL. The methods are described in Sec. II. In section III the results on electronic and dynamical properties are discussed. Finally, a summary and conclusions are given in Sec. IV.

II Computational methods

The ab initio calculations are based on Density Functional Theory (DFT) using different exchange-correlation (xc) functionals as implemented in the Quantum Espresso (QE) package qua . The electron-ion interaction is described by pseudopotentials for Ga and Se cores. Thereby the semicore Ga 3$d$ electrons are treated as valence electrons. Both the local density approximation (LDA) Kohn and Sham (1965) and the generalized gradient approximation (GGA) as proposed by Perdew, Burke, and Ernzerhof (PBE) Perdew et al. (1996) are applied. The influence of the vdW interaction between the covalently bonded Se-Ga-Ga-Se tetralayers is studied by adding a semiempirical dispersion (D) potential to the DFT total energy, through a pairwise force field according to Grimme’s DFT-D2 method Grimme et al. (2010).

For all considered GaSe structures the calculations of electronic properties are performed using periodic boundary conditions and plane-wave expansion of the electronic eigenfunctions. The bulk unit cell is represented in Fig. 1 with its side and top view. The nanosheets with two and three tetralayers are stacked using the single tetralayer (see Fig. 1 (b)) as building block and respecting the $z$-stacking of the bulk Bechstedt (2003). A vacuum spacing of 28 Å between the periodic images of mono-, bi-, and tri-tetralayers is chosen in order to diminish the electronic interaction between the nanosheets. To sample the Brillouin Zone (BZ) in the total energy and self-consistent electronic structure calculations, a Monkhorst-Pack (MP) special point mesh Monkhorst and Pack (1976) of 12$\times$12$\times$1 (6$\times$6$\times$1, 6$\times$6$\times$1, 6$\times$6$\times$6) is applied in the case of a GaSe 1TL, (2TL, 3TL, bulk $\varepsilon$ polytype). The structural optimizations are performed by relaxing the atomic geometries till forces are less than $10^{-6}$ a.u. The plane-wave energy cutoff and the total energy threshold are set to 300 Ry and $10^{-8}$ Ry, respectively. All DFT approaches lead to a significant underestimation of fundamental energy gaps Bechstedt (2015). Since GW self-energy calculations would be too expensive, we investigate the band structures of the GaSe geometries within an approximated quasiparticle approach. Explicitly, we apply the hybrid xc-functional HSE06 Heyd et al. (2003).

Dynamical properties, i.e., the phonon dispersions as well as Raman frequencies, are calculated within Density Functional Perturbation Theory (DFPT) Baroni et al. (2001) as implemented in the QE package. In particular, we use the qe-6.3 version because it implements a truncated Coulomb interaction in the direction parallel to the layered-system stacking orientation ($z$-direction for structures periodic in the $x$-$y$ plane). The implementation of the cutoff guarantees the correct treatment of the LO-TO splitting breakdown in polar two-dimensional materials Sohier-2Dcutoff .

As a compromise for a reasonable simulation of strong intrasheet and weak intersheet interactions we also apply here the LDA xc-functional Kohn and Sham (1965). A plane-wave cutoff of 90 Ry and robust convergence criteria (total energy/forces for the ionic minimization and self-consistence criteria) are necessary to get the correct interatomic force constants to be used in the phonons calculations. A mesh of 21$\times$21$\times$1 k-points is used for the BZ integration and sampling of the reciprocal space for 1TL-GaSe. For 2TL and 3TL cases, given the increasing number of atoms and the volume of the simulation cell, we adopted a less dense q-grid (12$\times$12$\times$1) to keep the calculations reasonable in terms of cpu time. The length of the unit cell along the $z$ direction was chosen to be larger than twice the thickness of the slab (including electrons) Sohier-2Dcutoff . For bulk $\varepsilon$-GaSe we used a mesh of 18$\times$18$\times$6 k-points and 9$\times$9$\times$3 mesh for the sampling of the reciprocal space.

III Results and Discussion

III.1 Atomic geometries

The results of the atomic geometry optimization of bulk GaSe and few-tetralayer systems are summarized in Table 1. Bulk $\varepsilon$-GaSe (hexagonal space group $P\bar{6}m2$ $(D^{1}_{3h})$ with point group $\bar{6}m2$ $(D_{3h})$) contains in the unit cell two Se-Ga-Ga-Se layers rotated against each other (Fig. 1a). The nanosheets with one and three tetralayers crystallize to the same space group as the bulk $\varepsilon$-GaSe, while the 2TL exhibits a reduced symmetry. It belongs to the trigonal space group $P3m1(C^{1}_{3v})$ with the point group $3m$ $(C_{3v})$.

The $a$ and $c$ lattice constants of $\varepsilon$-GaSe are underestimated with respect to the experimental values independently of the xc-functional employed. Surprisingly, the LDA bulk lattice constants are on the average the closest ones to the experimental data. Differences are less than 0.01 Å. For what concerns the dependence of in-plane lattice constants with the layer thickness we observe that with the vdW functionals, the $a$ value tends to increase. An homogeneity of the $a$ values with respect to the xc-functional is found in few-tetralayer systems. Going from LDA-D2 to LDA and finally to GGA-D2 the intra- and interatomic distances exhibit the average trend to increase. However, for a given xc-functional the trend with the number of layers is not unique. For LDA a decrease of both intralayer and interatomic distances with the decrease of the system thickness is observed. With vdW interaction the intratomic distances decrease whereas the Se-Se interlayer distances increase with decreasing material thickness. Independently of the xc-functional the Ga-Se interlayer distances remain the same. The found atomic geometries are reasonably well described also in LDA, in agreement with the literature Rybkovskiy et al. (2011); Zólyomi et al. (2013).

III.2 Electronic structures

In order to illustrate the electronic structure of the GaSe systems, we display their band structures in Fig. 2 for two different xc-functionals, GGA-PBE and LDA. Despite the neglect of quasiparticle corrections, all the studied systems represent semiconductors, whose gaps increase from bulk to a single tetralayer by more than 50 % (see Table 1) because of the electronic confinement effects. Within the HSE06 functional the gap gets increased by about 0.9 eV (also see Table 1), independent of the confinement. Comparing with the measured gap of $\varepsilon$-GaSe a gap increase of about 1.15 eV should be expected instead. A quasiparticle gap at $\Gamma$ of about 2.34 eV for bulk GaSe, and of $3.89$ eV for a GaSe 1TL, are predicted within GW calculations Rybkovskiy et al. (2011). The optical gaps, however, may be significantly smaller especially in the 1TL case, due to the large binding energies of the excitons in 2D systemsPulci et al. (2015, 2012). This is in agreement with the value of 2.58 eV experimentally obtained from optical-absorption spectroscopy Rybkovskiy et al. (2011). A larger optical gap of 3.3 eV as derived from the photocurrent spectrum has been explained by optical transitions from the second valence band to the conduction band minimum Lei et al. (2013).

Bulk $\varepsilon$-GaSe comes out in our calculations to be a direct gap semiconductor with the lowest optical transitions at the $\Gamma$ point between parabolic conduction and valence bands (see Fig. 2). The top of the valence band at $\Gamma$ possesses different masses along and perpendicular to the hexagonal axis. All the few-tetralayer systems are indirect semiconductors with gap values in agreement with previous calculations Rybkovskiy et al. (2011); Ma et al. (2013); Zólyomi et al. (2013). The uppermost valence band forms a Mexican-hat or camel-back dispersion near the $\Gamma$ point. As a consequence of the strong Rashba contribution to the spin-orbit interaction the indirect $\Gamma$K $\rightarrow$ $\Gamma$ gap is smaller than the direct $\Gamma$ $\rightarrow$ $\Gamma$ gap (see Table 1).

III.3 Phonon dispersions

The dispersion of the phonon branches, obtained by means of the DFPT and the LDA xc-functional for LDA-optimized atomic equilibrium geometries, are plotted along the wave-vector path connecting $\Gamma$-$K$-$M$-$\Gamma$ high symmetry points in Fig. 3, from single tetralayer to three tetralayers, and finally the bulk.

The phonon modes of bulk $\varepsilon$-GaSe can be classified by means of the symmetry of the zone-center modes. The 24 $\Gamma$-modes can be decomposed into four groups $(A^{{}^{\prime}}_{1}\oplus A^{{}^{\prime\prime}}_{2}\oplus E^{{}^{\prime}}\oplus E^{{}^{\prime\prime}})$. Among these, $A^{{}^{\prime\prime}}_{2}\oplus E^{{}^{\prime}}$ are acoustic and $4A^{{}^{\prime}}_{1}\oplus 3A^{{}^{\prime\prime}}_{2}\oplus 3E^{{}^{\prime}}\oplus 4E^{{}^{\prime\prime}}$ are optical modes. Except for $A^{{}^{\prime\prime}}_{2}$ modes which are only infrared active, all the others are Raman active modes. Out of them, $E^{{}^{\prime}}$ modes around 212 cm*-1* are also infrared active. The corresponding displacement patterns are displayed in Fig. 4.

1TL and 3TL have the same space group as the bulk, with 12 symmetry operations, while 2TL has a $C_{3v}$ $(3m)$ symmetry with six symmetry operations. Since it is not centrosymmetric, all optical modes are simultaneously Raman- and IR-active Allakhverdiev et al. (2006).

The bulk phonon dispersion has three acoustic modes $A^{{}^{\prime\prime}}_{2}\oplus E^{{}^{\prime}}$. The $E^{{}^{\prime}}$ mode vibrating in-plane (longitudinal and transverse acoustic) is two-fold degenerate and has a linear dispersion for vanishing wave vector. The $A^{{}^{\prime\prime}}_{2}$ mode is a flexural phonon mode (out-of-plane acoustic mode) which exhibits a peculiar quadratic dispersion near the zone center, typical of layered crystals Zabel2001 and which can be explained as a consequence of the point group symmetry Saito1998 . It is the easiest mode to be excited, because of its lowest frequency among all phonon modes and it may be then affected by numerical instabilities Zolyomi et al. (2014). Here the three acoustic branches converge smoothly to zero frequency for vanishing wave vector, demonstrating that numerical instabilities have been safely controlled: no soft modes or negative energies appear in the phonon branches of bulk and GaSe nanosheets. We conclude that not only the bulk $\varepsilon$-GaSe but also the layered systems are dynamically stable, at least, considering small 2D hexagonal unit cells.

Table 2 shows the bulk modes at $\Gamma$ together with their symmetry character. Modes of 1TL to 3TL are also listed, establishing the correspondence of few-tetralayer modes with bulk phonon symmetry representations. For each representation the number of phonon modes reflects the number of tetralayers. From 1TL to 3TL the phonon frequencies only vary of few cm*-1* because of the weak bonding between tetralayers.

GaSe is a polar material, hence the behavior of the longitudinal optical (LO) modes depends on dimensionality. The dipoles generated by the LO displacements interact with each other through long-range Coulomb interaction. This additional dipole-dipole interacting term causes the splitting between the LO and transverse optical (TO) phonon modes. In bulk this splitting is independent of the phonon wavevector $\vec{q}$, in particular in the limit of $\vec{q}\rightarrow 0$, we find the LO mode at 242.6 cm*-1* and the TO mode at 212 cm*-1*. In the isolated single tetralayer this splitting vanishes for $\vec{q}\rightarrow 0$ and in this limit the derivative of the LO phonon dispersion is discontinuous. Here the breakdown of the LO-TO splitting occurs at 213 cm*-1*. Increasing the number of tetralayers the slope of the highest LO mode becomes large and the LO frequency smoothly approaches the bulk limit lifting the degeneracy Sohier2017NL .

Turning to the analysis of 1TL phonon band structures, two linear branches, related to in-plane vibrations, a transversal acoustic (TA) and a longitudinal acoustic (LA) one, can be seen in Fig. 3(a). In addition, a flexural acoustic (ZA) branch is visible for out-of-plane vibrations. The two linear branches have large sound velocities of $\sim$2600 m/s (TA) and $\sim$4300 m/s (LA) in agreement with the experimental values Abutalybov1995 . Approaching the $\Gamma$ point the ZA branch has the quadratic dispersion already discussed for the bulk case Carrete et al. (2016). Immediately above the acoustic branches low-frequency $E^{\prime}$ optical modes are present in the bulk (13.6 cm*-1*), 2TL (22.4 cm*-1*) and 3TL (25.5 cm*-1*). These modes are interlayer shear modes, typical of layered materials, corresponding to the sliding of two adjacent tetralayers (see Fig. 4 upper panel). No shear modes are obviously present in the 1TL case.

III.4 Raman and IR spectra

In the lower frequency range of bulk phonon band structure (see Fig. 3(d) and Table 2) at the zone-center of the hexagonal BZ, a bunch of optical modes are found at 13.6 cm*-1* ($E^{{}^{\prime}}$, the shear mode), 54.7 and 55.7 cm*-1* ($E^{{}^{\prime\prime}}$) as in-plane modes and 26.8 cm*-1* ($A^{{}^{\prime\prime}}_{2}$) as an out-of-plane mode. An in-plane mode in the lower frequency region is experimentally found around 59 cm*-1* (E”) in Raman spectrum Lei et al. (2013). In the intermediate frequency range we find two Raman active $A^{{}^{\prime}}_{1}$ modes at 127.2 and 129.8 cm*-1* which correspond to a $A_{1}^{\prime}$ dominant peak in the measured Raman spectra at 132 cm*-1* Lei et al. (2013); Rahaman2018 - 134.4 cm*-1* Rodriguez et al. (2014). In the high energy region we find eight modes, two of them are $A_{1}^{\prime}$ modes at 306.3 and 306.8 cm*-1*, two Raman forbidden modes at 242.9 and 243.6 cm*-1* (LO modes) and 2 $E^{{}^{\prime\prime}}$, 2 $E^{{}^{\prime}}$ modes around 206 and 212 cm*-1* (the latter is a TO mode) respectively. The consequent LO-TO splitting amounts to 30 cm*-1*. Raman measurements reveal an intense peak at 305.2 Lei et al. (2013) - 306 Rahaman2018 - 307.8 cm*-1* Rodriguez et al. (2014) ($A_{1}^{\prime}$) and an $E^{{}^{\prime}}$ mode at 211 cm*-1* Rahaman2018 - 213 cm*-1* Rodriguez et al. (2014). A peak for $E^{{}^{\prime\prime}}$ mode around 208 cm*-1* appears in experimental Raman spectrum Lei et al. (2013) just for the bulk sample, while a peak at 230 cm*-1* is measured in single or double tetralayer samples accompained by a shift of about 3 cm*-1* to lower frequencies in thicker, few-tetralayer samples. The fact that we do not find this mode at such a frequency in few-tetralayer samples supports the hypothesis Lei et al. (2013) that this mode is due to the interaction of few-tetralayer samples with the substrate, i.e. it does not exist in the bulk.

Calculated one-phonon Raman and IR spectra are displayed in Fig. 5 for bulk GaSe and varying the nanosheet thickness. Only zone-center phonons contribute to the spectra. In the 1TL case only the two $A_{1}^{\prime}$ modes are clearly visible. As expected from group symmetry analysis, the twofold-degenerate $E^{\prime}$ mode at about 213 cm*-1* shows a finite IR strength.

Frequencies of Raman and IR modes at zone-center of the same symmetry are grouped in blocks in Table 2. The average of the frequencies weighted by the Raman intensities leads to the position of the Raman peaks in Fig. 5(a). Independently of the number of tetralayers two strong Raman peaks of $A_{1}/A_{1}^{\prime}$ symmetry appear in the range of $306-309$ cm*-1* and $127-132$ cm*-1*. The highest frequency mode is hardly layer-thickness-dependent, while the lowest one slightly decreases with decreasing layer thickness. Both trends are confirmed by recent Raman measurements Rahaman2018 . For what concerns the Raman intensities, in both cases the intensity increases with increasing layer thickness. In particular, a really intense peak is visible in the bulk Raman spectrum around 130 cm*-1*. The Raman peak position around $207-209$ cm*-1* related to $E^{\prime\prime}$ modes does hardly vary with the material thickness, in agreement with the fact that they are in-plane modes. The $A^{{}^{\prime}}_{1}$ mode as well as the IR-active $A^{{}^{\prime\prime}}_{2}$ mode are out-of-plane vibrations and, therefore, more influenced by the interlayer interaction. Thereby, the highest $A^{{}^{\prime\prime}}_{2}$ mode frequency increases with reduced slab thickness, while its IR intensity decreases. Whereas the absolute peak positions, in particular for bulk, agree well with measurements, e.g. in Ref. Rahaman2018 with deviations of about 1 cm*-1*, and the thickness variation of the lowest $A^{{}^{\prime}}_{1}$ mode are in agreement with the experimental findings, we report some discrepancies for the highest $A^{{}^{\prime\prime}}_{2}$ mode. We trace back the small differences between measurements Lei et al. (2013); Rodriguez et al. (2014); Rahaman2018 and our calculations to the presence of a substrate in experimental setups, whereas only freestanding nanosheets are theoretically investigated.

The IR spectra in Fig. 5b are less complex compared to the Raman ones in Fig. 5a. Only a strong in-plane $E^{\prime}$ peak is visible in the range of $213-215$ cm*-1*. In the bulk case this peak is slightly shifted toward $212$ cm*-1* close to the peak found in experimental GaSe spectra Giorgianni et al. (1977). Interestingly, we also predict a peak related to the out-of-plane $A_{2}^{\prime\prime}$ mode close to 243 cm*-1*. However, the observability of the two peaks have to be carefully discussed in terms of ordinary or extraordinary light polarization.

IV Summary and Conclusions

We have performed a theoretical study of 1TL, 2TL and 3TL GaSe nanosheets by means of density functional theory. For the purpose of comparison the same computations have been done for the bulk $\varepsilon$-GaSe crystal polytype. Three different exchange-correlation functionals have been applied in the ground-state calculations. We found the typical underestimation (overestimation) of the chemical bonding using GGA-(LDA-)based functionals. We demonstrated that the LDA functional delivers the lattice constants in a sufficiently good quality. It can therefore be applied also for DFPT studies of the lattice vibrations. The investigations of the electronic structures clearly demonstrate that the GaSe nanosheets are indirect semiconductors because of a Mexican-hat dispersion of the uppermost valence band, while the bulk system is a direct semiconductor. Quasiparticle effects, in general, and confinement effects in the nanosheets increase the fundamental gaps toward values, which are comparable with experimental ones.

The lattice-vibrational properties are dominated by twofold degenerate in-plane $E$-type modes and non-degenerate out-of-plane $A$-type modes. The phonon band dispersion versus high-symmetry lines in the 2D hexagonal BZ shows that the few-tetralayer nanosheets are also dynamically stable and that the phonon band dispersion of 1TL, 2TL and 3TL exhibits the LO-TO splitting breakdown typical of polar two-dimensional systems. The frequencies, Raman intensities, and IR oscillator strengths allow the construction of Raman and IR spectra in qualitative but also quantitative agreement with available experimental data. This holds for the most intense peaks, in general, but also for the thickness dependence of their peak positions.

Acknowledgements.

E.C. acknowledges support by the Programma per Giovani Ricercatori - 2014 “Rita Levi Montalcini”. O.P. acknowledges financial support from the EU MSCA HORIZON2020 project ’CoExAN’ (GA644076). CPU time was granted by CINECA HPC center. We thank T. Sohier and C. Attaccalite for useful discussions, A. Zappelli for the management of the computer cluster Rosa.