Spectroscopic properties of few-layer tin chalcogenides

Stable structures of layered SnS and SnSe and their associated electronic and vibrational spectra are predicted using first-principles DFT calculations. The calculations show that both materials undergo a phase transformation upon thinning whereby the in-plane lattice parameters converge to a pseudo-cubic phase, similar to the high-temperature behaviour observed for their bulk counterparts. The electronic properties of layered SnS and SnSe evolve to an almost symmetric dispersion whilst the gap changes from indirect to direct. Characteristic signatures in the phonon dispersion curves and surface phonon states where only atoms belonging to surface layers vibrate can also be observed for these materials.


INTRODUCTION
Chalcogenides are a remarkable family of layered materials, displaying an extensive range of optical, electronic, thermal and mechanical effects [1]. They are used as phase-change materials in rewritable data storage [2], as high performance thermoelectrics [3,4], and as absorbing layers in photovoltaic cells [5][6][7]. Lead (Pb) chalcogenides and their alloys have been heavily studied for their excellent thermoelectric properties [8][9][10][11] but the presence of toxic chemical elements is a major industrial disadvantage. Great efforts have been invested recently in a less toxic analogue: Tin chalcogenides (SnX with X=S, Se, Te).
Researchers are currently investigating emerging properties in 2D confined systems aiming at their exploitation in a wide variety of applications [12][13][14][15]. Monochalcogenides are naturally layered compounds, which can be grown in the form of flakes only a few mono-atomic layers thick [16]. The reduction of the dimensionality of the crystal (3D to 2D material) has an impact on the geometry of the layers and hence on elastic, electronic and vibrational properties, which vary with the number of layers. Recent theoretical studies of electronic structure changes in few-layer SnX have been reported, indicating that the band gap expands significantly (from 1.32 eV indirect to 2.72 eV direct) as the number of layers is decreased [17]. Deb and Kumar reproduced this study taking into account the relaxation of the atomic positions and the unit cell [18]. Mehboudi et al. also studied electronic and optical properties for mono and bilayer [19].
In what follows, the structural, vibrational and electronic properties of SnS and SnSe are studied to gain insight into the behaviour of these compounds in fewlayer form. In particular, we focus on the possible structural transformations, variations of the vibrational spectra and electronic structure resulting from reduction of the dimensionality of the crystal. Both Raman spec-troscopy and reflectivity measurements offer easy and non-destructive characterization methods: a comparison with the computationally predicted spectra can be used to assess the material thickness. We represent the results as a function of 1/n where n is number of layers. This gives us a easy way of including the bulk properties (n = ∞) and visually appreciate the convergence of the different properties towards this limit.

FEW-LAYER STRUCTURES
Tin chalcogenides exist in an orthorhombic (Pnma) bulk structure comprising weakly coupled layers of covalently bound Sn-X (X = S or Se) atomic bilayer units ( Figure 1). Both materials can be isolated in few-layer form [20]. In the following, "layer" will refer to the natural atomic bilayer unit. We consider the structural distortions appearing when this compound is isolated in free standing slabs, with thicknesses from 1 to 6 layers, and compare the geometry of the slabs with the bulk. The binding energy of SnS and SnSe slabs are 30 and 10 meV/Å 2 respectively (see supplemental information) which classifies these compounds as "easily/potentially exfoliable" according to the criteria used in Mounet et al. [21].
We perform Density Functional Theory (DFT) calculations where the exchange correlation terms are calculated within the Generalized Gradient Approximation (GGA) of Perdew et al. [22]. The structures of Sn-X bulk are obtained by fully relaxing the internal positions and the shape of the unit-cell. Relaxed    Few-layer structures are relaxed, and the lattice parameters are allowed to vary in the in-plane directions. Details of the calculations can be found in the Methods section. As can be seen in Fig. 2, when reducing the number of layers, the in-plane lattice parameters (a and b) converge, and are almost identical for the monolayer case of SnSe. This behaviour mirrors experimental stud-ies on heterostructures containing SnSe slabs [24]. The monolayer of SnS also rectifies, but its a/b ratio does not converge to 1: we find a small dynamical instability for square in-plane parameters (see Supplemental information). Condensing these unstable modes at Γ creates a distortion on the atomic coordinates along the a axis, breaks the symmetric 0.25/0.75 reduced coordinates, and leads to a/b = 1. The resulting symmetry of the crystal is triclinic. The difference in total energy between the two structures is 9 meV/f.u., and it is entirely possible that epitaxy or other substrate constraints will stabilize the square lattice of the SnS monolayer as well.  The internal coordinates and the interlayer distance of the atoms also evolve with the slab thickness. To visualize the evolution of the internal coordinates of the atoms, their projections on the YZ and XZ planes have been superimposed in Fig. 1. The positions of the upper left atoms have been aligned and only the innermost layer of the slab is shown. The YZ projection shows that the bond angles deviate from 90 • when the number of layers increases, to acquire the familiar accordion shape of bulk. The change in the XZ plane is much weaker. The internal coordinates show surface effects already in the 3 layer case, and inner layers converge towards the shape of the bulk (Fig. 1) as the thickness of the slab increases. Surface layers converge, instead, towards coordinates between those of the bulk and the monolayer.
Because of the attractive electrostatic interaction between the layers, the interlayer distance decreases with the number of layers following a 1/n law (Fig. 3). The interlayer gap is significantly smaller in the SnS case, and varies less. The interlayer distances are measured between closest atoms belonging to adjacent layers as depicted in the right panel of Fig. 3. Starting from 4 layers, differences can be seen between those on the edge and those in the center of the slab. The distances reported in Fig. 3 Figure 3. Interlayer distance between of SnX layers. Right panel: the distance is measured between the atoms closest to each other in adjacent layers (Sn atoms).
We notice a subtle difference between the two compounds: the puckering of the surface (Sn and X atoms are not on the same plane of the surface of a layer) slowly disappears when reducing the number of layers in SnSe while it stays constant in SnS. This difference shows distinct hybridization of the (on average) half-filled p orbitals at the surface, which is closer to pure p x,y,z in SnSe, leading to a quasi-cubic structure, as in α-Po [25,26] or high pressure Ca [27]. This may be linked to the size of the surface lone pairs, or to the relative alignment of the orbital energy levels.

