Identifying Raman modes of Sb2Se3 and their symmetries using angle-resolved polarised Raman spectra

Stephenson Institute for Renewable Energy Liverpool, Peach Street, Liverpool L69 7ZF, University College London, Department of Gordon Street, London WC1H 0AJ, UK School of Engineering, London South Bank Departamento de F́ısica and i3N, Universid Diamond Light Source Ltd., Diamond H Campus, Didcot, OX11 0DE, UK † Electronic supplementary informa 10.1039/d0ta01783c ‡ These authors contributed equally to th Cite this: DOI: 10.1039/d0ta01783c


Introduction
Antimony selenide (Sb 2 Se 3 ) is an anisotropic semiconductor with a room temperature band gap of 1.18 eV. 1 It is currently a highly attractive candidate for low-cost, large volume photovoltaic electricity generation. 2 Indeed, single junction power conversion efficiencies have reached 9.2% 3 from 2.1% 4 in only 6 years. Pushing the solar cell efficiency to over 10% would make Sb 2 Se 3 viable for industrial scale up. It is also of interest for water splitting, as a superconductor and as a thermoelectric material. [5][6][7] Nevertheless the state of knowledge of its fundamental properties has been slow to develop. Sb 2 Se 3 has an orthorhombic structure under atmospheric pressure and belongs to the Pbnm space group and D 2h point group (Schoenies notation) with lattice parameters a ¼ d 100 ¼ 11.6311Å, b ¼ d 010 ¼ 11.7808Å, c ¼ d 001 ¼ 3.9767Å. 8 The structure is comprised of 1D nanoribbons coordinated via van der Waals interactions. 9 The three main crystal planes are shown in Fig. 1. The Pbnm space group nomenclature is used in this paper, meaning the [001] direction is parallel to the ribbon growth direction. Both Pbnm and Pnma space groups can be found in the discussion of Sb 2 Se 3 which may lead to confusion (Pnma convention aligns the [010] direction with the ribbon axis).
Raman spectroscopy is a fast and usually non-destructive technique which can be used to characterise Sb 2 Se 3 and its impurity phases. 11 While gaps in the literature remain, some recent analysis of Sb 2 Se 3 via Raman spectroscopy has furthered the current understanding of this material. In 2018, Shongalova et al. resolved an issue of misinterpretation of one of the peaks in the analysis of previous Raman spectra. 11 While much of the Fig. 1 Schematic of (100), (010) and (001) crystal planes of Sb 2 Se 3 . Images generated using Vesta software. 10 literature indexed the peak at 250 cm À1 as characteristic of Sb 2 Se 3 , it was found to be an antimony oxide (a-Sb 2 O 3 ) peak. This occurred mainly due to sample degradation under high laser power during measurement but these peaks are indistinguishable from Sb 2 O 3 12,13 or elemental selenium 14 impurities already present in the sample from the growth/synthesis. The following year, the process of comprehensively indexing the Raman active vibrational modes of Sb 2 Se 3 was begun by Vidal-Fuentes et al. 15 including a low temperature study. However, due to the small size of the crystallite analysed, its orientation could not be determined and the symmetries of the vibrations could not be assigned fully.
Raman mode symmetries describe the nature of the vibrations in a material. Sb 2 Se 3 is known to be isostructural to Sb 2 S 3 which has 60 zone centre phonons of symmetry G ¼ 10A g + 5B 1g + 10B 2g + 5B 3g + 5A u + 10B 1u + 5B 2u + 10B 3u . Both A and B modes are non-degenerate with A modes showing symmetric behaviour and B modes being anti-symmetric. The subscript g indicates the mode is symmetric under inversion while u indicates it is anti-symmetric. In Sb 2 Se 3 , the A u modes are silent and the B u modes are IR-active. The set of Raman active modes have the following symmetries predicted by group theory: G ¼ 10A g + 5B 1g + 10B 2g + 5B 3g .
Angle-resolved polarised Raman spectroscopy has been used in the literature to assign vibrational symmetries to phonon modes, for example in ZnO, NbSe 3 , MoTe 2 and BaTiO 3 crystals. [16][17][18][19] For anisotropic crystals, birefringence has an impact on the scattering of polarised light described by Kranert et al. 20 However, this effect is negligible for highly absorbing materials such as Sb 2 Se 3 hence the more straightforward approach described in the results continues to hold.
Optimum Sb 2 Se 3 solar cell device performance is currently achieved by growing the ribbons perpendicularly to the contacting layers i.e. preferential growth in the [001] direction. 21 Angle-resolved polarised Raman spectroscopy has been used to determine the orientation and variation thereof across polycrystalline, thin lm absorber, CuInSe 2 . 22 The same approach could be used to identify the orientation of Sb 2 Se 3 thin lms and even identify the variation from grain to grain which may signicantly affect device performance. Vidal-Fuentes et al. also proposed this, however it has not been realised to date. 15 In this way rapid determination of the dominant orientation of lms of Sb 2 Se 3 would streamline the fabrication and characterisation process.
In this study, a combined experimental, group theory and density functional theory (DFT) approach is presented to comprehensively measure and identify the Raman bands in Sb 2 Se 3 . This study differs from earlier work in that the experiments are for a set of oriented (100), (010) and (001) crystallographic surfaces which has enabled the symmetry of the vibrational modes to be assigned reliably for the rst time. Furthermore, DFT calculations were used to determine the Raman active phonon mode energies of Sb 2 Se 3 and these corroborated the experimental results as well as suggesting additional vibrational modes. The results are expected to have relevance to both the basic science and technological communities.

