Calculation of the Raman G peak intensity in monolayer graphene: role of Ward identities

The absolute integrated intensity of the single-phonon Raman peak at 1580 cm^{-1} is calculated for a clean graphene monolayer. The resulting intensity is determined by the trigonal warping of the electronic bands and the anisotropy of the electron-phonon coupling, and is proportional to the second power of the excitation frequency. The main contribution to the process comes from the intermediate electron-hole states with typical energies of the order of the excitation frequency, contrary to what has been reported earlier. This occurs because of strong cancellations between different terms of the perturbation theory, analogous to Ward identities in quantum electrodynamics.


Introduction
In the past decades, Raman spectroscopy [1] techniques were successfully applied to carbon compounds, such as graphite [2] and carbon nanotubes [3]. Upon the discovery of graphene [4], Raman spectroscopy has proven to be a powerful and non-destructive tool to identify the number of layers, doping, disorder, strain, and to characterize the phonons and electron-phonon coupling [5,6,7,8,9,10,11].
The most robust feature in the Raman spectra of all conjugate carbon compounds is the so-called G peak, corresponding to the in-plane bond-stretching optical vibration of sp 2 -hybridized carbon atoms. In graphene this vibration has the E 2g symmetry, it is doubly degenerate, and its frequency ω ph ≈ 1580 cm −1 . While the G peak frequency is sensitive to external factors, such as doping level [6,7,8] or strain [12,10,11], its total frequency-integrated intensity I G is often assumed to remain constant under the change of many external parameters, depending only on the excitation frequency. Thanks to this robustness, I G is often used as a reference to which intensities of other peaks are compared [13,14,15], since measurement of absolute peak intensities represents a hard experimental task (the first use of I G as a reference for other Raman peak intensities in graphite probably dates back to 1970 [16]).
Given this popular role of I G as a reference, it is clear that having at hand a theoretical expression of the absolute intensity in terms of the basic parameters of the material (such as the electronic dispersion, the electron-phonon coupling, etc.), would be useful. Nevertheless, to the best of the author's knowledge, no such expression is available at present. I G was discussed in a recent paper by the author [17], where the conclusion was reached that even in the limit when the excitation frequency ω in is small compared to the energy scale t 0 characterizing the electronic dispersion (several eV), the intensity I G is contributed by the whole conduction and valence bands, not just by low-energy states in the vicinities of the Dirac points. As will be shown below, this conclusion is wrong, since Ref. [17] missed strong cancellations occurring as a consequence of a Ward identity. This affects the dependence of I G on ω in .
In the present work I G is calculated for a clean graphene monolayer suspended in vacuum. When ω in ≪ t 0 , the intensity is indeed determined only by electronic states with energies ∼ ω in , not the whole conduction and valence bands (however, these energies do not have to be close to ω in /2). In this limit I G can be expressed in terms of a few parameters, characterizing these low-energy states. As correctly mentioned in Ref. [17], trigonal warping of the electronic bands and the anisotropy of the electronhole coupling turn out to be crucial for the G peak. In the limit ω ph ≪ ω in ≪ t 0 the frequency dependence is I G ∝ ω 2 in . At higher frequencies, the calculation is performed using the nearest-neighbor tight-binding model. The main feature found is a strong enhancement of I G when the frequency matches the energy of electron-hole separation at the M point of the electronic first Brillouin zone, where a van Hove singularity in the electronic density of states occurs.

