Flexural phonons in supported graphene: from pinning to localization

We identify graphene layer on a disordered substrate as a system where localization of phonons can be observed. Generally, observation of localization for scattering waves is not simple, because the Rayleigh scattering is inversely proportional to a high power of wavelength. The situation is radically different for the out of plane vibrations, so-called flexural phonons, scattered by pinning centers induced by a substrate. In this case, the scattering time for vanishing wave vector tends to a finite limit. One may, therefore, expect that physics of the flexural phonons exhibits features characteristic for electron localization in two dimensions, albeit without complications caused by the electron-electron interactions. We confirm this idea by calculating statistical properties of the Anderson localization of flexural phonons for a model of elastic sheet in the presence of the pinning centers. Finally, we discuss possible manifestations of the flexural phonons, including the localized ones, in the electronic thermal conductance.

Most of the research in graphene emphasizes the relativistic character of its electron spectrum. However, graphene is also interesting due to its out-of-plane (flexural) vibrational phonon modes. Flexural phonons (FPs) are a unique addition that Van der Waals heterostructures have brought into microscopic physics 1 . Usually, FPs are considered in the context of the suspended graphene. Here we argue that graphene layer placed on the top of the supporting SiO 2 substrate gives an opportunity to observe Anderson localization 2 for the FPs.
To get an idea, let us recall the known facts about a long-wave acoustic wave scattering on a cylinder. The output depends drastically on the boundary conditions for the velocity potential Φ on the surface of the cylinder 3 . If the velocity component normal to the surface of the cylinder vanishes, i.e., ∂ Φ = = 0 r r a , the scattering cross-section σ R is proportional to a(ka) 3 , where a is the radius of the cylinder and k is the wave vector. This is the conventional Rayleigh scattering result 4 for two dimensional (2d) geometry. However, when pressure is constant, the boundary condition reads Φ(a) = 0, and this influences forcefully the scattering. Unlike the Rayleigh scattering, the zero angular harmonic is involved, and as a result the cross section diverges at small k as ∝ − ( ) k ln ka 2 1 1 . (The same takes place for an electro-magnetic wave scattering on a metallic cylinder).
In graphene, the substrate cannot scatter effectively the usual acoustic waves, longitudinal and transverse, because graphene itself is one of the most rigid substances. The other thing are the out of plane phonons. From the analysis of intrinsic and extrinsic corrugation of monolayer graphene deposited on SiO 2 substrate, it has been concluded that in this system the layer is suspended between hills of the substrate [5][6][7][8] . We have checked that scattering of the FPs from areas attached to the substrate is similar to the scattering from a rigid obstacle 9 . The zero harmonic is also involved, and the scattering cross-section diverges as σ fl = 4/k.
In this work we study statistical properties of the out of plane excitations for a pinned-suspended flexible sheet. Whether pinning centers are located in the vicinity of the maximal heights of the substrate where the interaction with the layer is the strongest, or there are charges on the substrate which interact strongly with their images, will be not important for our purposes. First of all, we are interested in the scattering rate of the FPs in the presence of the randomly located pinning centers with concentration n i . Taking into consideration that the spectrum of the FPs ω(k) = αk 2 is quadratic, i.e., velocity is linear in k, one obtains a scattering rate τ −1 = vσ fl n i , that is finite in the low-energy limit. Thus, for the FPs one may expect localization of the low-energy modes with Scientific  ω(k) < τ −1 . This is in striking contrast with localization of acoustic modes, which is known to happen only at high enough frequency [10][11][12][13][14] .
Let us comment upon the graphene layer deposited on the top of the corrugated substrate. Naively, the membrane-like layer either follows the substrate or hovers over the surface at some distance. Measurements with the use of cantilevers 15 indicate, however, towards a possibility of the detaching a graphene sheet from a substrate to relieve its strain by slipping. (This is manifest by straightening of the cantilever). In the case of the SiO 2 substrate, both experiment and theory agree that for typical magnitude for corrugations, the graphene layer is partially detached from the substrate. Moreover, the theoretical considerations 16,17 justify the use of a contact force that is finite when graphene is conforming to the substrate and zero otherwise. The basic experimental facts 7 which lead to the conclusion that graphene layer deposited on SiO 2 is partly freely suspended are as follows: The long-range corrugation of the substrate with the correlation length of about 25 nm is also visible on the graphene sheet, but with a smaller amplitude than on the substrate. Mesoscopic corrugations with smaller length of about 15 nm not induced by the substrate were also identified. These short range corrugations are similar in height and wavelength to the ones observed on suspended graphene 18,19 . In addition, the picture of partially suspended graphene, and the presence of the FPs in the graphene on the SiO 2 substrate, has been confirmed by the transport 20 and thermal measurements 21,22 .

