Quantitative Analysis of Temperature Dependence of Raman shift of monolayer WS2

We report the temperature-dependent evolution of Raman spectra of monolayer WS2 directly CVD-grown on a gold foil and then transferred onto quartz substrates over a wide temperature range from 84 to 543 K. The nonlinear temperature dependence of Raman shifts for both and A1g modes has been observed. The first-order temperature coefficients of Raman shifts are obtained to be −0.0093 (cm−1/K) and −0.0122 (cm−1/K) for and A1g peaks, respectively. A physical model, including thermal expansion and three- and four-phonon anharmonic effects, is used quantitatively to analyze the observed nonlinear temperature dependence. Thermal expansion coefficient (TEC) of monolayer WS2 is extracted from the experimental data for the first time. It is found that thermal expansion coefficient of out-plane mode is larger than one of in-plane mode, and TECs of and A1g modes are temperature-dependent weakly and strongly, respectively. It is also found that the nonlinear temperature dependence of Raman shift of mode mainly originates from the anharmonic effect of three-phonon process, whereas one of A1g mode is mainly contributed by thermal expansion effect in high temperature region, revealing that thermal expansion effect cannot be ignored.

Since the discovery of graphene, atomically thin two-dimensional (2D) layered materials have drawn intense attention due to their unique physical, chemical, and mechanical properties 1,2 . 2D transition metal dichalcogenides (TMDCs) such as MX 2 (M = Mo, W and X = Se, S) have been successfully synthesized and prove to be one of the most stable atomically thin 2D materials 1,3 . Few layer TMDCs are stacked by several monolayers interacted by weak Van der Waals forces. Monolayer TMDCs consist of a triangular or hexagonal plane of transition metal M atom sandwiched by two triangular layers of dichalcogenides X atom 4 . Unlike graphene, monolayer TMDCs don't have inversion symmetry of crystal space group so that TMDCs undergo a transition of band structures from an indirect band gap in bulk to a direct band gap in a monolayer material 5 , making them useful in nanoelectronic device applications [6][7][8] . For example, monolayer WS 2 has a direct band gap of 2.1 eV, while the bulk WS 2 has an indirect gap at 1.3 eV 9 . Among all TMDCs, monolayer MoS 2 was first synthesized and has been studied extensively. Its field-effect transistors (FET) have been demonstrated and exhibited large current on/off ratio 10 . In contrast to MoS 2 , WS 2 was synthesized later and studied much less. It was reported recently that monolayer WS 2 had intense photoluminescence demonstrating the characteristics of its direct band gap 11 and giant spin-valley coupling 12 , which implied potential applications in light emission, optical sensors and spin valleytronics applications 1,2 .
The integrated photoelectronic devices based on monolayer WS 2 will have extremely high degree of integration, and hence heat transfer and phonon behaviors are very important for heat management of nanointegrated devices. Raman spectroscopy is a powerful tool and has been used widely to non-destructively characterize the structure, symmetries, and optical phonon behaviors of nanomaterial. The position and width of Raman scattering peak can reflect vibrational frequency and dynamics of optical phonons, respectively, while the latter is directly related to the heat diffusion rate. In nanoelectronic device application, it is very important to understand the effect of phonons because the self-heating of the device can significantly affect the performance. The temperature dependence of Raman shift of monolayer and few layer MoS 2 has been extensively studied over a heating process 13,14 , a cooling process 15 , and a wide temperature range from 77 K to 623 K 16 . Thermal conductivity of MoS 2 was given out based the temperature dependence of Raman shifts 16,17 . The effect of different-type substrates on Scientific RepoRts | 6:32236 | DOI: 10.1038/srep32236 temperature dependence of Raman shifts was also reported from room temperature to 500 °C 18 . In contrast, temperature dependence on Raman shifts of WS 2 was reported sparsely 3,[19][20][21][22] , so that thermal conductivity of monolayer WS 2 was reported sparsely too 20 , and thermal expansion coefficient has not been reported yet. Thripuranthaka et al. reported the first experimental investigation on temperature-dependent Raman shifts of mechanically exfoliated monolayer WS 2 transferred onto a Si substrate over a wide temperature range from 77 K to 623 K 19 , and observed an obvious nonlinear temperature dependence of Raman shifts of E g 2 1 and A 1g modes. However, authors ignored the nonlinear dependence, and gave out a small first-order temperature coefficient of − 0.006 cm −1 /K for both E g 2 1 and A 1g modes by simply linear fitting to temperature-dependent Raman shifts 19 . Peimyoo et al. studied temperature dependent Raman shifts of mono-and bi-layer WS 2 grown on Si substrates by chemical vapor deposition (CVD) over a low temperature range from 80 to 380 K 20 , but observed a good linear temperature dependence of Raman shifts of E g 2 1 and A 1g modes, and gave out a large first-order temperature coefficients of − 0.0125 and − 0.0149 cm −1 /K, respectively for E g 2 1 and A 1g modes 20 . Su et al. studied the temperature dependence of Raman shifts of CVD-grown monolayer WS 2 films onto SiO 2 /Si and sapphire substrates as well as transferred on SiO 2 substrates over a high temperature range from 25 to 500 °C 21 , and found good linear temperature dependences of Raman shift of E g 2 1 mode of all samples except for the WS 2 grown on SiO 2 /Si substrates, but complex nonlinear temperature dependences of Raman shift of A 1g mode of all sample, showing strong dependence of Raman shifts on substrate types and bonding between WS 2 and substrates. Gaur et al. studied the temperature dependence of Raman shifts of CVD-grown monolayer WS 2 on Al 2 O 3 substrates over a wide range from 83 to 573 K 22 , and observed a weakly nonlinear temperature dependence of Raman shift of A 1g mode. Those reports mentioned above presented diverse temperature dependences on Raman shift of E g 2 1 and A 1g modes, or inconsistent temperature dependences to each other. On one hand, the diverse temperature dependences may show strong dependence of Raman shift on sample-prepared methods, substrate types and bonding between WS 2 sample and substrates. Meanwhile, it was also implied that the reliability of temperature dependences reported may need to be confirmed further, and hence more experimental studies are necessary very much to extract the accurate temperature coefficient of Raman shift of Raman active modes because the temperature coefficient directly reflects the strength change of Raman vibration bond with varying temperature. It is also an important parameter to differentiate layer number of layered films. Moreover, quantitative analysis of temperature dependence of Raman shifts using physical models is also absent so that physical origin of nonlinear temperature dependence is not clear. Even the experimental value of the thermal expansion coefficient of monolayer WS 2 has not been available so far so that some theoretically calculated values of the in-plane modes for monolayer WS 2 23-25 could not be verified experimentally. Moreover, experimentally the only reported values of thermal expansion coefficients of bulk 2H-WS 2 are also puzzled because the value of in-plane mode was as twice more as one of out-plane mode 26 .
In this article, we investigate temperature dependence on Raman shift of monolayer WS 2 27 , directly CVD-grown on a gold foil and then transferred onto quartz substrates over a wide temperature range from 84 to 543 K. To our knowledge, temperature dependence of Raman shift of monolayer WS 2 with the combination of such a sample and substrate is studied for the first time. The nonlinear temperature dependence of Raman shifts for both E g 2 1 and A 1g modes has been observed. A physical model, including thermal expansion and three-and four-phonon anharmonic effects, is used quantitatively to analyze the observed nonlinear temperature dependence. Thermal expansion coefficient of monolayer WS 2 is extracted from the experimental data for the first time. It is found that thermal expansion coefficient of out-plane mode is larger than one of in-plane mode, being more reasonable physically. It is also found that the nonlinear temperature dependence of Raman shift of E g 2 1 mode mainly originates from the anharmonic effect of three-phonon process, whereas one of A 1g mode is mainly contributed by thermal expansion effect in high temperature region, but still by three-phonon anharmonic effect in low temperature range. However, thermal expansion effect was ignored in the most of current reports 20-22 , obviously being not justified.