Calculation
The Bloch form of the electronic wave function in the tight-binding model is where ψ a (r) is the atomic wave function localized near the origin, α = A, B labels the sublattices, and the two-dimensional integer vector n labels the unit cells, each containing two atoms from the two sublattices. The wave vector k runs over the first Brillouin zone. The Bloch amplitudes U αk satisfy the Shrödinger equation (we neglect the non-orthogonality of the atomic orbitals) where the Hamiltonian matrix is given bŷ Here t 0 ≈ 3 eV is the nearest-neighbor coupling matrix element, and d 1,2,3 are the vectors connecting an A atom to its three nearest neighbors, as shown on Fig. 1. Their length is |d 1,2,3 | = a ≈ 1.42Å. The energy eigenvalues are ǫ = ±ǫ k , where ǫ k ≡ |H k |. Note that H k is not a periodic function in the first Brillouin zone; its values on the boundaries related by a reciprocal lattice vector differ by a constant phase factor e ±2πi/3 . Suppose that all atoms belonging to the sublattice α are displaced by u α , and u A = −u B . Such a uniform displacement corresponds to optical phonons with the wave vector q = 0. These are the phonons responsible for the G peak, as fixed by the momentum conservation (the wave vectors of both incident and scattered photons are negligibly small). In the tight-binding model the mechanism of the electron-phonon coupling is the change of the nearest-neighbor electronic matrix element t 0 due to the change in the bond length a, ∂t 0 /∂a ≈ −6 eV/Å. Then, to linear order in the atomic displacements, the electron-phonon interaction Hamiltonian is given bŷ The electromagnetic field of the incident and scattered light is described by the vector potential A. The Hamiltonian describing interaction of the electrons with the field can be obtained by the Peierls substitution k → k − (e/c)A in the electronic Hamiltonian. Since we are interested in a two-photon process, we should expand the Hamiltonian to the second order in A: Here i, j, l = x, y label the in-plane Cartesian components, and the summation over repeated indices is assumed. Note that for the description of resonant multiphoton processes it is sufficient to keepĤ (1) k term only, as it was done in [17], since the contribution of other terms is (i) off-resonant and (ii) smaller by a factor ∼ ω in /t 0 . However, for the one-phonon process studied here the contributions from all terms turn out to be of the same order; in fact, they almost cancel each other.
To calculate the Raman scattering probability, we follow the general steps of Ref. [17]: pass to the second quantization, calculate the scattering matrix element, and substitute it into the Fermi Golden Rule. The zero-temperature electronic Green's function is a 2 × 2 matrix, Upon quantization of the phonon field the atomic displacements corresponding to the optical phonons are expressed in terms of the phonon creation and annihilation operatorsb † q,µ ,b q,µ aŝ Here µ = x, y labels the two degenerate phonon modes, e (µ) is the unit vector in the corresponding direction, M is the mass of the carbon atom, N is the total number of the carbon atoms in the crystal, and "h.c." stands for the hermitian conjugate. A small finite wave vector q has to be introduced in order to treat the in-plane momentum conservation properly; only the q = 0 modes will contribute to the final result. The transverse electromagnetic field is quantized inside the large volume L x L y L z in the Coulomb gauge:Â Here Q is the three-dimensional wave vector, ℓ = 1, 2 labels the two transverse polarizations along the two unit vectors e (Q,ℓ) ⊥ Q.
The probability for an incident photon with the polarization e in and frequency ω in to be scattered into an element of the solid angle do out with the polarization e out , and with the frequency ω out = ω in − ω ph fixed by the energy conservation, is given by Here M ijl is the transition matrix element (the factor of 2 explicitly takes care of the two spin projections of the electron). It is given by the sum of the following contributions: The factor e iǫ0 + prescribes closing of the ǫ-integration contour in the upper halfplane. In fact, it is important only for the last term, Eq. (9f) which corresponds to the sum over the filled valence band. Upon integration over ǫ, Eq. In fact, D 1k =D 1k because the Green's function and all the vertex matrices satisfŷ G T k (ǫ) =σ xĜk (ǫ)σ x andĜ k (−ǫ) = −σ zĜk (ǫ)σ z (the T superscript denotes the matrix transpose, andσ x,z are the Pauli matrices).
If one integrates each of the terms in Eqs. (9a)-(9f) separately, the main contribution to the k-integral comes from the "bulk" of the first Brillouin zone, rather than the vicinities of the Dirac points (as it was pointed out in Ref. [17]). Let us formally consider, however, the whole expression for M ijl at ω in = ω out = 0. Then the sum of the terms (9a)-(9f) is a total derivative, ∂ 2 [F l kĜ k (ǫ)]/(∂k i ∂k j ). Thus, the k-integral vanishes, i. e., at ω in = ω out = 0 all terms contributing to the matrix element, cancel. This cancellation is far more general than the tight-binding model, used here. Indeed, consider linear response of the electronic current j to a static homogeneous vector potential A. The corresponding response function χ ij must vanish, since an observable (current) cannot respond to a pure gauge ‡. This holds for any configuration of the atomic positions, hence its derivative with respect to ‡ Strictly speaking, it is the limit lim q→0 lim ω→0 χ ij (q, ω) which vanishes, while the calculation done in this work corresponds to the opposite order of limits. The two limits commute for undoped graphene, when the Fermi surface has zero length. At finite doping, the contribution of the filled valence band states far below the Fermi energy still cancels out. the displacements, ∂χ ij /∂u l , must also vanish. As follows from the Kubo formula, the matrix element M ijl at ω in = ω out = 0 is equal just to this derivative, up to a constant factor. This fact is analogous to the Ward identity in the quantum electrodynamics. The cancellation of different terms contributing to M ijl was overlooked in Ref. [17], and it was wrongly concluded that M ijl is determined by the whole of the first Brillouin zone. In fact, for ω in ≪ t 0 it is sufficient to focus on the vicinities of the two Dirac points, as will be seen below.
The rest of the calculation is tedious, but straightforward. Let us focus on the limit ω in ≪ t 0 first. In the vicinity of the K point we write k = K + p and expand: The expansion around the K ′ point can be obtained by flipping e x → −e x , p x → −p x . Note that going one order beyond the Dirac approximation is necessary, since the leading contribution to M ijl vanishes. (Indeed, the Dirac Hamiltonian has a continuous rotation symmetry, so any third-rank tensor must vanish; the continuous rotation symmetry is lifted by the trigonal warping term.) In the nearest-neighbor tight-binding model one has Still, in Eqs. (10a), (10b) we prefer to keep all four parameters v, α 3 , F 0 , F 1 independent, since going beyond the nearest-neighbor tight-binding model would invalidate Eq. (11), while the form of Eqs. (10a), (10b) is fixed by the symmetry §. It is instructive to write down explicitly the expression for the matrix element after the integration over the directions of p (we change the integration variable from p § In fact, the C 6v crystal symmetry allows a term α 0 (p 2 x + p 2 y ) in the expansion of H k (electron-hole asymmetry), which appears already in the second-nearest-neighbor approximation. However, it does not contribute to the matrix element calculated here because of its full rotational symmetry.
to ǫ k and denote by A C = L x L y /N = √ 27a 4 /4 the area per carbon atom): The integral over ǫ k corresponds to the sum over the intermediate states whose energies are 2ǫ k (an electron with the energy −ǫ k in the valence band is promoted to the conduction band where its energy is ǫ k ). The poles are bypassed by adding an infinitesimal imaginary part ǫ k → ǫ k − i0 + , as follows from Eq. (4). At frequencies ω in not small compared to t 0 the k integral in M xxy should be evaluated numerically. We focus in the limit ω ph ≪ ω in . In this limit the overall Raman efficiency, i. e., the absolute probability for an incident photon with the polarization e in and frequency ω in to be scattered in the full solid angle 4π with any polarization, accompanied by the emission of a 1580 cm −1 optical phonon, is given by The dimensionless electron-phonon coupling constant is defined as λ Γ = ( √ 27/M ω ph )(t −1 0 ∂t 0 /∂a) 2 , which coincides with the definition of Ref. [17] in the tightbinding model. The dimensionless function f (ω in /t 0 ) is plotted in Fig. 3. In the low-energy limit, ω in ≪ t 0 , its value is f (0) = 1, as follows from Eq. (12a). When the frequency matches the van Hove singularity at the M point of the first Brillouin zone, ω in ≈ 2t 0 , the intensity diverges, f (ω in /t 0 ) ∝ 1/(ω in /t 0 − 2) 2 . This divergence is cut off at the scale |ω in − 2t 0 | ∼ max{ω ph , γ}. However, in the absence of a reliable information on electronic relaxation processes in this region of the spectrum, we prefer not to study this issue in detail, leaving the divergence as it is.