Experimental and computational details
Sb 2 Se 3 crystals were prepared via the Bridgman technique using a single-zone vertical furnace. The details can be found in the ESI (ESI Section I †). Crystals were grown up to 4 mm in diameter and 1 cm in length. The crystals demonstrated bright cleavage planes extending over their full width and displayed narrow Xray diffraction (XRD) peaks (see Table S1 †). The crystals could be easily manually cleaved to expose the (010) plane and reveal reective facets. 8 This is consistent with the expectation that the weakest van der Waals interactions lie along the [010] axis. 23 Straight, parallel lines observed using optical microscopy on the surface of the cleaved planes (see Fig. S1 †) result from cleavage steps and are expected to run parallel to the 1D nanoribbons ([001] axis), as this is the axis with the most strongly bound layers. These were used as guides to orient the crystal optically on a goniometer and cut to reveal the (100) and (001) planes without the need for Laue equipment.
The cleaved (010) plane and cut (001) plane can be seen in Fig. 2a. Cutting the (001) plane was found to be most challenging as it was mechanically less stable than the other orientations. Crystals of this orientation were therefore mounted in epoxy to prevent unwanted cleaving during polish/ storage. The surfaces of the cut (100) and (001) plane crystals were polished with 0.3 mm alumina particles in de-ionised water. The cleaved (010) faces did not require polishing. To verify the surface orientations of the crystals, XRD (q-2q scan) was performed using a Rigaku Smartlab diffractometer in parallel beam geometry using a monochromated Cu-Ka 1 X-ray source.
Angle-resolved polarised Raman spectroscopy (backscattering geometry) was used to assign the vibrational symmetries to the Raman modes. The setup is shown in Fig. 2b. Linearly polarised light impinges on the sample. The scattered light is passed through a polariser oriented either parallel (solid arrow) or perpendicular (dashed arrow) to the polarisation of the incident beam. The geometry can be described in Porto notation 24 where Z(XX) Z implies the photons are incident in the Zdirection and directly backscattered light ( Z) is detected. The labels in parenthesis indicate the incident beam polarisation is in the X-direction and X-polarised light is detected (parallel polariser geometry). The crossed polariser geometry, Z(XY) Z was used also.
The crystals were rotated in the polarisation plane being studied while the Raman peak intensities were monitored. The sample angle (q) was varied between 0 and 180 in 10 increments and a spectrum was acquired at each yielding 19 spectra in each dataset. A dataset was taken for each of the three crystal planes in both parallel and perpendicular polariser congurations. A total of 114 spectra were recorded.
Room temperature Raman spectroscopy was carried out using a Renishaw in-Via microscope with a 532 nm laser. The maximum laser power (at the sample) was limited to 0.1 mW with maximum power density of 3.2 kW cm À2 to prevent sample damage. No degradation effects (peak at 250 cm À1 or higher) were observed in any of the spectra acquired. The power density at which oxidation occurs varies across studies 11,15 and was observed to vary even for different crystal planes. Analysis of this behaviour requires further work and is beyond the scope of this study but may be dependent on the wavelength of the illuminating laser and/or surface termination of the sample. A lowpass lter with a cut off at 100 cm À1 relative to the laser excitation wavelength was used. The scan time was 30 s with 20 accumulations for each spectrum. A silicon reference sample was used to calibrate the instrument and the Raman peak at 520 cm À1 had a FWHM of 4.3 cm À1 . A 1800 lines per mm grating was employed and the error in position of the tted peaks was estimated to be AE1 cm À1 . Calibration spectra with the polariser in place were taken before and aer the crystal rotational measurements and a consistent upward dri of $1 cm À1 was noted between the rst and nal measurement.
Lorentzian peaks were tted to the spectra in each dataset using a global tting method. Values for the peak positions, FWHM and peak intensities were extracted by a self-consistent tting methodology using the "global t" feature of the Origin 2019 soware. A linear background was subtracted from all spectra before tting. Seed values for peak position, FWHM and peak intensity were input with the positions and FWHM set as shared parameters. For each equivalent peak across the spectra, a single value of position and FWHM was calculated per dataset (i.e. 19 spectra) while the peak intensity was allowed to vary from spectrum to spectrum. The peak position values were not constrained to a given range during the tting. Constraining the parameters would force consistency across the datasets. The positions were kept as shared parameters within a dataset but free to vary between datasets. The resulting values were compared across datasets to determine the accuracy and reliability of the tting procedure.
Group theory calculations were performed to evaluate the sample angle dependence of the intensities of the Raman active modes. These dependences were calculated using a matrix method for both the parallel and perpendicular polariser geometry. The full description can be found in ESI Section II † and the outcomes are described in the results.
Periodic DFT calculations provided theoretical prediction of the Raman vibrational energies using scalar relativistic potentials within the Vienna Ab initio Simulation Package (VASP). [25][26][27][28] The Perdew-Burke-Ernzerhof generalised gradient approximation (GGA) exchange-correlation functional revised for solids (PBEsol) 29,30 was used for both structural relaxation and phonon supercell calculations. Interactions between the valence and core electrons were described via the projector augmented wave method. 31 The Phonopy 32 package was used to generate twenty 2 Â 4 Â 2 (320 atom) nite displacement supercells from the PBEsol relaxed structure to obtain the phonon frequencies of Sb 2 Se 3 . In all calculations, the total energy was converged to within 1 Â10 À8 eV and the forces per atom were converged to within 1 Â10 À4 eVÅ À1 . A plane wave cut off of 500 eV was used throughout and a G-centred k-point mesh of 2 Â 2 Â 2 was used for the supercell calculations. Plotting of the phonon dispersion curve was performed through the use of the sumo package. 33 Previous computational studies on Sb 2 Se 3 have either used the PBE or HSE06 functionals (with or without dispersion corrections) to describe its structural properties. [34][35][36] PBEsol optimises the structure to within 0.5% of experimental values for the a and c lattice parameters while underestimating the b parameter by just over 2%. While HSE06 gives an optimised geometry to within 2% of experiment, 37 PBEsol is much more cost effective for the large supercells required to converge the phonon frequencies. PBEsol has also shown improved performance over PBE in the calculation of vibrational properties in other chalcogenide materials. 38,39 While electronic structure calculations of Sb 2 Se 3 have been conducted by numerous groups and calculation of the phonon dispersion has also been reported, 40,41 tabulation of predicted phonon mode energies via DFT calculations and comparison to experimental data has not been published to date. Fig. 3 shows XRD patterns of the cut (100) and (001) and cleaved (010) crystal planes. In all cases, the rst structure-factorallowed reection is present followed by the higher order peaks. This conrms that the use of the cleaving method to identify the (010) plane and cutting relative to the [001] axis is a reliable means of orienting Sb 2 Se 3 on its orthorhombic faces. This method may also be successful for other van der Waals crystals.