Results
Monolayer WS 2 sample studied here was grown on gold foil substrates by CVD, and then transferred onto a quartz substrate for Raman measurement and a SiO 2 /Si substrate for good optical contrast (See Fig. 1(a)). The details of the synthesis and crystal-quality characterization of monolayer WS 2 can be found in ref. 27. Large area high-quality monolayer WS 2 was grown by this CVD method on gold foils. As shown in Fig. 1(a), single crystal monolayer WS 2 on SiO 2 /Si presents homogenous triangular blue domain and has a size over 100 μ m. Here Raman spectroscopy is used further to characterize the quality and layer number of WS 2 on quartz. The monolayer WS 2 sample on a quartz substrate is mounted in a cryostat cooled by liquid nitrogen for measurement of temperature-dependent Raman shift.
A representative Raman spectrum over a wavenumber range of 100-800 cm −1 at room temperature is taken using 514.5 nm laser and shown in Fig. 1(b). The Raman spectrum consists of many first-order and second-order peaks 15,20 . The first-order peaks mainly include ones of E g resonance characteristic of 2LA(M) mode. Raman frequency difference between E g 2 1 and A 1g modes is found less than 61 cm −1 , indicating that the sample is monolayer 28,29 .
The temperature-dependent Raman spectra of monolayer WS 2 are taken over a wide temperature range from 83 K to 543 K, and are plotted in Fig. 2(a) in the form of the waterfall graph to show the variation of peak intensity clearly with the temperature, where only Raman peaks of phonon modes that we are interesting in, are displayed in the range of 340-540cm −1 including E g 2 1 , A 1g and Si-related peaks. One can find two obvious features. One is that the combined peak of E g 2 1 and 2LA modes is much stronger than A 1g peak in amplitude. The other is significantly non-monotonous variation of the amplitude of the combined peaks at ~354 cm −1 with temperature. Similar features were also reported 21,22 . the first feature originates from the significant enhancement of 2LA(M) mode under the excitation of 514.5 nm laser whose energy is in the vicinity of B exciton of WS 2 30,31 , as shown similarly in refs 21 and 22. The second feature originates from temperature dependence of B exciton energy of monolayer WS 2 32,33 . The energy of B exciton is almost resonant to photon energy of 514.5 nm laser as the temperature of WS 2 is set at ~223 K, so that Raman scattering intensity approaches maximum. As temperature deviates from 223 K, the energy of B excitons becomes non-resonant to the photon energy of incident laser so that Raman scattering intensity weakens. The more far the deviation of temperature is from 223 K, the larger the detune is between the energies of B exciton and laser photon, and hence the weaker the Raman scattering intensity is. However, there is a difference of 70 K between this work and ref. 32. for the temperature that strongest Raman scattering peaks appeared. It may just reflect the effect of different substrates and bonding between film and substrates 31 .
In order clearly to display the shift of Raman peak position with temperature, the Raman spectra normalized by the intensity of overlapped peak of E g 2 1 and 2LA(M) modes are plotted in Fig. 2(b). One can easily discern the red-shift of all Raman peaks with increasing temperature, similar to what was observed for MoS 2 monolayers [13][14][15][16][17][18] . One can also find that the contribution of 2LA(M) mode in the wide overlapped peak of E g 2 1 and 2LA(M) modes becomes weak with the increase of temperature, while one of E g 2 1 mode is enhanced.