Discussion
First, let us see which electron-hole states contribute the most to the G peak intensity. In principle, for ω in ≪ t 0 one can imagine three situations. They are illustrated Fig. 4 where the π electron dispersion is shown schematically, together with the contributing states. In the situation (a) these are the states whose energies are close to the half of the excitation frequency, the difference ǫ k −ω in /2 being small by some parameter. This was shown to be the case for the 2D peak, the second-order overtone of the D peak at 2700 cm −1 , which is fully resonant and |ǫ k − ω in /2| ∼ γ, where 2γ ≪ ω in is the electron inelastic scattering rate [18,17]. This is also the case for the doubly-resonant D peak at 1350 cm −1 , where |ǫ k − ω in /2| ∼ ω ph ≪ ω in [19,20]. In the situation (b) the difference |ǫ k − ω in /2| ∼ ω in itself, and no small parameter enters. Finally, in the situation (c) the Raman process is contributed by the whole Brillouin zone, ǫ k ∼ t 0 . A natural guess for the G peak would be the option (a): indeed, there is a contribution from intermediate states that violate energy conservation by the energy ∼ ω ph . However, in Ref. [17] it was found that the main contribution to the integral in Eq. (8) with only two terms (9a), (9b) came from the states ǫ k ∼ t 0 [option (c)]. In the present work we have seen that this contribution, in fact, is canceled by other terms, missed in Ref. [17], due to the Ward identity. To see that the correct picture for the G peak in fact corresponds to Fig. 4(c), let us consider doped sample, i. e., with the Fermi energy ǫ F detuned from the Dirac point. Then transitions involving states with ǫ k < |ǫ F | are simply blocked by the Pauli principle, so the lower limit of the integral in Eq. (12a) must be set to |ǫ F |. Let us analyze the first term in Eq. (12a) (the second term is almost equal to the first, the third one is smaller by a factor ∼ ω 2 ph /ω 2 in ): The value of the integral is determined entirely by the ratio |ǫ F |/ω in , and no other energy scales enter, whether small (ω ph ) or large (t 0 ). That is, the value of the integral in Eq. (12a) is determined not just by the vicinity of the pole, but by the whole range of energies from ǫ k = 0 to ǫ k ∼ ω in . Numerically, when ǫ F is raised from zero to ω in /4 (half of the distance to the pole), |M xxy | 2 is changed by about 12%. In other words, the uncertainty in the energy of the electron-hole states, contributing to the process, is of the order of their energy itself (∼ ω in ). By virtue of the energy-time uncertainty principle, the duration of the process (the typical lifetime of the virtual electron-hole pair) is ∼ 1/ω in . The relevant length scale is thus v/ω in ∼ 3Å for v = 10 8 cm/s ≈ 7 eV ·Å and ω in = 2 eV. The Raman process giving rise to the G peak is thus extremely local in space, involving just a few carbon atoms, as it was noted earlier [21]. Its locality can be also understood from the quasiclassical Figure 4. Schematic view of the π electron band structure (solid curves), with the states effectively contributing to a Raman process (pink areas)for ω in ≪ t 0 : (a) states within the narrow region |ǫ k −ω in /2| ≪ ω in ; (b) a broad region of states with |ǫ k − ω in /2| ∼ ω in , but still ǫ k ≪ t 0 ; (c) the whole Brillouin zone, ǫ k ∼ t 0 .
real-space picture of Raman scattering in graphene [20]. This short length scale should be contrasted to much longer scales responsible for Raman peaks produced according to Fig. 4 (a), for example, the doubly-resonant D peak at 1350 cm −1 [23,24,22,20].
Another consequence of the cancellation of the contribution from ǫ k ∼ t 0 is the frequency dependence of I G . The standard textbook dependence I ∝ ω 4 (for simplicity we assume ω in ≈ ω out ) is obtained for systems whose excited states have energies much higher than ω in [25]. The ω 4 dependence, suggested in Ref. [17], was essentially of the same origin: it was assumed that the main contribution to the process came from states in the bulk of the Brillouin zone, and thus with high energies. As we have seen, this conclusion is wrong, and the relevant states have energies ǫ k ∼ ω in . This remains true even at low frequencies, as the electronic spectrum in graphene has no (or almost no) gap. As a result, at low frequencies the dependence is ∝ ω 2 in , as seen from Eq. (13). Nevertheless, a recent measurement of Raman peak intensities performed on graphite nanocrystallites gives a ω 4 in dependence [15,26]. For comparison, we plot I G from Eq. (13) together with two curves ∝ ω 4 in with two different coefficients in Fig. 5. Let us compare I G to the intensities of other peaks which have been calculated theoretically and measured experimentally. The D peak at 1350 cm −1 is due to emission of the scalar A 1 phonons with wave vectors near the K point. Due to momentum conservation, the D peak is activated by defects or edges. For a defect-free graphene flake of the size L a its intensity is given by [20] Here ω D ≈ 0.17 eV and λ K are the frequency and the dimensionless electron-phonon coupling constant for the corresponding phonons, v ≈ 7 eV ·Å is the electron velocity (the slope of the Dirac cones), 2γ is the electron inelastic scattering rate, L x L y is the area of the excitation laser spot, and C is a numerical coefficient which depends on the shape of the flake and the character of the edges. We take λ Γ = 0.03, as obtained from the doping dependence of the G peak frequency [6,7], and λ K /λ Γ = 3 as extracted from the ratio of intensities of the intensities of the 2D peak at 2700 cm −1 and the 2D ′ peak at 3250 cm −1 [17]. Let us assume the inelastic scattering rate to be dominated by electron-phonon scattering (which is a lower bound), then it can be estimated as γ = (λ Γ + λ K )ω in /8 [18]. For a small flake, Eq. (13) should be weighted by the ratio of the area of the flake to the excitation spot area. If we take an ideal hexagonal flake with armchair edges (which definitely represents an upper bound for I D ), then its area is L 2 a √ 27/8 and C = 4. For a round flake of the diameter L a with atomically rough edges the area is πL 2 a /4 and C = π 2 /18, i. e., I D is about 10 times smaller than for the ideal flake. Taking ω in = 2eV and using the low-frequency limit of Eq. (13), f (ω in /t 0 ) = 1, we obtain I D /I G ≈ (300 − 3000 nm)/L a . The experimentally found intensity is I D /I G = (560 nm · eV 4 )/(L a ω 4 in ) [15,26], which at ω in = 2 eV gives an almost 10 times smaller value than the calculated one for a flake with rough edges.
We can also consider 2D peak at 2700 cm −1 , whose intensity was calculated to be [18,17] For the same values as above we obtain I 2D /I G ≈ 150. The experimental value for graphene samples on a substrate is typically I 2D /I G ≈ 5 [5,12]. For suspended samples the ratio can be I 2D /I G ≈ 9 [27] and even I 2D /I G ≈ 17 [28]. Still, the experimental ratio is about 10 times smaller than the calculated one. Thus, the theory overestimates both I 2D /I G and I D /I G by about an order of magnitude. It seems more reasonable to assume that Eq. (13) for I G should be blamed for this discrepancy, rather than expressions for I D and I 2D . Indeed, as discussed above, I D and I 2D are determined by just a few material parameters (electron-phonon coupling constants, electronic velocity, etc.) which have been calculated by different methods and measured in many independent experiments. Moreover, corrections to Eq. (16) due to the trigonal warping and electron-hole asymmetry have been estimated in Ref. [17] and turned out to be small. At the same time, the parameters determining I G (electronic trigonal warping, electron-phonon coupling anisotropy) are not really known, the validity of the tight-binding model at sufficiently high energies is questionable, and even the calculated frequency dependence of I G does not seem to agree with the experiment.
Still, what is the main source of such a strong discrepancy? Is it possible that ω in = 2 eV is already sufficiently close to the singularity at ω in = 2t 0 , so that the dependence of I G on ω in is noticeably steeper than ω 2 , and the value of I G is higher than that predicted by the low-frequency asymptotics? On the one hand, the nearestneighbor tight-binding model with t 0 ≈ 3 eV seems to describe the band structure reasonably well in the energy range of interest, as is seen from its comparison with the results obtained by angle-resolved photoemission spectroscopy [29]. Then the energy of electron-hole separation at the M point is 2t 0 ≈ 6 eV and ω in /t 0 = 0.6 is quite close to the low-energy limit, as seen from Fig. 3. On the other hand, a recent ab initio calculation of the band structure of graphite gives the electron-hole separation at the M point about 4 eV [30]. Also, the anisotropy of the electron-phonon coupling may be stronger than that obtained from the nearest-neighbor tight-binding model [31]. In addition, the spectrum of electron-hole excitations can be modified by excitonic effects: the Coulomb attraction between the electron and the hole lowers the energy of the excitation with respect to its non-interacting value, so that the singularity in I D should be at a lower frequency. The author is aware of only one study of excitonic effects in the optical response of monolayer graphene; within the Dirac model it was shown that excitonic effects result in a weak feature at the border of the electron-hole continuum [32]. Their importance for optical excitation of electrons near the M point of the first Brillouin zone remains to be studied. Also, a direct measurement of the energy corresponding to the singularity by an excitation in the ultraviolet frequency range would shed light on this issue.

Conclusions
We have calculated the frequency-integrated intensity of the Raman G peak in monolayer graphene using the nearest-neighbor tight-binding model for electrons. The resulting intensity is determined by the trigonal warping of the electronic bands and the anisotropy of the electron-phonon coupling, and is proportional to the second power of the excitation frequency. Comparison to the intensities of other peaks, which have been measured experimentally and calculated theoretically, suggests that the present calculation underestimates I G by about an order of magnitude.