Results
Here, we demonstrate that the FPs in a pinned-suspended flexible sheet are more similar to disordered electrons rather than to acoustical phonons. The main point here is that the pinning centers are effectively rigid obstacles for the FPs. As we have already explained, this leads to the non-vanishing scattering rate at small energies. Let us touch upon this point in more detail.
Pinning potential as a barrier for FPs. Usually by a rigid obstacle one understands an inclusion with the Young's modulus much higher than that in surrounding area. For graphene, which itself is very rigid, this is not an issue. However, the pinning potential introduces an energy barrier of a finite height for the flexural modes: Here h(r, t) is displacement in the out of plane direction; the term describing the barrier is on the right-hand side of the above equation. As a result, a flexural phonon with an energy smaller than ω 0 cannot enter the area of pinning. We have checked that such an inclusion is equivalent to a rigid obstacle. For not small ka, the cross-section σ fl (k) = 4f(ka)/k. Interestingly enough, for  ka 1, f(ka) ≈ ka, for a discussion see Section I in METHODS. Therefore, the limiting cross-section is ≈4a, i.e., twice larger than the width of the obstacle 3 . (This form of σ fl (k) is, of course, valid only when the energy of the phonon is much less than the pinning potential, i.e., ω ω  k ( ) 0 . Relying on the existing experimental data 23 we assume that the pinning potential ω 0 is about few meV).
Let us now touch upon a subtle question of the openness of the ensemble of the FPs which can be relevant for their localization. In the discussed model, the low-energy FPs cannot penetrate into the areas of strong contact with the substrate. Therefore, the most relevant channel connecting the FPs with the other degrees of freedom is the interaction with conducting electrons. The effect of this interaction can be estimated by comparing the amount of heat stored in the FPs with the rate of cooling of electrons. At low temperatures, most of the heat in the graphene layer is stored by the FPs, as they are the softest modes. The amount of this heat is ∝T 2 . The rate of energy exchange between electrons and FPs, as discussed below in section "Thermal transport", is proportional to T 3 . Therefore, the openness of the FPs, cannot effectively destroy their localization at low temperatures.
Numerical study. A general question has been addressed: If to compare with the electrons propagating in a disordered lattice, will the statistical properties of the eigemodes of the pinned elastic layer be the same or different? The question makes sense because for phonons in a pinned-suspended sheet there is no analogue of the on-site disordered potential W. Instead, there is concentration of the pinned sites. Furthermore, the FPs are described by the square of the Laplacian, rather than by the Laplacian in the case of electrons. To understand the general properties of the FPs in the presence of random pinning scatterers, we solve the equation of motion for the out-of-plane displacements using finite difference method on a 2d square lattice. We study a model in which the graphene sheet is completely attached at the pinning centers. For that, we use discretized LHS of Eq. 1 with condition h = 0 at randomly chosen pinned sites, so that one pinned site represents an attached area of the size ≈2a (See Fig. 6 in METHODS for illustration of the model).
We estimate the size of an attached area to be  a 2 7 nm. Typical distance between the pinning centers a i is around 20 nm; we will assume that = − n a . Correspondingly, we studied "samples" with 5-20% of the pinned sites to determine statistical properties of the eigenmodes and eigenvalues. In doing so, we considered samples with periodic boundary conditions of the size up to 200 × 200 sites. In what follows, we measure the energy eigenvalues E in units of α/(2a) 2 which approximately equals 0.08 K for the parameters mentioned above. We found out here that typical energy scale for strong localization, ω loc str , for the discussed concentration of pinned sites is a fraction of 1 K. We believe that for graphene layer on the top of the SiO 2 substrate ω ≈ . Phononic "conductance". As is well known, localization is a quantum critical phenomenon. The peculiarity of 2d is that the critical point is at  = g 1/ 0, where  g is electrical conductance per square measured in units e 2 /(2πħ). At a finite g  , statistical properties are determined by the localization length l loc , which is analogue of a correlation length at a quantum phase transition. A sample of size < L l g ( ) loc  is in the regime of criticality which may take place in a very broad range of the sample sizes, because for small  g the localization length is exponentially large 24 . A consequence of strong fluctuations of wave function amplitude in the critical region is the multifractality (that is when an eigenstate is extended but the occupied volume is noticeably smaller than the volume of the sample).
Turning back to disordered FPs, the first question that needs to be answered is: Is there a transition to delocalized states at a certain energy (i.e., the "metal-insulator" transition with the mobility edge), or there is a crossover from strong to weak localization (WL)?
To figure this out, we first studied the dependence of the Inverse Participation Ratio (IPR) on the sample size L for various phonon energies. A discussion of the IPR is given below in METHODS, section "Weak multifractality of eigenfunctions". (For more details the reader is referred to ref. 25 ). From our simulations, it is clear that low-energy modes are localized: the IPR scales with the sample size to a finite value. For higher energies, the behavior of the wave functions changes, see Fig. 1, because the localization length l loc starts to exceed the sample size. We are particularly interested in studying the FPs in this region when  l L loc . Note that although the 2d Anderson model does not constitute a truly critical system, thanks to exponentially large (but still finite) localization length at   g 1, the criticality takes place in a very broad range of the system sizes,  L l loc . Therefore, 2d electrons at large g  share many common properties with systems at the critical point of the metal-insulator transition 2,26 . As we shall see, similar physics holds also for our system of the FPs.
We proceed as follows: From the size dependence of the IPR at a given energy as shown on the Fig. 2, we extracted an energy-dependent fractal dimension. For disordered electronic system of a given symmetry class, the fractal dimension is determined by the conductance g  . For example, in the case of the Gaussian Orthogonal Ensemble (GOE), the size dependence of the IPR is described by the fractal dimension 27,28 equal to on L is due to the WL corrections. Thus, for each concentration of the pinned sites we can prescribe for different energies ω the corresponding value of the phonon "conductance" g ph (%, ω), using the expression for the fractal dimension D 2 for electrons. We defer the discussion of the dependence of D 2 on the sample size, D 2 (L), to the end of section "Weak multifractality of eigenfunctions".  For disordered electrons in 2d, the well developed theory connects the behavior of various physical quantities with the value of the conductance 25 . We have calculated numerically the same quantities for the FPs, using the values of g ph extracted from IPR, and found a very good agreement with the theoretical predictions existing for the disordered electrons in the case of the GOE. Below we present some results of such calculations.
Statistical properties. First, we have checked (see Fig. 7 in METHODS, section "Energy level statistics") that the distribution function of the level spacing P(s) for localized states has almost Poissonian statistics, while for metallic states it is of the Wigner-Dyson form. Next, in 2d it becomes especially interesting to study the variance ω Σ Ω ( , ) 2 , which is a two-level correlation function characterizing the fluctuations of the number of levels N in a strip of width Ω around the energy ω: 2 . The reason why it is of particular interest is that, in contrast to d = 1 and 3, in two dimensions this quantity is directly related to the WL corrections 29 . Figure 3 demonstrates the level number variance as a function of the ratio Ω/Δ, where Δ is the average level spacing. It starts with the ergodic behavior described by the Random Matrix Theory (RMT). The ergodic regime holds up to Ω about the Thouless energy. For larger Ω there is a noticeable deviation: the variance starts to increase rapidly. The numerical results presented in Fig. 3 are in full accord with the theoretical expression obtained by us for d = 2; for details the reader is referred to section "Energy level statistics".
To the best of our knowledge, this is the first demonstration of the mesoscopic fluctuations of the number of levels in d = 2, while for the 3d Anderson model, the function ω Σ Ω ( , ) 2 was studied long ago 30 . Another important statistical property is the distribution of the amplitudes of the eigenmodes, ψ 2 , which is called the wave function intensity distribution  y ( ), where in our case y is ∝h 2 . For metallic granulas, in the ergodic regime described by the RMT, the intensity distribution is given by the Porter-Thomas distribution  y ( ) RMT . Owing to the fluctuations in the diffusive motion, there appear deviations from the ergodic behavior. When calibrated with respect to  y ( ) RMT , the function y ( )  yields a curve with a very specific non-monotonous shape. As Fig. 4 shows, an excellent agreement with the theory of ref. 31 is found.
To summarize, we have calculated numerically a number of quantities characterizing statistical properties of FPs using the values for g ph extracted from the data for the IPR, and found a very good agreement with the theoretical predictions existing for the disordered Anderson model electrons in the case of the Orthogonal Class of Universality. Furthermore, the theoretical expressions for the number of variance ω Σ Ω ( , ) 2 and the wave function intensity y ( )  , both are intimately connected with the effects of the WL originating from the Cooperons. The excellent agreement demonstrated in Figs 3 and 4 justifies that in the discussed model the regime of WL is the same as in the Anderson model in 2d. We believe that the reason for the observed universal behavior is that the FPs in the lattice with pinned sites are eventually described with the same Non-Linear σ-model as disordered electrons in the Orthogonal Class of Universality.
Estimate of the scales. Let  . So far, we didn't consider the effect of strain. The strain u, ignoring anisotropy, is known to add the term ρuv k 2 . One has to keep in mind that scattering of a FP from a high enough barrier doesn't depend on details, and the cross-section remains 4/k, if k < a −1 . Then, for the linear spectrum the condition for strong localization is  k n 4 i 2 , which is similar to what we have got above. Typically, u is ~10 −4 , and the two terms in ω(k) are of comparable strength for the discussed scales.
So far, we have discussed point-like pinning centers only. In reality, the size of the attached areas can be comparable with the distances between them. Then, owing to the factor f ≈ ka > 1, the energy of the strongly localized FPs can be a few times larger. As we have already mentioned, a reasonable estimate for ω loc str for graphene layer on SiO 2 is ≈0.3 K. Furthermore, effects of the weak localization noticeably expand localization of the FPs. One may easily show . As a result, the energy of localized FPs may increase up to few K.