ELECTRONIC BANDSTRUCTURES
Electronic band structures were calculated for bulk and few-layer materials using the relaxed lattice parameters reported above. Comparison with experimental bulk values shows an underestimation of the fundamental gap, which is expected for DFT calculations, but can often be corrected with a constant factor. The fundamental gap was measured by Albers et al. at 1.08 eV [28] for SnS and 0.90 eV [29] for SnSe, while our results yield 0.85 eV and 0.57 eV respectively. Fig. 4 shows the fundamental and optical gaps of the SnX compounds as a function of the number of layers. The band gap increases for thinner layers, due to quantum confinement effects. Nevertheless, the fundamental gap remains indirect in our calculations, except for the monolayer were the gap is almost direct (the difference between fundamental and optical gap is less than 0.03 eV). The reduced energy difference between indirect and direct band gaps will result in a sharper onset of optical absorption.
A comparison between our band structures for relaxed lattice parameters and fixed bulk lattice parameters [17], shows that the X-M-Y path presents a higher degree of

Gap [eV]
(Bulk) symmetry in the relaxed case. This is directly related to our prediction that the in-plane lattice parameters converge as the number of layers is reduced, and the associated atomic rearrangement tends towards a more cubic structure. Globally, we are able to tune the band gap over a large range ( 0.6 and 0.4 eV for SnS and SnSe, respectively) by changing the thickness of these compounds, which is very useful in optoelectronic applications. The electronic dispersions for both compounds can be found in supplemental information.