Results
A summary of the results of group theory analysis is presented in Table 1. The constants A-F represent Raman tensor elements. The Raman peak intensity variation, with respect to sample rotation (I-q), of each vibrational mode (A g , B 1g , B 2g and B 3g ) was determined. This peak intensity variation depends on the crystal plane being analysed (I 100 , I 010 or I 001 ) and whether the polariser was parallel or perpendicular to the incident laser polarisation (I k or I t ). In the parallel polariser geometry, all A g mode intensities in each of the three crystal planes have a periodicity of 180 . The B g mode intensities will repeat aer 90 due to their 2q dependence. In perpendicular polariser geometry, both A g and B g modes have 90 periodicity, however, A g modes show sine-like behaviour and B g modes show cosinelike behaviour. A g and B g modes can therefore easily be distinguished in both polariser orientations. In this way, knowing the crystal plane as well as the expected peak intensity variation of the individual Raman peaks allows their symmetries to be identied. Fig. 4 shows an example (the (100) plane, parallel polariser) of the angle-resolved datasets collected and the tting procedure used for analysis. The black line indicates the raw data while the coloured peaks show the tted Lorentzian lineshapes. Similar plots for the (010) and (001) planes are shown in Fig. S2 † and a more detailed example t of a single Raman spectrum can be found in Fig. S3. † For each peak in the gure, the peak intensity (I) vs. sample angle (q) was extracted and plotted. This data was compared and tted to the expected I-q behaviour from Table 1. Fig. 5 shows the data points for the 185 cm À1 peak I-q plot. The solid lines are ts to the data points of the relevant functions from Table 1. For example, the 185 cm À1 peak has maximum intensity at 0 rotation, minimum at 90 and returns to its maximum at 180 rotation of the crystal. The 180 period is characteristic of an A g symmetry vibration. Therefore the function tted to the data was for an A g mode in the (100) plane with parallel polariser geometry. The tensor elements were used as tting parameters. Due to the reliable and consistent t of the expected A g function for all crystal orientations and both polariser geometries the symmetry can be conrmed as A g . These plots were prepared for each peak in each of the six datasets in order to identify the mode symmetries. The remaining plots can be found in the Fig. S5-S7. † These plots show excellent consistency in peak position with only small variations (within the expected error of 1 cm À1 ) in values  Table 1 Relations between sample angle (q) and Raman peak intensities (I) for Sb 2 Se 3 crystals with exposed (100), (010) or (001) plane in parallel (I k ) or perpendicular (I t ) polariser geometry. obtained from the different planes and polariser congurations. This further validates the rigorous analysis procedure. Fig. 6 shows the DFT-calculated dispersion of the phonon bands along the high symmetry path within the Brillouin zone. Due to the fact that Raman spectroscopy probes only zone-centre phonons the corresponding theoretical values can be noted from the intersection of the phonon bands with the G point. 30 phonon energies were predicted in this way which is consistent with the expected 10A g + 5B 1g + 10B 2g + 5B 3g modes of Sb 2 Se 3 . The peaks determined using DFT were found to be shied from the experimental positions in a systematic manner (À11 cm À1 ) which has been discussed previously for phonon energies calculated using the PBEsol functional with a similar shi found (À10 cm À1 ). 42 Such disagreement can result from a combination of differences in computational lattice parameters versus experiment. This in turn may be a consequence of the athermal limit and inaccuracies in the DFT approach in modelling the interatomic forces precisely. For a material such as Sb 2 Se 3 , where strong and weak interactions are present in different directions, modelling the interatomic forces is challenging. Nevertheless a consistent linear shi of the DFT computed peaks upwards by 11 cm À1 resulted in excellent agreement with experimental results as seen in Table 2.
A complete description of the determined vibrational modes of Sb 2 Se 3 is presented in Table 2. The experimental peaks, their symmetries, theoretically predicted modes and previously reported peak positions 15 are compared and show very good agreement. Seven A g modes and seven B g modes were assigned experimentally. Several assignments deviate from those seen by Vidal-Fuentes and are further discussed below. Vidal-Fuentes et al. 15 could not resolve the peaks found at 212, 187 and 156 cm À1 with a randomly oriented crystal. For the same reason, the B g symmetries could previously not be assigned. All other peaks seen by Vidal-Fuentes were corroborated by the experimental peak positions determined here. In addition, the partial symmetry assignments Vidal-Fuentes report are consistent with the symmetries determined in this study.
While Vidal-Fuentes report a single peak at 213 cm À1 , the experimental results here from the (001) plane indicated two closely spaced peaks (213 and 212 cm À1 ). The B 1g modes show large relative intensity in the perpendicular polariser geometry for this plane (see Fig. S7 †). Due to intensity differences, the 213 cm À1 A g mode cannot be resolved in the perpendicular polariser case and the 212 cm À1 B 1g mode is of too low intensity to be seen in the parallel polariser geometry. The existence of two peaks in this region is also consistent with the closely spaced modes predicted via DFT. Data from the (010) plane suggests two A g peaks in this region, however this second peak of A g symmetry was not observed in any other plane. It has been included for completeness but could not be condently labelled (Fig. S6 †).
For the (001) plane, the global t for the A g mode around 185 cm À1 was shied to 187 cm À1 for both the parallel and perpendicular polariser case (see Fig. S7 †) over numerous repetitions of the experiment. This is a strong indication of the existence of another mode and is again corroborated by the theoretical calculations. This also holds true for the mode at 155 cm À1 . The peak at 155 cm À1 is clearly visible in the (010) plane spectra (see Fig. S6 †) and a 156 cm À1 feature is visible in the (100) plane (Fig. S5 †). The latter could be an artefact from the collection cone angle of the objective lens which does not uniquely collect light scattered normal to the sample surface. 19 However the consistency of the results with theory and literature suggest this is not the case.
Once a shi of +11 cm À1 was applied, the DFT calculations predicted all of the experimentally observed modes to a high degree of accuracy. Some modes which were not observed experimentally were also predicted. In contrast to the calculations, no modes around 140 cm À1 were seen in the Raman spectra which is consistent with the experimental data from Vidal-Fuentes et al. This discrepancy may be due to the low intensity of these vibrational modes. The remaining predicted modes show very good agreement with those experimentally observed.  5 Raman peak intensity plots, I-q, for the 185 cm À1 peak for each crystal plane in parallel and perpendicular polariser configuration. The points represent the experimental data while the solid line is a fit to the data using the relevant dependences in Table 1. The 110-140 cm À1 region was difficult to analyse in the (100) and (001) planes as there are many overlapping modes. The broad nature and low intensity of the features in this region prohibit accurate global ts to be calculated. Individual spectra in each crystal orientation show a clear peak at 117 cm À1 (see ESI Fig. S2-S4 †). This mode was assigned as A g as only these modes are expected to be visible in spectra for all three crystal orientations. However, this assignment should be considered tentative on the present evidence, even though it is consistent with both Vidal-Fuentes and DFT predictions. Nevertheless, the additional peaks predicted by DFT in this region are expected to be present although they cannot be resolved in the data measured at room temperature.
The peak at 100 cm À1 lies at the edge of the Raman spectrometer lter cut off and thus the position may be articially high. This mode showed strong A g vibrational symmetry in all parallel polariser congurations and was therefore included. The cut off using the perpendicular polariser case was slightly higher so that the mode at 100 cm À1 was no longer visible and could not be analysed. The lter cut-off generally limited the analysis of the vibrations at wavenumbers below 100 cm À1 . The mode at 82 cm À1 is included in the table, with a tentative assignment of symmetry, even though its intensity was reduced by the cut off of the low-pass lter. Nevertheless this peak remained visible in the spectra and displayed a peak intensity variation consistent with an A g mode (Fig. S7 †). Further work could be done with a Raman spectrometer congured to measure below 100 cm À1 using the same analysis technique to identify and characterise the low wavenumber vibrations. While the orientation of the ribbons is difficult to unequivocally deduce from the Raman spectra, it is clear that the device relevant (001) plane shows little to no Raman scattering intensity around 155 cm À1 . Minimisation of this peak may provide a fast approach to testing the alignment of the ribbons in thin lms and crystals. In the crossed polariser geometry, altering growth parameters to maximise the 205 cm À1 peak intensity (see Fig. S7 †) can also ensure growth in the [001] direction. Raman spectroscopy already provides a fast, non-destructive approach to mate †rial characterisation and this could be further exploited.

Conclusions
van der Waals crystals such as Sb 2 Se 3 may easily cleave, revealing cleavage steps. These may be used in optical orientation of crystals for cutting without the need for X-ray methods. Angle-resolved polarised Raman spectroscopy of individual crystal planes was used to identify the Raman modes of Sb 2 Se 3 and assign vibrational symmetries to each peak in the spectrum. A total of ten A g , ve B 1g , ten B 2g and ve B 2g modes are expected from group theory and DFT calculations. Experimentally, seven A g , two B 1g , three B 2g and two B 3g modes were identied. Results are consistent with those of Vidal-Fuentes, 15 but additional peaks were identied and concrete symmetry assignments made through the use of oriented crystals in this work. DFT calculations supported and veried the phonon mode energies as well as predicting some low intensity and lower wavenumber modes. Some low energy modes were not observable experimentally, but are expected to be measured in further studies. A methodology was reported to quickly identify the technologically important (001) orientation using the 155 cm À1 Raman peak. This study both strengthens the fundamental understanding of Sb 2 Se 3 as a material and paves the way for rapid determination of lm orientation and grain to grain variation of growth.

Conflicts of interest
There are no conicts to declare.