Effects of anharmonicity. In writing Eq. (1) we have neglected effects of interaction of the FPs with the
in-plane phonons (anharmonicity) and with ripples [32][33][34][35][36] . For graphene on a random substrate, FPs not only scatter from the points of contact with the substrate, but also excite acoustic phonons. One may check, following the calculations of ref. 37 for disorder-assisted scattering, that at low temperatures the effect is vanishingly small.
Next, as it is well known anharmonicity yields a strong effect on the bending rigidity in graphene. The point is that the in-plane rigidity of graphene is extremely high (i.e., it has a very large Young's modulus, Y 0 ), while a bending rigidity κ 0 is relatively modest. The anharmonicity transfers the strong in-plane rigidity into the bending one. The effect is controlled by the temperature, and is of the infrared origin. For small wave vectors  q q th , the bending rigidity driven by thermal fluctuations is scale dependent: κ R (q) ~ κ 0 (q/q th ) −η with the scaling exponent η ≈ 0.8-0.85 38,39 . The transition scale is given by = πκ q th k TY 3 16 B 0 0 2 . For graphene, at room temperature q th ≈ 1.6 nm −1 . As a result, effective bending rigidity of an atomically thin graphene ribbon that are 10-100 micrometers in size at room temperature can be thousands times larger than at T = 0 [40][41][42][43] . (At T = 0, the quantum non-linear effects lead to only logarithmic corrections, hence generally much smaller than the power-law renormalization produced by thermal fluctuations 44 ).
Strong effects caused by the renormalizations, which have been mentioned above, correspond to the vanishing wave vector q → 0. However, above the transition vector q th thermal fluctuations caused by anharmonicity are no longer significant and κ R (q) ≈ κ 0 . In graphene, the energy of a FP with the wave vector q th is equal to Thus, there is a substantial energy gap between the thermal phonons and phonons, for which the effects of the anharmonicity are relevant. For our estimate of ω ≈ . 0 3 loc str K, temperature should be about 10 K or higher to influence our analysis of statistical properties of the FPs on the SiO 2 substrate. Furthermore, in our analysis, we were mostly interested in FPs with the frequency exceeding ω loc str , so that they can propagate between pinned regions colliding randomly with them.

Discussion
We shall discuss now the implications of the FPs on the thermal transport of a graphene layer placed on a corrugated SiO 2 substrate. We argue that traces of localization of the FPs may have been observed in the experiments at low temperatures. Thermal transport. The temperature behavior of overheating in graphene on the top of SiO 2 at low temperatures, refs 37,45-49 , has not been fully understood, yet. Theoretically, the heat flux from electrons to the FPs is known 37 to be ∝T el 3 (correspondingly, the reversed flux from the FPs to electrons is ∝T ph 3 ). We consider the heat exchange with the FPs as realistic explanation of the exponent δ = 3 observed in refs 45,46 at low temperatures and far away from the neutral point. At concentrations of the electric carriers, electrons or holes, n ~ 10 12 cm −2 the alternative explanation of δ = 3 with the use of the result obtained for the case of the unscreened deformation potential 49 is not realistic; a discussion of this point is given below in section "Electron-phonon interaction at low temperatures".
Moreover, we believe that the existence of localized FPs may explain the experimental result of ref. 46 for cooling rate of graphene at the lowest temperatures  .
T 0 85 K. In this regime, the cooling is dominated by the electron heat diffusion along the sample. However, the Lorenz number  estimated in this way, was found to be 35% above its nominal value 0  . It has been shown in refs 50-52 that neither the Fermi liquid nor renormalization-group corrections in disordered 2d electron systems can modify the Lorenz number  and, therefore, the result requests First of all, we recall that interaction of an electron with the FPs is described by the two-phonon processes. Correspondingly, the heat exchange between the electrons and FPs, which is proportional to square of the two-phonon amplitude, contains four powers of the FP-momenta. However, the localized FPs are not goldstone modes anymore. For localized FPs, the momenta that enter into the matrix elements of the electron-FP interaction should be substituted by the inverse of the localization length. As a result, two powers of frequencies in the expression for the heat flux P saturate at ω ω  loc str . This, however, leads to a dramatic consequences. The factor describing the dependence on the occupation numbers in the case of the two-phonon processes for temperatures larger than the energy of the phonons diverges like ω −2 . The diverging integration should be cut-off at energies typical for the localized FPs. As a result one gets a contribution to the cooling rate of the order δ ω The obtained correction to the heat flux has just the form of the electron heat diffusion. In order to obtain experimentally observed magnitude of the deviation of  from 0  , one has to suggest ω ∼ . However, a layer placed on the top of a corrugated substrate, which has been discussed in the present work, is very different from the disordered membranes considered so far. First, disorder here is external rather than internally quenched. Next, disorder pins rigidly the height of randomly chosen points of the layer, rather than acting on the metric and curvature tensors describing the deformation of the membrane. Questions of localization of the FPs to the best of our knowledge have not been addressed previosly.
Graphene layers on top of SiO 2 substrates are expected to play an important role in applications related to thermal transport and for ultrasensitive bolometry. While the graphene sheet is pinned at random points, the FPs may exist in between due to the corrugation typical for this substrate surface. We have argued that in this system Anderson localization of low-energy FPs develop. We showed that the ensemble of flexural phonons in a pinned-suspended flexible sheet is statistically identical to an ensemble of disordered electrons, despite the very different underlying mathematical descriptions. Localization of flexural phonons should be important for thermal transport in such hybrid systems even at not very low temperatures. Traces of localization of the FPs may have been already observed in the experiments. Let us note that we have already shown that FPs give a significant contribution to dephasing rate of electrons in graphene 60 . Here, we have argued that localization of the FPs opens interesting perspectives for thermal transport.

Methods
Scattering of a flexural phonon by a rigid obstacle. According to 9 , the scattering cross-section of the where f(z) is given by the following expression: . The function f(z) with asymptotes f(0) = 1 and is shown in Fig. 5. For short wave-lengths the limiting cross-section is twice larger than the width of the obstacle. Note that the factor f is important for the effectiveness of the weak localization in samples of large size.
Weak multifractality of eigenfunctions. The spatial distribution of wave functions is conveniently characterized by inverse participation ratios: After sample average, 〈P q 〉 shows the scaling behavior with the system size L: Obviously, in the insulating state D q = 0, while in a metal D q = d. At a critical point, D q is a fractional which leads to anomalous scaling behavior in 〈P q 〉. This is a manifestation of the wave function multifractality, which is a consequence of the spatial correlations of the wave function.
In 2d, the inverse participation ratios scale as Here, the deviation of D q from dimension 2 is determined by a small parameter 1/πg. The above result 6 was first obtained by Wegner 27 via the renormalization group calculations for a system of disordered (non-interacting) electrons; see also 28 . The dimensionless parameter g is a conductance of a sample measured in quantum units. In analogous to weak localization, the phenomena is coined "weak multifractality".
We study a model of random pinning centers, in which the graphene sheet is completely attached at the pinning centers; see Fig. 6. In our analysis of multifractality of the FPs we used q = 2. The corresponding inverse participation ratio, P q = 2 , was denoted as IPR. This quantity allowed us to extract the phononic "conductance", which we used for the statistical analysis of our system. Note that, because of the absence of the genuine critical point in 2d, g becomes size-dependent, i.e., g(L) = g 0 − (1/π) ln L/l where l is the mean-free path of FPs at a given energy, and g 0 = g(l). This implies that for each scale L one can use the standard formula given above, but with slowly varying g(L) in the exponent. This is possible because corrections to g are not large in a finite size sample,  and owing to the slow dependence of g on spatial scale L. With this procedure, we have obtained an excellent agreement between the theory of the logarithmic corrections to the conductivity and our numerical results as it is shown in Fig. 2 of the main text.
Energy Level Statistics. Several quantities are introduced to measure the fluctuations of energy levels ω n , such as the distribution function of level spacing P(s), and the level number variance ω Σ Ω ( , ) 2 . Random Matrix Theory (RMT) could be used to describe these quantities in ergodic systems (e.g., for electrons in metallic granules). Here, we will focus on the case of Gaussian Orthogonal Ensemble (GOE). Then the distribution function of level spacing is well described by the Wigner surmise: and Δ is mean level-spacing. In the localized phase the level correlations are absent, and the distribution function of level spacings is Poissonian: P(s) = exp(−s). In our system of the FPs the crossover from the localized to delocalized behavior is illustrated on Fig. 7.
In the RMT, the level number variance ω Σ Ω ( , ) 2 increases logarithmically with Ω. For Ω Δ  , it varies as π , where 〈N〉 = Ω/Δ, and c β is a known constant. This, however, is far from the true behavior (as one can see from Fig. 3 of the main text). For a further analysis it useful that function ω Σ Ω ( , ) 2 is closely related with the two-level correlation function Ω R( ), which is defined as is the Density of States (DOS), and 〈ρ〉 is the average DOS, which is related to the mean level-spacing Δ as ρ = ΔV 1 . The connection between the two correlation functions can be presented as Now, instead of splitting this expression into two parts as it was done in 29 ), we use Eq. 8 to calculate the level number variance. After integration in Q and s, this yields . We use the expression determined by formula 11 for fitting the numerical data presented in Fig. 3 in the main text.