Discussion
Multiple-peak Lorentzian sum function is used to best fit each Raman spectrum. The peak positions of multiple modes, including E g   modes that we are focusing on are plotted in Fig. 3 as a function of the lattice temperature of monolayer WS 2 . One can see that the temperature dependence of Raman shift of E g 2 1 mode is weakly nonlinear, while one of A 1g mode is strongly nonlinear.
The temperature dependence of Raman shifts of E g 2 1 and A 1g modes is preliminarily analyzed by a linear approximation. The following linear equation is used best to fit the temperature dependence in Fig. 3, as done in refs [19][20][21].
where ω 0 is the extrapolated peak position at zero Kelvin, and χ is the first-order temperature coefficient. The best linear fitting is plotted in Fig. 3 by solid lines. It gives out the first-order temperature coefficients of monolayer WS 2 as − 0.0093 (cm −1 /K) and − 0.0122 (cm −1 /K), respectively for E g 2 1 and A 1g modes, which are well in the range of the minimum (− 0.0066 cm −1 /K 19 ) and maximum (− 0.0155 cm −1 /K 21 ) reported for E g 2 1 mode, and one of the minimum (− 0.006 cm −1 /K 19 ) and maximum (− 0.0149 cm −1 /K 20 ) reported for A 1g mode, respectively. Therefore, our results provide new reference data for first-order temperature coefficients of monolayer WS 2 .
One can see from Fig. 3 the temperature dependence actually is nonlinear, similar to previous reports 19,21,22 . However, main origin of the nonlinear dependence has been unknown yet because previous reports either ignored the nonlinear dependence 19 or simply analyzed it with cubic polynomial 21 . Gaur et al. made an only physical analysis with anharmonic model, but ignored thermal expansion effect, so that the dominant origin of the nonlinear dependence was still unclear, being thermal expansion or anharmonic effects, or both of them cooperatively? Here we make the first quantitative analysis using a full model including both thermal expansion and anharmonic effects. As a result, thermal expansion coefficient of monolayer WS 2 is obtained for the first time.
Before the quantitative analysis is started, it is very necessary to discuss the temperature dependence of A 1g mode shown in Fig. 3 because it looks quite strange. The Raman shift reduces and the reduction rate increases progressively with increasing temperature in the range of below ~380 K. However, such a change trend stops at ~380 K. Then a nearly linear slow decrease starts from 380 up to 460 K. Finally, the dependence becomes almost unchanged in the range of 460 to 550 K. Such temperature dependence of A 1g mode is very similar to ones observed by Su et al. on mechanically exfoliated and CVD-grown monolayer MoS 2 18 and CVD-grown WS 2 21 transferred onto SiO 2 /Si substrates, where anomalous temperature-dependent change occurred at ~100 °C (corresponding 373 K, agreeing very well with 380 K here). The anomalous temperature dependence of A 1g mode starting at ~380 K may be explained by possible forming of wrinkles or ripples in monolayer WS 2 because thermal expansion coefficient of WS 2 is about one order of magnitude higher than that of SiO 2 21 , whereas the wrinkles or ripples can lead to significant strain in monolayer WS 2 and weakening of bonding between the sample and the substrate. It was reported that the out-of-plane mode (A 1g ) was much more sensitive to and stronger affected by the bonding between the film and the substrate than the in-plane mode (E g 2 1 ) 21 , so that wrinkles or ripples results in anomalous temperature dependence of Raman shift of A 1g mode in high temperature range, but not one of E g 2 1 mode. To avoid the effect of significant strain on the temperature dependence, we analyze quantitatively the temperature dependence of A 1g mode only in the range of below 380 K, but one of E g 2 1 mode in whole range of 83-543 K.
To understand the physical origin of these nonlinear temperature dependencies, a physical model, including thermal expansion and three-and four-phonon anharmonic effects 18,34 , is used to quantitatively analyze the nonlinear temperature dependence of Raman shifts of E g 2 1 and A 1g modes. The model can be expressed as 18 , where ∆ω E and ∆ω A are Raman shift change induced by lattice thermal expansion and pure temperature effects, respectively. Volume expansion-induced contribution to the change of Raman shift can be described by Grüneisen constant model 18 , where n is the degeneracy, 1 for A 1g mode and 2 for E g 2 1 mode, γ G is the Grüneisen parameter, and α is the thermal expansion of the material. The integration denotes the decrement of the vibrational frequency resulting from the expansion of volume. As the Grüneisen parameter (γ G ) and thermal expansion coefficient (α ) of monolayer WS 2 material is unknown experimentally in our wide temperature range, we write their product as a polynomial expression, where p 0 , p 1 and p 2 are constants, and will be obtained as fit parameters by best fitting to the nonlinear temperature dependence. The contribution from pure temperature effect mainly considers the anharmonic effects of three-and four-phonon processes. According to the viewpoint of Klemens 34 , the light scattering process can be viewed as involving the absorption of a photon, the emission of a photon, and the creation of an optical phonon which then decays via anharmonicity into two phonons, three phonons, etc. The production of two and three phonons is called three-phonon processes and four phonon processes, respectively. The pure temperature effect including three-and four-phonon processes can be described by a semi-quantitative simple model developed by Klemens 34 , where ω ω = = ћ ћ kT y k T x /2 , /3 , and coefficients A and B are constants as fit parameters, representing the contributions of three-and four-phonon processes to the frequency shift, respectively.
Equations (2)(3)(4)(5) are used best to fit the nonlinear temperature dependence of Raman shifts shown in Fig. 3. The fits are plotted in Fig. 4 by red solid lines and agree very well with experimental point data. Meanwhile, individual contribution of thermal expansion, three-and four-phonon effects is also plotted in Fig. 4. All fit parameters are extracted as A = − 0.902 ± 0.047, B ≈ 0., p 0 = (1.619 ± 2.467) × 10 −6 , p 1 = (3.523 ± 2.893) × 10 −8 and p 2 = − (4.751 ± 6.051) × 10 −11 for E g 2 1 mode and A = − 1.200 ± 0.025, B ≈ 0., p 0 = (3.471 ± 2.462) × 10 −6 , p 1 = − (9.952 ± 3.922) × 10 −8 and p 2 = (5.527 ± 1.103) × 10 −10 for A 1g mode. It is worth noting that actually only four free fit parameters exist in our fit model because parameter B is found very small so that it may be set to zero and ω 0 can be determined priorly by nonlinear fitting. The uncertainties of parameters A and (p 0 , p 1 , p 2 ) are given only if one of the two group parameter is fixed to the mean. The allowed uncertainty of parameters, p 0 , p 1 and p 2 , is larger because error mainly occurs in high temperature range, and is not distributed homogenously in the whole experimental temperature range.
One can see that the dominant contribution to nonlinear temperature-dependent Raman shift of planar mode E g 2 1 is from the three-phonon anharmonic process, while thermal expansion contributes weakly and the contribution of four-phonon anharmonic effect is completely negligible. In contrast, for A 1g mode, the contribution of thermal expansion and three-phonon anharmonic effects competes with each other though four-phonon process is negligible. In low temperature range, the latter is dominant, while the former is dominant in high temperature region. Our analysis reveals the quantitative contribution of thermal expansion, three-and four-phonon effects for the first time. The contribution of thermal expansion cannot be ignored for either A 1g or E g 2 1 mode.
The parameters, p 0 , p 1 and p 2, have been extracted in last section. Instituting them into Eq. (4), the temperature dependence of the product of thermal expansion coefficient and Grüneisen parameter, γ G α , can be achieved. If  Fig. 5(a)), whereas one of A 1g mode also agrees with three reported values and have the same temperature-dependent trend although it is significantly larger than reported ones (see Fig. 5(b)). More important is in our results that the values of E g 2 1 mode is smaller than ones of A 1g mode above room temperature, but in unique experimental report in ref. 26 the values of E g 2 1 mode is larger than ones of A 1g mode. We believe our results may be more justified physically than ones in ref. 26 because in layered materials the out-of-plane direction is confined very weakly. Furthermore, our results show fully positive values above 130 K, which is reverse fully with ones in ref. 35.
In summary, we have reported temperature dependent Raman study of the first-order E g 2 1 and A 1g mode in monolayer WS 2 sample directly CVD-grown on a gold foil and then transferred onto quartz substrates over a wide temperature range from 84 to 543 K. The nonlinear temperature dependence of Raman shifts for both E g 2 1 and A 1g modes has been observed. The first-order temperature coefficients of Raman frequency shifts are given as − 0.0093 (cm −1 /K) and − 0.0122 (cm −1 /K) for E g 2 1 and A 1g peaks, respectively. A physical model, including thermal expansion and three-and four-phonon anharmonic effects, is used quantitatively to analyze the observed nonlinear temperature dependence. Thermal expansion coefficient of monolayer WS 2 is extracted from the experimental data for the first time. It is found that thermal expansion coefficient one of out-plane mode is larger than one of in-plane mode, being more reasonable physically. It is also found that the nonlinear temperature dependence of Raman shift of E g 2 1 mode mainly originates from the anharmonic effect of three-phonon process, whereas one of A 1g mode is also mainly contributed by three-phonon process in low temperature range but by thermal expansion effect in high temperature region. Our results are useful for further experimental and theoretical studies on the thermal properties of two dimensional materials and the development of nano devices of WS 2 . They are also helpful to get deep insight into heat transfer and photon dynamics of monolayer WS 2 .

Methods
Monolayer WS 2 sample studied here was grown on gold foil substrates by CVD, and then transferred onto a quartz substrate. The details of the synthesis and crystal-quality characterization of monolayer WS 2 can be found in ref. 27. Large area high-quality monolayer WS 2 could be grown by this CVD method on gold foils. Single crystal monolayer WS 2 is triangular, and has a size over 100 μ m.
The micro-Raman system used in this study is Renishaw inVia with a Linkam TS1500 heating system. Renishaw inVia Micro-Raman system has a spectral resolution smaller than 1 cm −1 and a 50 × long working-distance lens which can focus laser beam to a spot less than 1 μ m in diameter. A 514.5 nm laser is used. The power incident on the sample is 1.97 mW, low enough to avoid heating of the sample. The Linkam TS1500 heating system has a temperature control accuracy of 1 °C, heating the sample with a step of 10 °C at a rate of 10 °C/min. To stabilize the sample temperature, ten minutes' delay is applied at each temperature step till a Raman spectrum is taken, ensuring sufficient time to reach thermal equilibrium.