Frequency-dependent ab initio Resonance Raman Spectroscopy
FFrequency-dependent ab initio
Resonance Raman Spectroscopy
G. Kukucska, V. Z´olyomi, and J. Koltai Department of Biological Physics, ELTE E¨otv¨os Lor´and University,P´azm´any P´eter s´et´any 1/A, 1117 Budapest, Hungary National Graphene Institute, University of Manchester, Oxford Road, Manchester M13 9PL, UK. (Dated: April 9, 2019)We present a new method to compute resonance Raman spectra based on ab initio level calcu-lations using the frequency-dependent Placzek approximation. We illustrate the efficiency of ourhybrid quantum-classical method by calculating the Raman spectra of several materials with differ-ent crystal structures. Results obtained from our approach agree very well with experimental datain the literature. We argue that our method offers an affordable and far more accurate alternativeto the widely used static Placzek approximation.
Raman spectroscopy is one of the most versatile non-destructive characterization methods for molecules andsolid state systems. The Raman effect allows one togauge the structural properties of materials through thefrequencies of vibrations which can be determined by thedifference between the incoming and outgoing photon en-ergy. In addition, resonance effects in the Raman spectracan reveal details of the electronic structure and opticalproperties of the examined material. The distributionof spectral weights between the different peaks in theRaman spectra can also carry information about pertur-bations in the material such as strain or doping, or evenlattice defects. Theoretical modeling of Raman spectra is an excep-tionally challenging task. Resonant processes are usuallydescribed at the semiempirical level, e.g. using the tight-binding model.
In contrast, when ab initio methods areemployed, the calculations are limited to the static ap-proximation where the matter-light interaction is approx-imated with the response to a static external electric fieldwithin the static Placzek approximation.
Ramanpeak intensities predicted by the static Placzek approxi-mation are fairly accurate for wide gap semiconductors,i.e. when the laser excitation energy is small comparedto the optical gap. However, when the gap is comparableto or smaller than the laser excitation energy, it cannotproduce accurate relative intensities any more. More-over, Raman spectra of metallic or semi-metallic systemscannot be calculated in this way, since the response to astatic external electric field in defect-free metals is diver-gent.If resonance effects are taken into account in the calcu-lation of the Raman spectrum, theory can make accuratepredictions regardless of the electronic properties.
Furthermore, by calculating laser energy dependent Ra-man intensities, resonance effects can be studied in theexcitation profile. However, available commercial ab ini-tio codes only offer to calculate Raman intensities basedon the static Placzek approximation, which limits the ex-tent to which Raman intensities can be predicted for thereasons discussed above.A few recent works employed many-body theory to I n t eg r a t ed R a m anpea k i n t en s i t y Excitation energy Time dependentperturbation theoryFrequency-dependentPlaczek approximationStatic approximation ε A ε B + ħω ν I n c o m i ng r e s onan c e O u t go i ng r e s onan c e I n t eg r a t ed R a m anpea k i n t en s i t y Excitation energy Time dependentperturbation theoryFrequency-dependentPlaczek approximationStatic approximation ε A ε B + ħω ν I n c o m i ng r e s onan c e O u t go i ng r e s onan c e A A a) b)c) FIG. 1. Feynman diagram of elastic light scattering (a) andthe Stokes Raman process (b). Solid lines indicate electronand hole propagation, vertical dashed lines with capital let-ters denote intermediate virtual states, while wavy solid andstraight dotted lines indicate photon and phonon propagation,respectively. Solid and empty vertices denote electron photonand electron phonon interaction, respectively. In (c) we showa schematic representation of how the laser energy dependenceof the Raman peak intensity changes when comparing timedependent perturbation theory with the frequency-dependentPlaczek approximation and the static approximation (c). compute the frequency-dependent Raman spectrum taking excitonic effects into account through the Bethe-Salpeter equation (BSE). While these methods are ableto provide very good accuracy for theoretical predictions,they are limited to small systems due to the extremelyhigh computational demand of many-body calculations.In this work we demonstrate a method that enablesthe efficient computation of resonance Raman spectra atthe ab initio level using existing computational software.Specifically, we demonstrate how to combine the clas- a r X i v : . [ c ond - m a t . o t h e r] A p r sical frequency-dependent Placzek approximation with ab initio level calculations as implemented in the Vi-enna ab initio Simulation Package (
VASP ). To val-idate our method we compare our results to Ramanspectra calculated using the static approximation imple-mented in Quantum Espresso (QE) and experimentaldata. As we will show, our Frequency-dependent abinitio
Placzek approximation is applicable to any mate-rial and provides a fairly accurate description of Ramanintensities, without the expense of many-body theory.To begin, we briefly introduce the theoretical basics ofour method. The frequency-dependent Placzek approx-imation is based on the linear response to a monochro-matic electromagnetic field within the long wavelengthlimit (i.e. neglecting the spatial dependence of the elec-tric field): P = ˆ α ( ω l ) E cos( ω l t ) , (1)where E is the polarization vector of the incident electricfield, P is the induced polarization and ˆ α ( ω l ) is the socalled polarizability tensor at the laser energy ω l . In theclassical description, Raman intensity is proportional tothe derivative of the polarizability with respect to phononnormal modes, resulting in the following expression forthe Stokes branch of the Raman spectrum: I ( ω, ω l ) = (cid:88) ν ω s ω ν (cid:12)(cid:12)(cid:12)(cid:12) ∂ ˆ α ( ω l ) ∂Q ν (cid:12)(cid:12)(cid:12)(cid:12) δ ( ω − ω ν )( n ( ω ν ) + 1) , (2)where ω = ω l − ω s is the Raman shift, Q ν , ω ν are thephonon normal mode and frequency, ω s is the frequencyof the scattered light, δ ( x ) is a normalized Lorentzianfunction and n ( ω ν ) is the Bose-Einstein distribution atroom temperature.In the following we compare the laser energy dependenceof amplitudes calculated in this way, with those calcu-lated from time-dependent perturbation theory. Let usconsider the expression of the polarizability correspond-ing to the Feynman diagram of elastic light scatteringdepicted in Fig. 1a:ˆ α ( ω l ) = (cid:88) A (cid:104) i | ˆ H e − p | A (cid:105)(cid:104) A | ˆ H e − p | f (cid:105) ¯ hω l − (cid:15) A − iγ e − p , (3)where the sum goes over all intermediate virtual statewith energy (cid:15) A and wavefunction | A (cid:105) , ˆ H e − p is theelectron-photon interaction and γ e − p is the electroniclinewidth. The initial | i (cid:105) and final | f (cid:105) states are con-sidered to be the electronic and vibrational ground statewith zero energy.To obtain the formal expression of Raman intensitieswe need to calculate the derivative of this quantity withrespect to the normal modes. By displacing atoms ac-cording to normal modes, one can see that energies areonly perturbed in the second order, whilst the wavefunc- tions are already perturbed in the first order. The per-turbed state | A (cid:48) (cid:105) can be expressed as | A (cid:48) (cid:105) = | A (cid:105) + (cid:88) B (cid:68) B (cid:12)(cid:12)(cid:12) ∂V e − I ∂Q ν Q ν (cid:12)(cid:12)(cid:12) A (cid:69) (cid:15) A − (cid:15) B − iγ e − ph | B (cid:105) , (4)where H e − ph = ∂V e − I ∂Q ν is the electron-phonon interaction(i.e. the derivative of the electron-ion potential with re-spect to phonon normal modes) and γ e − ph is the electron-phonon linewidth. Finally, the derivative of the polariz-ability can be obtained by substituting Eq. (4) into Eq.(3): ∂ ˆ α ( ω l ) ∂Q ν = (cid:88) A,B (cid:104) i | ˆ H e − p | A (cid:105)(cid:104) A | ˆ H e − ph | B (cid:105)(cid:104) B | ˆ H e − p | f (cid:105) (¯ hω l − (cid:15) A − iγ e − p )( (cid:15) A − (cid:15) B − iγ e − ph ) . (5)For comparison one can also derive the expression of Ra-man scattering from time dependent perturbation theoryas depicted in the Feynman diagram in Fig. 1b: K = (cid:88) A,B (cid:104) i | ˆ H e − p | A (cid:105)(cid:104) A | ˆ H e − ph | B (cid:105)(cid:104) B | ˆ H e − p | f (cid:105) (¯ hω l − (cid:15) A − iγ e − p )(¯ hω l − (cid:15) B − ¯ hω ν − iγ e − ph ) . (6)The main difference between Raman intensities in ex-pression (5) and (6) manifests in their denominators. Inexpression (6) excitation energy dependence appears inboth energy denominators. This means that by tuningthe laser excitation energy, the Raman peak intensitywill have two maxima, at (cid:15) A and (cid:15) B + ¯ hω ν , which arethe incoming and outgoing resonance, respectively. Theschematic representation of the incoming and outgoingresonances in the integrated Raman intensities can beseen in Fig. 1c.In the expression (5) describing the frequency-dependent Placzek approximation only one of the de-nominators contains the excitation energy. This impliesthat only the incoming resonance will be found in theRaman spectra at (cid:15) A . Nevertheless, the amplitude ofthis resonance is approximately correct, because if the¯ hω l ≈ (cid:15) A condition is satisfied, the second denominatorcan be written as ¯ hω l − (cid:15) B , thus expression (5) and (6)are approximately equivalent. The outgoing resonance,however, is not present in the approximate formula (5),therefore some differences can still be expected in theexcitation profile as depicted in Fig. 1c.In practice, inclusion of the outgoing resonance in thesecond denominator is not possible within the frequency-dependent Placzek approximation, as the derivative iscalculated after the sum over the virtual | A (cid:105) states isperformed in Eq. (3). Therefore, the proper treatmentof both energy denominators would require the calcula-tion of both electron-photon and electron-phonon matrixelements, which is currently unavailable in most DFTcodes. However, our method can be applied on top ofjust about any DFT software and is far more affordable Frequency Raman intensity Symmetry
Frequency-dependent ab initio
Placzek StaticPlaczek
Experiment
Frequency-dependent ab initio
Placzek StaticPlaczek
ExperimentSiC (fcc) 970 977 960 T (LO)797 801 784 T (TO)ZnO (hcp) 448 444 439 E
415 414 412 E
397 390 379 A SiO (quartz) 455 458 465 A
255 259 263 E225 222 207 A
128 131 129 EAnatase (bcc) 646 680 640 E g
505 518 515 A g
503 516 515 B g
371 373 396 B g
152 137 147 E g Monolayer MoS
400 407 406 A (cid:48)
374 389 382 E (cid:48) Black Phosphorene 453 453 471 A g
434 433 440 B g
364 364 363 A g Blue Phosphorene 610 550 N/A 1.00 1.00 N/A A (cid:48)
439 439 N/A 0.92 0.43 N/A E (cid:48)
Armchair GraphiticNanoribbon (N=6) A A A A
457 539 451 A TABLE I. Comparison of frequencies (cm − ) and normalized Raman intensities of experimentally observable vibrational modesbetween different computational approaches and experiments. than a full time dependent perturbation theory calcula-tion of the same, while delivering very good accuracy forthe incoming resonances as we show below.The frequency-dependent polarizability tensor wasevaluated using the built-in linear response algorithms of VASP, within the local density approximation of den-sity functional theory. Raman intensities were calculatedwith our own Python code. Our code displaces atomsaccording to phonon normal modes symmetrically forboth positive and negative directions and calls VASP tocalculate the frequency-dependent dielectric tensor in thedisplaced geometries. After numerical differentiation ourcode computes the Raman spectra according to Eq. (2).Convergence test of this method and technical details ofthe calculations can be found in the Supplementary Ma-terial in sections S1 and S2.To demonstrate the versatility of our method weconsidered several crystals with different lattice struc-tures: face centered cubic (SiC), hexagonal close packed(ZnO), quartz (SiO ), body centered cubic (anatase),two-dimensional (MoS , black and blue phosphorene)and quasi one-dimensional structures (armchair graphiticnanoribbon). We calculated vibrational frequencies and Raman intensities both with the Frequency-dependent ab initio Placzek approximation (using VASP) and thestatic Placzek approximation (using QE) as shown in Ta-ble I. Experimental data of frequencies and Raman in-tensities are also shown in Table I. The irreducible rep-resentations corresponding to each normal mode are alsonoted in the last column.Comparing vibrational frequencies calculated by thetwo DFT codes, generally a good agreement can befound. Note, that since several vibrations are forbiddenby symmetry in the Raman spectra, experimental valuesfor the frequencies are limited to modes with measurableRaman intensity.The overall good agreement of vibrational frequenciesbetween theory and experiments does not propagate intothe Raman intensities as presented in Table I. Since ab-solute Raman intensities are usually difficult to compare,both theoretical and experimental peak intensities werenormalized to the intensity of the highest peak. Com-paring the results of the static Placzek approximationwith experiments one can see that this approach mostlypredicts which vibrational mode will have the highest in-tensity, but relative intensity ratios of smaller peaks arenot accurate for most materials.One possible explanation of the inaccuracy could beattributed to the fact that polarizations of incident andscattered light have a major effect on intensity ratios. Toexclude this effect we took experimental data recorded onpowdered samples containing multiple grains with vari-ous crystallographic orientation. In the case of the two-dimensional systems we compare to Raman spectra mea-sured with unpolarized light. During the theoretical cal-culations, the Raman intensities were averaged over alldirections for bulk crystals and parallel to the plane ofthe crystal for two-dimensional systems.Apart from polarization, Raman intensities are funda-mentally dependent on the excitation energy, as differentelectron-hole pairs are excited at different laser energies.These electron-hole pairs are coupled to vibrations withvarious amplitudes, leading to different peak intensities.In off-resonance conditions (i.e. when the optical gapis large compared to the excitation energy), the Ramanspectra can be modeled in the static limit. MeasurableRaman signal can be detected due to Raman scatteringthrough virtual electron-hole pairs with which resonancescannot occur, however, the absolute peak intensities areseveral orders of magnitude smaller compared to excita-tion energies where the resonance condition can be ful-filled. In the investigated systems shown in Table I thegap is usually comparable to the energy of visible light,thus the difference in the intensity ratio between exper-iments and the static Placzek approximation can be at-tributed to the method neglecting the excitation energydependence.During our calculations using the Frequency-dependent ab initio
Placzek approximation, thepolarizability can be obtained as a function of excitationenergy within the whole visible spectrum. Using thepreviously described method, we calculated the Ramanintensities as shown in Table I and in Fig. 2. Comparingthese results with experimental values, very good agree-ment can be seen. Some small differences can still befound between the calculated and measured intensities,e.g. in the case of black phosphorene the intensity ratioof the two A g modes does not perfectly reproduce theexperimentally observed values. These minor differencescan be attributed to the inaccuracies of the applied abinitio methods.Treating the exchange and correlation energy usingthe local density approximation (LDA) always results inunderestimated gap values, resulting in inaccurate reso-nance positions. A relatively simple method to correctthis is the application of a scissor correction, that is,stretching the gap to the experimental value, but thiswould not necessarily alter the relative intensities of dif-ferent vibrations. Alternatively, the electronic structurecan be improved by taking into account many body cor-rections using the GW method.An additional source of inaccuracy is that the polariz- ability is calculated within the independent particle ap-proximation (IPA), which excludes excitonic effects. Thepolarizability can be calculated more accurately using theBethe-Salpeter equation (BSE), which is typically per-formed on top of a many-body GW calculation. In prac-tice, however, the positions of electronic resonances areusually reproduced within margin of error by treatingthe polarizability at the LDA+IPA level, even for ma-terials with large exciton binding energies. This seem-ingly contradictory behavior is the result of the cancel-lation of two errors, as the difference between the LDAand GW quasiparticle gap usually matches the excitonbinding energy. As a result, the peak energies in theoptical absorption can be approximately reproduced onthe LDA+IPA level, whilst the many-body correctionsonly change the amplitude of these peaks. In Ramanspectroscopy the error introduced by the IPA is expectedto be even less significant, as recent works show neg-ligible difference between Raman spectra calculated onthe LDA+IPA and GW+BSE level. Note, finally, thatthe theory behind the method we have presented in thiswork does not assume that the polarizability is calcu-lated on the IPA level. Therefore, many-body effectscan be included in our approach by replacing the IPApolarizability with the solution of the BSE. While thisupgrade to our method presents an extremely high com-putational challenge, theoretical prediction of Raman in-tensities taking many-body effects into account shouldbecome feasible in the near future as high performancecomputing facilities improve.In conclusion, we presented a hybrid classical-quantummodel of resonance Raman spectroscopy by applying thefrequency-dependent Placzek approximation to ab initio quantum theory. We showed that this approach providesvery good agreement with measurements for the relativeintensities of the Raman modes. While the method islimited to describing incoming resonances, it is more af-fordable than a full time dependent perturbation theorycalculation would be, and can be readily applied usingany DFT code that has the built-in functionality of calcu-lating the frequency-dependent polarizability matrix andthe vibrational modes.Support from the Hungarian National Research, De-velopment and Innovation Office (NKFIH, Grant No.K-115608) is acknowledged. We acknowledge [NIIF]for awarding us access to resource based in Hungaryat Debrecen. This research was supported by the Na-tional Research Development and Innovation Office ofHungary within the Quantum Technology National Ex-cellence Program (Project No. 2017-1.2.1-NKP-2017-00001). This work was completed in the ELTE Excel-lence Program (1783-3/2018/FEKUTSTRAT) supportedby the Hungarian Ministry of Human Capacities. G. K.acknowledges support from the New National ExcellenceProgram (UNKP) of the Ministry of Human Capacities inHungary. V.Z. acknowledges support from the Graphene T E E A A E A E E g A B B E g A ' E' A g B g A g A E g A A A A A R e l a t i v e i n t en s i t y r a t i o Static Placzek Frequency-dependent ab initio
Placzek ExperimentSiC ZnO SiO Anatase MoS BlackPhosphorene BluePhosphorene Armchair GraphiticNanoribbon (N=6)
FIG. 2. Relative intenstiy ratio (compared to the peak with highest intensity) of peaks belonging to different materials andirreducible representations, results are equivalent of those shown in Table I.
Flagship Project and the Computational Shared Facilityat the University of Manchester. J. K. acknowledges theBolyai and Bolyai+ program of the Hungarian Academyof Sciences. [1] C. V. Raman and K. S. Krishnan, Nature , 501(1928).[2] B. R. Carvalho, L. M. Malard, J. M. Alves, C. Fantini,and M. A. Pimenta, Physical Review Letters , 136403(2015).[3] G. Kukucska and J. Koltai, Physica Status Solidi (b) , 00184 (2017).[4] L. G. Canado, A. Jorio, E. H. M. Ferreira, F. Stavale,C. A. Achete, R. B. Capaz, M. V. O. Moutinho, A. Lom-bardo, T. S. Kulmala, and A. C. Ferrari, Nano Letters , 3190 (2011).[5] P. Venezuela, M. Lazzeri, and F. Mauri, Physical ReviewB , 035433 (2011).[6] G. Kukucska, V. Z´olyomi, and J. Koltai, The Journal ofPhysical Chemistry C , 1995 (2019).[7] D. Porezag and M. R. Pederson, Physical Review B ,7830 (1996).[8] M. Lazzeri and F. Mauri, Physical Review Letters ,036401 (2003).[9] G. Placzek, Z. Phys. , 585 (1929).[10] G. Placzek, Z. Phys. , 84 (1931).[11] G. Kukucska, V. Z´olyomi, and J. Koltai, Physical ReviewB , 075437 (2018).[12] H. P. C. Miranda, S. Reichardt, G. Froehlicher,A. Molina-Snchez, S. Berciaud, and L. Wirtz, Nano Let-ters , 2381 (2017).[13] Y. Gillet, S. Kontur, M. Giantomassi, C. Draxl, andX. Gonze, Scientific Reports , 7344 (2017).[14] Y. Wang, B. R. Carvalho, and V. H. Crespi, PhysicalReview B , 161405(R) (2018). [15] E. E. Salpeter and H. A. Bethe, Physical Review , 1232(1951).[16] G. Kresse and J. Furthm¨uller, Physical Review B ,11169 (1996).[17] G. Kresse and D. Joubert, Physical Review B , 1758(1999).[18] P. Giannozzi, S. Baroni, N. Bonini, M. Calandra, R. Car,C. Cavazzoni, D. Ceresoli, G. L. Chiarotti, M. Cococ-cioni, I. Dabo, A. D. Corso, S. d. Gironcoli, S. Fabris,G. Fratesi, R. Gebauer, U. Gerstmann, C. Gougoussis,A. Kokalj, M. Lazzeri, L. Martin-Samos, N. Marzari,F. Mauri, R. Mazzarello, S. Paolini, A. Pasquarello,L. Paulatto, C. Sbraccia, S. Scandolo, G. Sclauzero, A. P.Seitsonen, A. Smogunov, P. Umari, and R. M. Wentzcov-itch, Journal of Physics: Condensed Matter , 395502(2009).[19] U. Balachandran and N. G. Eror, Journal of Solid StateChemistry , 276 (1982).[20] P. Gillet, A. L. Clac’h, and M. Madon, Journal of Geo-physical Research: Solid Earth , 21635 (1990).[21] H. Harima, Microelectronic Engineering The SymposiumK Proceedings of the 3rd International Conference onMaterials for Advanced Technologies (ICMAT 2005), ,126 (2006).[22] H. Li, Q. Zhang, C. C. R. Yap, B. K. Tay, T. H. T. Ed-win, A. Olivier, and D. Baillargeat, Advanced FunctionalMaterials , 1385 (2012).[23] T. Sander, S. Eisermann, B. K. Meyer, and P. J. Klar,Physical Review B , 165208 (2012).[24] W. Lu, H. Nan, J. Hong, Y. Chen, C. Zhu, Z. Liang,X. Ma, Z. Ni, C. Jin, and Z. Zhang, Nano Research ,853 (2014).[25] H. Kuzmany, L. Shi, S. Cambr´e, M. Martinati,W. Wenseleers, J. K¨urti, J. Koltai, G. Kukucska, K. Cao,U. Kaiser, T. Saito, and T. Pichler, “Unpublished,”(2019).[26] M. Gajdos, K. Hummer, G. Kresse, J. Furthm¨uller, andF. Bechstedt, Physical Review B , 045112 (2006).[27] G. Kukucska, V. Z´olyomi, and J. Koltai, “RamPy pack- age,” https://github.com/gkukucska/RamPyhttps://github.com/gkukucska/RamPy