RAMAN AND REFLECTIVITY SPECTRA
Phonon dispersion curves were calculated for the bulk and the few layer compounds. In the monolayer case, we find a small unstable phonon mode at Γ that breaks the 1/4-3/4 symmetry of the reduced coordinates in the X axis. In the other structures, the absence of any imaginary values confirms the dynamical stability of our predicted structures. The band structures split between lowand high-frequency manifolds. This is a signature of binary compounds where the masses of the two atoms are sufficiently different [30]. In SnS, the gap in the phonon dispersion curves is more pronounced than in SnSe, due to the larger difference in mass between the two atoms.
By projecting the mode eigenvectors over the atoms and identifying them with a color code (Fig. 5), we identify surface modes which are not present in the bulk compound. Indeed, the red curves in the phonon spectra represent modes from the S atoms at the surface. We clearly identify the red curves as modes where only S We identify a surface state at 120 cm −1 with lower energy than the corresponding bulk mode. Also, the two modes of highest energy involve only surface atoms.  Figure 6. Comparison of theoretical bulk results for SnS (left) and SnSe (right) with experimental work of Chandrasekhar et al. [31].
A useful input for experimental characterisation is the prediction of Raman and reflectivity spectra. Fig. 6 shows the comparison of our simulated Raman spectra for bulk compounds and the experimental results of Chandrasekhar et al. [31], showing good agreement with calculated frequencies being within 13% of the experimental ones. We also find agreement within 1% between the thickest SnS flake presented in the work of Li et al. [16] and our theoretical study of the bulk compounds. B 1g and B 2g phonon modes do not appear in the experimen-tal work of Li et al. [16]. In our work, their relative Raman intensities is lower than the intensity of other modes. They might be hidden on a unified graph. Figure 7 shows the evolution of the Raman spectra with the thickness. The Raman spectra change drastically between the 1-and 2-layer cases. There are fewer active peaks in the monolayer, due to its symmetry which is almost cubic. The surface states are also visible in the spectra as they are Raman active. These modes are indicated by a black arrow in Fig. 7. The frequency of the surface modes appear to be converged already for 3and 4-layer slabs. However, the relative intensity of the mode will most likely decrease with thickness, as can already be observed between 3 and 4 layers. We can relate the surface modes B 3g of SnS and SnSe at respectively 103 cm −1 and 173 cm −1 by comparing their eigenvectors. The equivalence of the induced atomic displacements are shown in the supplemental information. This is also the only mode whose intensity in perpendicular polarization is much larger than for parallel polarization, giving an easy way to identify it, with a depolarization ratio ρ = I ⊥ I is larger than 0.75. In supplemental information, we show the reflectivity of IR frequency waves normal to the surface, with electric field along the two optical axes of the crystal (as defined in Ref [32]). The different peaks correspond to polar phonon modes in resonance with the incident electric field. In all spectra, the general shape is dictated by the resonance with one mode at about 225-250 cm −1 . More polar modes resonate when the number of layers increases. Also, the difference between the in-plane components of the IR reflectivity tensor evolves with the number of layers. For the 1-layer slab, the in-plane components are almost identical due to the square in-plane lattice parameters. However, with increasing number of layers, the difference between the X and Y component spectra becomes significantly noticeable as phonon modes of SnS (SnSe) with frequencies around 100 cm −1 (75 cm −1 ) resonate only with the electric field polarized Y axis. Tancogne-Dejean et al. demonstrated that standard ab-initio formalisms fail for the out-of-plane optical response in a 2D system [33], and we do not present here the Z component of the reflectivity.

BORN EFFECTIVE CHARGES, BADER CHARGES AND DIELECTRIC TENSOR
We now turn to the dielectric response and the electronic density distribution of the SnX slabs. The Born effective charge quantifies the variation of a material's polarization when the atoms are displaced and is defined as: with E the total energy, E the electric field, τ an atomic displacement, P the electrical polarization, and F the force on the atom. The Z * also govern the frequency split between longitudinal and transverse optical modes. We find that the Born effective charge tensors stay roughly constant as a function of the number of layers (see SI). This implies that the local electronic configuration of the atoms and the fundamental bonding nature do not change significantly when the number of layers varies.
To complete our analysis of the electronic density distribution, we also calculate the Bader charges [34], which confirm that the (static) electronic charge is not redistributed either by nanostructuring. The Bader charges increase slightly with the number of layers in the slab, showing a more ionic character in the bulk (see SI). As S is more electronegative than Se, the transfer of charge is larger in SnS compared to SnSe.
Dielectric constant measurements are easy to perform and give crucial information on the electronic response of materials. Our bulk values (Fig. 8) compare well with previous experimental works [31]. The computation of dielectric properties in 2D systems requires care as our calculations are performed with periodically repeated slabs separated by vacuum. In this periodic approach, the calculated dielectric tensor contains both the contribution of the slab and of the vacuum. For 2D materials of geometrical thickness t computed in a cell of cross-plane lattice parameter c, an effective dielectric constant 2D can be derived as 2D = 1 + ( DF T − 1) c/t from the dielectric constant computed in the periodically repeated approach DF T [35][36][37][38], to be able to compare the different 2D and 3D results. We choose the geometrical thickness t geom as the distance between the two outermost atoms in a slab, plus one bulk interlayer distance -a common choice in Transition Metal Dichalcogenides (TMDs), where 2D is relatively insensitive to t. For both SnX compounds, however, the variation with t is much stronger: the calculated effective electronic dielectric tensor 2D , presented in Fig. 8, increases with 1/n. This behaviour is counter-intuitive, as should vary as 1/ E g [39], and E g increases with thinning.
This effective dielectric constant model has been applied successfully on numerous systems including TMDs [37] but has also been found to fail for monochalcogenides by Gomes and Carvalho in Ref. [36]. In our case the gap decreases and the DFT bare dielectric constant increases with thickness, as expected, but the variations do not follow a square root law with a constant prefactor, and the model dielectric constant can be smaller in the bulk than in the monolayer. In practice, this means that the effective dielectric thickness is not a simple function of the geometric thickness, and is super-linear for thin slabs. We have verified that this is not simply correlated to the extension of the electronic density outside the surface, which we find to be very similar in all slabs. In Fig. 9 (a), we present a comparison of the geometrical thickness t geom and an effective thickness t ef f defined as the minimum thickness of the different layers required to retrieve the linearity between 2D and 1/ E g . In Fig. 9 (b), we also present the resulting electronic dielectric constant rescaled with t ef f , which recovers (by construction) the intuitive relation between and E g . For the same reasons as for the reflectivity, we do not report the out-of-plane component of .