Statistical Properties of the Wave Functions.
Porter and Thomas were first who studied the distribution of eigenfunction amplitudes within the RMT framework. Their result shows simply Gaussian distribution, which leads to the following distribution of the intensities y: RMT y/2 is normalized in such a way that 〈y〉 = 1. The supersymmetric field theory was applied to the study of the eigenfunction statistics in a d-dimensional disordered system. The eigenfunction intensity ψ = | | y V r ( ) 2 0 in a point r 0 is distributed as: For d = 2, and not too large y, one can calculate perturbatively the deviations from the RMT distribution P(y). The corrected distribution function was found by Fyodorov and Mirlin in (31) . Notice that the expression in the square brackets is non-monotonous. A peculiar behavior of this expression as a function y can be seen in Fig. 4 of the main text.

Electron-Phonon Interaction at Low Temperatures.
Some recent experiments interpreted the cooling rate in terms of the weakly screened electron-phonon (e-ph) interaction down to very low temperatures. In particular, in 46 the T 3 law has been observed down to 0.5 K for sample D3. Such a behaivior corresponds to the dirty regime with negligible screening of the e-ph deformation potential according to the Table 1 presented in 49 . We believe, however, that applicability of the screenless approximation for the e-ph interaction at such low temperatures and concentration of carriers n ~ 10 12 cm −2 is highly questionable.
To illustrate our point, let us estimate, following 49 , temperature above which screening becomes irrelevant in the dirty regime  < = T T slk / dis B . Comparing the expressions for the energy flux at strong and weak screening 49 , one finds that crossover from the T 5 to the observed T 3 behavior should happen at > ⁎ T T which satisfies 2 for the Bloch-Gruneisen temperature and κ for an effective dielectric constant (for graphene on SiO 2 substrate, κ ≈ 3). Here we expressed the result in terms of T BG using the fact that in graphene the inverse of the screening radius is of the order of the Fermi momentum. The unscreened T 3 behavior may occur at temperatures higher than ⁎ T assuming that it is less than T dis . (For samples D1 and D3 in (46) one has T dis ≈ 40 K).
For the aforementioned sample D3, T BG ≈ 80 K. Thus, one may expect for the crossover temperature ≈ ≈ κ ⁎ T 80 250 K. It is clear that pushing ⁎ T down to 1 K region for the discussed densities of charge carriers would require an unrealistic value of κ. We, therefore, believe that this mechanism may be discarded for explanation of the T 3 scaling of the electron-phonon heat flux observed at low temperatures.

Data Availability
All data needed to evaluate the conclusions in the paper are present in the paper. Additional data related to this paper may be requested from the authors.