SUMMARY & CONCLUSIONS
We calculate electronic and phononic properties of SnS/SnSe slabs from 1 to 4 monolayers in order to pinpoint spectroscopic signatures of 2D material thickness. We propose several simple non-invasive techniques to dis-criminate mono and few layer Sn chalcogenide samples. The behaviour of these two iso-electronic compounds is globally similar upon nanostructuring, but presents subtle differences. We identify structural distortions of the unit cell, as the thickness of the slab is reduced, leading to an almost cubic symmetry for the monolayers. The structural evolution of the unit cell with nanostructuring drives the transition from indirect to direct band gap as the number of layers is reduced, and the optical band gap expands. Surface phonon modes are identified by projecting the phonon bands over the atoms. They can be associated to a B 3g bulk mode through a comparison of the eigenvectors of the respective phonon modes. The unique depolarization ratio of these modes allows us to identify them in the Raman spectra. Furthermore, the monolayer spectra shows a specific feature : the number of active modes are reduced compared to thicker layers, which can then be used to distinguish between a monolayer, a fewlayer slab or a thicker sample. Reflectivity spectra show a similar evolution and can also be used to determine the thickness. The local electronic environment of the atoms are almost independent of layer number, as quantified by the Bader (static) and Born (dynamical) charges. Finally, the electronic gap and phononic properties have a strong thickness dependence, and common models for the effective 2D dielectric constant break down, due to a super-linear variation of the effective dielectric thickness. The results have the potential to enable fast recognition of ultrathin chalcogenide samples, and we hope to stimulate experimental work on the dielectric properties of these systems.

METHODS
Density Functional Theory (DFT) calculations are performed using the ABINIT package [40,41], which implements the plane-wave methodology (here using normconserving pseudopotentials). The exchange-correlation energy is given by the generalized gradient approximation (GGA) of Perdew, Burke and Ernzerhof [22]. Normconserving Troullier-Martins type pseudo-potentials generated with fhi98PP code are used to describe interactions between atomic cores and valence electrons of SnSe and we use an ONCVPSP [42] generated pseudopotential for SnS which produces lattice parameters that compare well with experiment [3,31].
We have checked that including Van der Waals interactions within the Grimme D3 approximation [43] does not affect the interlayer distance significantly, and it is not employed for the results above.
The wave functions are represented in a plane-wave basis with a cutoff energy of 30 Ha for SnSe and 40 Ha for SnS. The reciprocal space of SnSe bulk is sampled with a 4 × 4 × 4 Monkhorst-Pack-type grid [44], whereas an 8 × 8 × 8 unshifted grid is used for SnS bulk. The total energy is converged to within 3 meV per unit cell with respect to the k-point grid and cutoff kinetic energy of the plane waves used as a basis set. Atomic positions and lattice parameters are relaxed using a Broyden [45] algorithm until the maximal absolute force on the atoms is less than 10 −6 Ha/Bohr. The phonon band structure along high symmetry lines is obtained by standard methods based on response function calculations and DFPT [46]. Ten irreducible q-points from an unshifted 4 × 4 × 4 grid are used for the calculation of dynamical matrices. Electron band structures are calculated on a finer 24 × 24 × 24 k-point grid. For few layer calculations, convergence with respect to the size of the vacuum in the unit cell is also performed, to within 1 meV with a vacuum gap of 20 Bohr. Also, the dispersion along the vacuum is suppressed by considering k-point and q-point grids with only one point along the Z axis.
Raman intensities are computed using perturbation theory to calculate the third derivative of the energy with respect to two electric fields and one atomic displacement [47]. In our simulations, the energy of the incident light is chosen to be 2.41 eV (514,8 nm) and the temperature 300 K [48]. We use a common approximation of calculating the third derivative only within the local density approximation. The width of the Raman peaks is mainly determined by anharmonic scattering, which limits the phonon lifetime. We do not consider this effect here, and broaden the peaks with a Lorentzian function having a fixed width of 5 10 −6 Ha. The Raman tensor is averaged to represent a powder sample as in Ref. [49].