Dynamic dielectric function and phonon self-energy from electrons strongly correlated with acoustic phonons in 2D Dirac crystals

The unique structure of two-dimensional (2D) Dirac crystals, with electronic bands linear in the proximity of the Brillouin-zone boundary and the Fermi energy, creates anomalous situations where small Fermi-energy perturbations critically affect the electron-related lattice properties of the system. The Fermi-surface nesting (FSN) conditions determining such effects via electron–phonon interaction require accurate estimates of the crystal’s response function (χ) as a function of the phonon wavevector q for any values of temperature, as well as realistic hypotheses on the nature of the phonons involved. Numerous analytical estimates of χ(q) for 2D Dirac crystals beyond the Thomas–Fermi approximation have been so far carried out only in terms of dielectric response function χ(q,ω) , for photon and optical-phonon perturbations, due to relative ease of incorporating a q-independent oscillation frequency (ω) in calculation. Models accounting for Dirac-electron interaction with acoustic phonons, for which ω is linear to q and is therefore dispersive, are essential to understand many critical crystal properties, including electrical and thermal transport. The lack of such models has often led to the assumption that the dielectric response function χ(q) in these systems can be understood from free-electron behavior. Here, we show that, different from free-electron systems, χ(q) calculated for acoustic phonons in 2D Dirac crystals using the Lindhard model, exhibits a cuspidal point at the FSN condition. Strong variability of ∂χ∂q persists also at finite temperatures, while χ(q) tend to infinity in the dynamic case where the speed of sound is small, albeit non negligible, over the Dirac-electron Fermi velocity. The implications of our findings for electron-acoustic phonon interaction and transport properties such as the phonon line width derived from the phonon self-energy will also be discussed.


Introduction
Dirac crystals are a broad class of zero band-gap solids in which the electronic band structure is linear in the crystal momentum, k, instead of quadratic, as commonly observed in metals and semiconductors. [1,2,3] This leads to dispersionless electrons with a behavior reminiscent of photons on the Dirac light cone in relativity. [4] In twodimensional (2D) Dirac crystals, the best known of which is graphene, the density of electronic states near the Dirac point also linearly tends to zero, thus enabling additional remarkable properties, including extreme carrier mobility and quantum Hall effects, and the possibility of topologically insulating characteristics. [5,6,7,8] The experimental discovery of a plethora of new 2D Dirac crystals in recent years [9,10,11,12,3,13] compels their deeper understanding.
In undoped Dirac crystals, the valence and conduction bands meet at the Brillouin-zone boundary Dirac K-point. The Fermi energy E F coincides with the energy of this point, with the Fermi surface area collapsing to zero. Thus, the crystal's Fermi surface degenerates into a single point of the electronic band structure. [1,2] Such a degeneracy enables us to tune the density of electronic states at E F , as well as the system's electrical and thermal conductivity, via external electric fields or tunable doping-an effect that leads the carrier density to increase by orders of magnitude upon small fluctuations of E F . [14,15] These Fermi-level shifts are also expected to dramatically alter the strength of the interaction between charge carriers and lattice phonons, [16,17] which may have profound effects on the applicability of the Born-Oppenheimer approximation, leading to highly specific dielectric properties from electrons strongly correlated with acoustic phonons in 2D Dirac crystals. [18,19,20] A key parameter to understand the effects of strong electron-phonon correlation is the crystal's dielectric response function, χ(q). χ(q) differs from the dielectric susceptibility χ(q, ω) in that it takes into account the dispersion relationship ω q of the phonons for which it is calculated and cannot simply be obtained from replacing ω q into χ(q, ω), with which processes violating the conservation of energy or momentum would be inappropriately considered. For phonons of wavevector, q, χ(q) originates from the superposition of all of the inelastic scattering of electrons and holes by such phonons. The customary approach for calculating χ(q) with appropriate selection rules relies on the Thomas-Fermi (also known as) Debye-Huckel [21,22,23] approximation, which assumes long wavelength and a short wave number from the scattered phonons. Such an assumption is specifically designed for metals with large (> 10 22 cm −3 ) electron densities and large Fermi surfaces in the proximity of E F and, therefore, q much shorter than the dimensions of the Fermi surface, implying short screening length. [24] On the contrary the Fermi surface area in undoped 2D Dirac crystals collapses to zero, as already stated. The screening length may tend to infinity and the Fermi wavenumber k F corresponding to the Fermi surface radius, is expected to remain very small also in the presence of external electric fields and doping, always leading to free-electron densities 10 19 cm −3 orders of magnitude below the expected applicability of the models [23] suitable for metals. This discounts the applicability to 2D Dirac crystals of early models based on the polarizability of 2D gases in the context of short-range screening. [24] Because 2D Dirac crystals are zero-band gap semiconductors they dramatically amplify the effects of any phonon disturbance for which q ≈ 2k F . This effect, known as Fermi-surface nesting (FSN) [25,26] and potentially leading to Kohn anomalies [27,28] has been often considered using the Lindhard model. The Lindhard model is a method for calculating the effects of electric field screening by electrons in solids based on first-order quantum perturbation theory and it accurately predicts a common limitation of most of these calculations. For example, Kohn anomalies and FSN conditions for free electron gases of any dimensionality are correctly predicted using a static Lindhard model at 0 K [29]. In 2D Dirac crystals, the Lindhard model of electron screening by phonons has been considered by many authors such as [30,31,32,33,34,35,36,27], however, because of the use of a q independent from the disturbance energy ω 0 in these reports none of them are well suited to estimate the Lindhard response function for acoustic phonon branches in 2D Dirac crystals. This is because a common limitation of most Lindhard model calculations is the use of the random phase approximation (RPA) in which the contribution to the dielectric response function from the total electric potential is assumed to average out so that only the potential at wave vector k contributes. This considers only relatively weak electron screening potential and does not take into account the dispersion relationship ω q of the phonon mode for which it is calculated, and cannot simply be obtained from replacing ω q into χ(q, ω).
Acoustic phonons, for which the energy ω q = c · q is proportional to the wavenumber via c, which is the speed of sound written in energy units are present in any crystalline lattice regardless of the number of atoms within their basis. Because acoustic phonons are responsible for the long-wavelength, low-energy end of the vibrational density-of-states, they are critically important for a host of measurable properties including but not limited to thermal transport [37,38], electric conductivity [39,40], and Brillouin scattering [41,42], particularly at low or room temperatures where the optical phonons are not excited. Such properties are often calculated in 2D Dirac crystals by considering weak electron-phonon interaction and, consequently, short-range charge screening effects at the level of Thomas-Fermi. [43,44] This is unreliable in the very frequent case in which small Fermi-energy fluctuations by doping or impurities produce Fermi level shifts bringing k F within the range of acoustic phonon wavenumbers q. This may lead to FSN and anomalous electron-phonon interaction. For example, Ramezanali et al [45] found that the accuracy of specific-heat calculations in graphene can be remarkably improved by introducing a Lindhard-based correction still assuming random-phase approximation and relatively weak electron-phonon coupling, which limits the generality of their approach. Furthermore, by deriving the dielectric response function as a function of q using the Lindhard model we can improve on the phonon self-energy calculation performed in 2D Dirac crystals. The phonon line width derived from the phonon self-energy provides a way to gain experimental information about the electron-phonon coupling strength. The calculations for the phonon self-energy have been mostly carried out at T → 0 [46] and q = 0 [47] which limits the application of such results.
The objective of our work is to analytically calculate the dielectric response function of electrons strongly correlated with acoustic phonons by using the Lindhard model beyond the RPA in 2D Dirac crystals. Using scaling laws and introducing reduced phonon and electron wavevectors, respectively ψ and ξ, as well as a reduced temperature τ and (where v F = E F /k F is the Fermi velocity) we show that a universal scaling law for the Lindhard dielectric response of acoustic phonons only depending on the dimensionless quantities in Eq.(1) can be established. Such an expression will be general enough to describe the electron-lattice interaction at any Fermi-level shifts and temperatures even for cases where the small-τ Sommerfeld approximation [22] customarily used in solid state physics is not valid. Further, using the derived dielectric response function we calculate the phonon line width from the phonon self-energy and compare our results with experimental data presented in the literature.

Methodology
Using the normalized quantities in Eq.(1) the Lindhard dielectric response function [22,24] takes the form: where f ξ±ψ/2 are the occupation function of the single particle state following the Fermi-Dirac distribution and ω ψ is the phonon energy. The integral is over the momentum space of the Fermi sphere and E ξ±ψ/2 are the energy of the created and annihilated electrons, respectively. We do not integrate over the electron spins and have simplified the problem by multiplying Eq.(2) by a factor of 2. As shown in Fig. 1(a) the π-π * electron energy spacing in a 2D Dirac crystal at the zone center E Γ [48,49] is too high relative to the energy of acoustical phonons [50] E ph = ω ψ for any electron-phonon interactions to take place without violating the conservation of energy. The electron energy spacing becomes comparable to the energy of phonons at the K-point thus making the electron-phonon interactions possible. Assuming a honeycomb lattice structure as shown in Fig. 1(b) the electron-phonon interactions in the reciprocal lattice shown by the shaded grey region occur at the two-zone boundary points K1 and K2 depicted by the Dirac cones. Figure 1: (a) The π-π * electron energy spacing along a straight line from the center to the bounds of a 2D Dirac lattice. The electron energy spacing at the zone center is much larger than the phonon energy making the electron-phonon interactions impossible. The electron energy spacing decreases and becomes comparable to the phonon energy as we move toward the lattice zone boundaries making the electron-phonon interactions possible. (b) The reciprocal lattice of a honeycomb structured lattice is shown as the grey-shaded region where we have two Dirac cones in which electron-phonon interactions take place at.
To derive an analytical expression for the Lindhard dielectric response function we must solve the difference between the occupation number of the levels above and below the Fermi energy, which is provided by the Fermi Dirac statistics in the numerator of the integral in Eq.(2) by applying a suitable approximation. The most used approximation for simplifying the Fermi-Dirac distribution function is the "single step" function used at T ≈ 0 K where the Fermi-Dirac distribution has a value of 1 for energies below the Fermi energy, and a value of 0 for energies above. [51] However, at finite temperatures the "single step" function approximation is no longer accurate since the distribution gets smeared out, as some electrons begin to be thermally excited to energy levels above E F . We, therefore, approximate the Fermi-Dirac distribution with a "double step" function which has three values 1, 1/2, and 0 and covers the distribution of electrons at different energies more accurately, shown in Fig. 2(a). Furthermore, while at T = 0 K the chemical potential µ is equal to E F this is no longer the case at finite temperatures where µ becomes dependent on both E F and T. Therefore, to find the values of the Fermi-Dirac distribution function at different energies using the "double step" function we must find a relation between µ, E F and T.
In metals where we have large free electron density in the proximity of the Fermi energy E F , we use the Sommerfeld expansion to write the chemical potential, µ, as a function of E F and T [22]. However, Sommerfeld expansion cannot be reliable for 2D Dirac crystals which have a limited concentration of free electrons and we need to express small fluctuations of µ over E F . Therefore, to write µ in terms of E F and T for 2D Dirac crystals we have to develop a new relation. To this end, we define ∆E ± to be the range of the singly occupied energy levels in the Dirac cone situated above and below the Fermi energy, respectively, shown in Fig. 2(b). The integrated density of states over ∆E ± are equal and are defined as ∆g ±(T) . By analyzing Fig. 2(b) we can write µ as a function of E F and ∆E ± in the following manner: To write ∆E ± in terms of E F and T we first write the sum of the singly occupied energy levels in the Dirac cone situated above and below the Fermi energy level, Fig. 2(b), as follow: We further know that the integrated density of states over ∆E + and ∆E − are equal. This results in the surface area of the two red and blue regions on the Dirac cone in Fig. 2(b) to be equal to one another and result in the following equation: By inserting Eq.(4) into Eq.(5) we find a relation between ∆E ± and E F and T. We can therefore write the chemical potential as: To find the difference between f ξ+ψ/2 and f ξ−ψ/2 versus the energy we use a "double step" function approximation as shown in Fig. 2(a). By computing the area under the curve for the two occupation levels shown in Fig. 2(c) we have: Where the energy bounds in Eq. (7) for which f ξ+ψ/2 − f ξ−ψ/2 = 1/2 is shown by the purple shaded region in Fig. 2(c). The difference between the two, sets the bounds of integration for calculating the dielectric response function at non-zero temperature. By analyzing Eq.(7) at T ≈ 0 K we observe that the two purple regions in Fig. 2(c) overlap with one another, and the "double step" function turns into a "single step" function further confirming our results. The difference between f ξ+ψ/2 and f ξ−ψ/2 as a function of energy in a 2D Dirac crystal shown in the two purple squares which set the bounds of integration for calculating the dielectric response function.
By raising the temperature, the Fermi energy of the 2D Dirac crystal shifts leading the bounds of the integration in Eq.(2) for the dielectric response function to also change. While the bounds of the integration for ξ(ψ, ω) at T = 0 K are between 1/2 ± ψ/2, using the "double step" function to derive the difference in the Fermi-Dirac distribution for T = 0 K we have: Where we have written the dielectric response function in dimensionless coordinates as χ(ψ, ω) = χ(q, ω)/2k F . We see that by approximating the occupation function of the electron with a "double step" function we split the integral in Eq.(2) into two integrals where the bounds of each integral is one of the purple shaded regions depicted in Fig. 2

(c).
A disadvantage of the Lindhard model over more simplistic approximations such as Thomas-Fermi, is the requirement of carrying out integrations that often cannot be performed analytically, and may become cumbersome where: i) the model is dynamic-i.e. one assumes that phonons not only carry momentum, but also energy ω 0 ; or ii) nonzero temperature is considered-a necessary requirement to compare χ(q) with experiments involving, for example, T-dependent electrical or transport measurements; or iii) the phonon energy ω q is wavevector-dependent in a dispersive system-a respect in which it is worthwhile noting that the phonon energy is not only required in dynamic Lindhard-model calculations but also static ones, as it affects the electron Fermi-Dirac distribution, also at 0 K. Of course, the complexity of the calculation will increase if more than one of the phenomena (i − iii) are considered. To simplify the integral we assume the extreme carrier mobility in 2D Dirac crystals results in their Fermi velocity to be much larger than the speed of sound, v F >> c. This has been experimentally confirmed in various 2D Dirac crystals such as graphene [52], borophene [5], Weyl semimetals [53], and silicene [54] where c < 0.01v F . Defining w F = c/v F we can expand the dielectric response function around w F with w F → 0 in the following manner: The first addend in Eq.(9) is the static dielectric response of the 2D Dirac crystal with the phonon energy being equal to zero. The next addends are the higher order terms which are the dielectric response at nonzero phonon energy. Due to computational convenience in the calculation, we shall derive the static term in the Cartesian coordinates and the dynamic terms in polar counterparts. In the end, by adding the static term with the dynamic terms we get the dynamic dielectric response function of a 2D Dirac crystal. Figure 3: (a),(b) The shift in the momentum sphere respectively for q < 2k F or ψ < 1 and q > 2k F or ψ > 1 written in the Cartesian coordinate system. We write the shift in the momentum sphere in Cartesian coordinates to calculate the static dielectric response function. (c), (d) The shift in the momentum sphere respectively for q < 2k F or ψ < 1 and q > 2k F or ψ > 1 written in the polar coordinate system. We write the shift in the momentum sphere in polar coordinates to calculate the 2nd to last addends of the dielectric response function which contain the phonon energy.

Results
To derive the dielectric response function, we assume we have a shift in the momentum sphere [55,56] equal to q = (2k F )ψ, caused by a disturbance in the lattice. Fig. 3 reports the shift in the momentum sphere written in the Cartesian and polar coordinate systems. Two different trends for the dielectric response function are obtained in the two regions at q < 2k F or ψ < 1 and q > 2k F or ψ > 1, respectively.

Static dielectric response
The integral in Eq.(10) expresses the static dielectric response, χ 0 (ψ), function of a 2D Dirac beyond the RPA in the two regions of ψ < 1 and ψ > 1. Following Eq.(2) χ 0 (ψ) of a 2D Dirac crystal at T = 0 K is equal to: Since we integrate over the two momentum spheres Fig. 3, we multiply χ 0 (ψ) in Eq.(2) by a factor of 2. Furthermore, since the bounds along the y-axis change when writing the integral in the Cartesian coordinate system for ψ < 1, we divide the integral into two, with each integral covering a separate region as shown in Fig. 3(a) by the two blue and green colors. So, the first two integrals in Eq.(10) are for the region ψ < 1 and the third integral is for the region ψ > 1.
Analyzing each of these two regions separately, for ψ < 1 we have: where we have integrated over ξ y . The employment of numerical methods in different branches of physics and science has become more common in recent years [57,58,59,60,61]. Using a MATLAB routine, we solve Eq.(11) numerically by sums. The designated solid black line in Fig. 4(a) numbered (I) is the integral solution of Eq. (11). To further analyze the integral and come up with an analytical solution we partition the different parts of Eq. (11). The red dashed lines numbered (i) and (ii) are the two positive and two negative terms containing the arc hyperbolic sinus. The sum of these two terms is (i) + (ii) = (I 3 ) shown in Fig. 4(a) as a solid red line.
It can be seen in Fig. 4(a) that the designated line (I 3 ) is approximately a constant equal to (I 3 ) ≈ 0.39. This can be proven by Taylor expanding (I 3 ) up to the first order about ψ = 0. We further separate the remaining terms in Eq. (11) in the following manner: Where (I 1 ) is shown in a purple solid line and (I 2 ) is shown in a blue solid line in Fig. 4(a). By analyzing Fig. 4(a) we observe that (I 2 ) can be estimated as (I 2 ) ≈ 0.2 + (I 1 ) / (2ψ). This can also be proven by Taylor expanding (I 1 ) and (I 2 ) up to the second order about the midpoint of the integral bounds which are (1 + ψ)/4 and (1 -ψ)/4. We can therefore write Eq. 11 as: Where Eq. (15) is the analytical expression of the static dielectric response of a 2D Driac crystal in the region of ψ < 1.
For ψ > 1 we have: The designated solid black line numbered (I ) in Fig. 4(a) is the integral solution of Eq.(16) solved numerically. We again partition the different parts of the integral to come up with an analytical solution. The red dashed lines in Fig. 4(a) numbered (i ) and (ii ) are the first and second terms containing the arc hyperbolic sin term in Eq. (16). The sum of these two terms is (i ) + (ii ) = (I 3 ), shown in Fig. 4(a) as a solid red line.
This can be proven by Taylor expanding (I 3 ) up to the first order about the midpoint of the integral bound which is ψ/2. We further separate the remaining terms in Eq. 16 as following: Where (I 1 ) is shown by a purple solid line while (I 2 ) is shown by a blue solid line in Fig. 4(a). By analyzing Fig. 4(a) we observe that (I 2 ) can be estimated as (I 2 ) ≈ [0.4 + (I 1 )] / (2ψ). This can also be proven by Taylor expanding (I 1 ) and (I 2 ) up to the first order about their midpoint integral bound which is ψ/2. We can therefore write Eq.(16) as: Where Eq. (20) is the analytical expression of the static dielectric function of a 2D Driac crystal in the region of ψ > 1.
The comparison between the analytical and numerical solution of the static dielectric response using the Lindhard model is shown in Fig. 4(b) in which the two solutions are in good agreement with one another. By studying both graphs we see their cuspidal points are at ψ = 1 or q = 2k F which is the Fermi wave number. We, therefore, have the absolute value of the static dielectric response function increase up to ψ = 1 and then decrease as we go further away.

Dielectric response at nonzero phonon energy
The integral in Eq.(21) expresses the second addend of the dielectric response function in Eq. (9). To derive ∂ χ0(ψ) ∂w F of a 2D Dirac crystal we perform the calculations in a polar coordinate system. Therefore, following Eq.(2), ∂ χ0(ψ) ∂w F of a 2D Dirac crystal at T = 0 K is equal to: The first integral is for the region ψ < 1 and the second integral is for the region ψ > 1. ψ M is the maximum distance we place ourselves from the Fermi sphere to do the calculations for ψ > 1. Performing the calculations we observe that the second addend of the dielectric response function goes to zero for ψ > 3 so ψ M can be confidently put equal to 3.
To solve the integral in Eq.(21) analytically we use the binomial expansion: Applying the binomial expansion up to the first order of magnitude to the first addend of Eq. (22) designating the solution for ψ < 1 we get: Where Eq.(A.4) is the analytical expression of the first integral in Eq. (21). We have derived the steps to calculate Eq.(A.4) in the appendix for completeness. Applying the binomial expansion up to the first order of magnitude to the second addend of Eq.(21) designating the solution for ψ > 1 we get: Where Eq.(24) is the analytical expression of the second integral in Eq. (21). We solve Eq.  Fig. 4(c). As is shown in Fig. 4(c) the two results are in good agreement with one another. By analyzing Fig. 4(c) we see that the second-order addend of Eq.(9) of the dielectric response function diverges at ψ ≈ 1 making the study of this point important. Higher order terms of the dielectric response function present in Eq.(9) can be derived using the same analytical technique. These higher order terms will give us similar plots with divergence at ψ ≈ 1. However, due to the presence of the term w F n in the numerator where w F << 1 these higher-order terms can be neglected. By knowing the first and second addend of Eq.(9) which are respectively the static dielectric response function and the dielectric response function at nonzero phonon energy we can write the total dielectric response function in the two regions of ψ < 1 and ψ > 1. For the region ψ < 1 we insert the values derived for χ 0 (ψ) Eq. (15), and ∂ χ0(ψ) ∂w F Eq.(A.4), into Eq. (9). Therefore, the analytical solution of the total dielectric response function of a 2D Dirac crystal for ψ < 1 is equal to: And for the region ψ > 1 we insert the values derived for χ 0 (ψ) Eq. (20), and ∂ χ0(ψ) ∂w F Eq.(24), into Eq.(9). Therefore, the analytical solution of the total dielectric response function of a 2D Dirac crystal for ψ > 1 is: We have therefore derived the analytical expression for the dielectric response function of a 2D Dirac crystal in the two regions of ψ < 1 and ψ > 1 at T = 0 K. We can derive the dielectric response function of a 2D Dirac crystal for T = 0 K by employing Eq.(8) and following the same procedure.
For deriving the dielectric response function at T = 0 K we approximate the Fermi-Dirac distribution by a "single step" function. However, to derive the dielectric response function at higher temperatures we employ Eq.(8) where we approximate the Fermi-Dirac distribution with a "double step" function. This results in the dielectric response function having two integrals with their bounds shown by the two purple-shaded regions depicted in Fig. 2(c). Each of these two integrals can be solved analytically by following the same procedure we used to solve the integral of the dielectric response function at T = 0 K. To study the dielectric response function for T = 0 K we also define the variable τ = (k B T/2E F ) presented in Eq.(1). The static dielectric response function, χ 0 (ψ), of a 2D Dirac crystal derived analytically using the "double step" function for T = 0 K is shown in Fig. 5(a-d) by the solid red lines. We observe that although χ 0 (ψ) at T = 0 K has only one cuspidal point, for T = 0 K the number of cuspidal points increases to two with each integral giving us a separate cuspid. By analyzing Fig. 5(a-d) we observe that as we increase the temperature for a given Fermi energy the distance between the cuspidal points increases. This is because as we increase τ the regions covered by the two integrals in Eq.(8) get further apart from each other which results in the cuspidal points of the two integrals to also get further apart from one another. To study the nature of the two cuspidal points we derive χ 0 (ψ) numerically without applying any approximations to the Fermi-Dirac distribution function and compare it with the analytical solution. The numerical solution is shown in Fig. 5(a-d) by the black dashed lines. By analyzing the numerical solution we observe that we only have one cuspidal point for τ << 1 and the cuspidal point starts to vanish as we increase τ . This is not the case for the analytical solution were although as we increase τ the value of the cuspidal points decrease but they do not vanish. We, therefore, conclude that the two cuspidal points derived analytically are an artifact of our solution which is caused by using a "double step" function to approximate the Fermi-Dirac distribution. By analyzing Fig. 5(a-d) we observe that the analytical expression using the "double step" function gives us accurate estimations for τ << 1, and as we increase τ the analytical solution becomes less accurate, however, it still remains within a good approximation range of the numerical value. This is because although the two cuspids in the analytical solution don't vanish but their values decrease tending towards the numerical solution.
By deriving the analytical dielectric response function of 2D Dirac crystals as a function of the phonon wave number using the Lindhard model we are able to calculate the phonon line width of these class of condensed matter systems more accurately. To this end, we first write the first order phonon self-energy as a function of ψ as follow: [46] where M ψ is the electron-phonon matrix element [21], and N ψ is the number of points in the summation over ψ. The finite line width or the inverse lifetime of a phonon mode is connected to the imaginary part of the phonon self-energy and can be written as: To calculate the phonon line width we write the electron-phonon matrix element M ψ , as a function of the screening potential in the following manner: [21] where Ω is the unit cell area of the 2D Dirac crystal, m is the ion mass, and φ s (ψ) is the screening potential. The screening potential is a function of the dielectric response function and is equal to: [22] φ s (ψ) = 1 2k F The analytical solution is done by assuming the Fermi-Dirac distribution function to be a "double step" function while the numerical solution makes no such approximations. We observe that although the analytical solution always has two cuspidal points the numerical solution only has one cuspidal point for τ << 1 and as we increase τ the single cuspidal point in the numerical solution vanishes. Therefore, the two cuspidal points derived analytically are an artifact of our solution which is caused by using a "double step" function to approximate the Fermi-Dirac distribution.
where Q is the ion charge. By inserting Eq.(29) and Eq.(30) into Eq.(28) we get the phonon line width as a function of ψ in the following manner: We observe that the role of the dielectric response function in the denominator of Eq.(31) is essential in determining the phonon line width accurately. By knowing χ ψ , we can write the phonon line width as a function of ψ for acoustic phonons with dispersive energy wavelength relationships. At zero temperature we solve the difference between the Fermi-Dirac distributions using a single-step function and inserting the dielectric response function derived from Eq. (25) and Eq. (26). For nonzero temperatures, we use the double-step function and insert the value of the dielectric response function derived from Eq.(8) for different temperatures. In the next section, we analyze the dielectric response function and the phonon line width of 2D Dirac crystals extensively.

Discussion
We have analytically derived the dynamic dielectric response function of electrons strongly correlated to acoustic phonons in 2D Dirac crystals as a function of the phonon wave number, Eqs. (25,26). Our expression of the dielectric response function is general enough to describe the electron-lattice interaction at any Fermi-level shifts and temperatures even for cases where the Sommerfeld approximation is not valid. Also, as can be seen from Eqs. (25,26), our expression of the dielectric response function has the advantage of being applicable to various 2D Dirac crystals with different Fermi velocities, v F , and not only to a specific type of 2D Dirac crystal. The analytical expression of the dielectric response function derived in our paper for acoustic phonons where ω q = cq is more accurate compared to methods that use the RPA method. For example, using the RPA Wusch et al. [30] states that at q = 2k F or ψ = 1 the dielectric response function has a discontinuity only in the second derivative while we show that in the static regime, the dielectric response function has a discontinuity in the first derivative and in the dynamic regime the dielectric response function diverges itself at ψ = 1, Fig. 6(a). A literature review of the work done on the dielectric response function in 2D Dirac crystals is shown in Table 1. Furthermore comparing our results to the Lindhard-Mermin dielectric function [62] we observe that using RPA Mermin derives the Lindhard dielectric function in the relaxation time approximation by assuming the relaxation to occur not towards the thermal equilibrium but to a more general distribution characterized by the chemical potential µ. The advantage of the Mermin approach is its simplicity, however, it demands further systematic elaboration of the quantum statistical theory of the dielectric response function and a more comprehensive formulation of the chemical potential beyond the Sommerfeld approximation for 2D Dirac crystals. Also, our analytical expression of the dielectric response function written as a function of the phonon wave number is necessary to better understand some of the features of 2D Dirac crystals such as the phonon self-energy. The customarily used model to derive the phonon self-energy is the jellium model [21]. Jellium models disregard the lattice structure and only consider free electrons. Furthermore, the ions in the jellium are approximated into a uniform background of positive charges resulting in a constant electrostatic potential. Such a model is suitable for studying electron-electron and electron-phonon interactions in metals but not in 2D Dirac crystals where the electrostatic potential undergoes strong fluctuations. Specifically, in 2D Dirac crystals the electronic dispersion is linear (i.e. E = v F k as opposed to E = 2 k 2 /2m as in the jellium). In order to consider this, the electrostatic potential is assumed to be Lindhard-screened, as in Eq. (30). This approach, which also considers the specific crystalline structure of 2D Dirac crystals leading to Fermi-level electrons at the Brillouin-zone K-point, has been used by us to derive the phonon line width in our model as in Eq. (31). It is worthwhile noting that the non-uniformity of the electrostatic potential plays a critical role in determining the phonon line width at wave numbers q approaching Fermi surface nesting conditions, where we anticipate it to increase linearly with the temperature in agreement with the experimental results [63]. This would not be the case if the electron-phonon coupling would have been considered at the approximation level of the jellium model.

Literature Reported Results Limitations
Hwang [31] Bahrami [32] Have calculated the dielectric response function by assuming phonons to have non-dispersive energy-wavelengths Have neglected phonons with dispersive energy wavelengths.
Iurov [34] Have calculated the dielectric response function for silicene illuminated by circularly polarized light.
The disturbances studied are, photons with fixed energy and not phonons and the study on silicene cannot be generalized to other Dirac crystals.
Lu [35] Have calculated the zero-temperature, static, Lindhard response function for single-layer and bilayer graphene.
Only assumes the static dielectric response function at zero temperatures. Also, the study on graphene cannot be generalized to other Dirac crystals.
Calandra [36] Have calculated the e-ph coupling in electron-doped graphene using e-ph matrix elements extracted from density functional theory simulations.
A direct analytical calculation has not been made and the results can not be extended beyond the specific crystal under study.
Wunsch [30] Have used RPA to calculate the response the function of graphene in the two scenarios of q → 0 relevant for photon spectroscopy and ω q = 0 applicable for the static case.
A direct analytical calculation has not been made, phonon wave functions where q = 0 have not been considered and the results cannot be extended beyond graphene.
Lazzeri [27] Have derived the dielectric response function in doped graphene as a function of the charge doping for q = 0 Phonon wave functions where q = 0 have not been considered, and the results cannot be extended beyond graphene. We will now analyze the dynamic dielectric response function along the ψ axis and discuss the physical meaning of the observed cuspid and divergence near the FSN region. As shown in Eq.(9) the dielectric response function of a 2D Dirac crystal is the sum of the static term with zero phonon energy, w F = 0, and higher order terms where the phonon energy is nonzero, w F = 0. In Fig. 6(a) we plot the dielectric response function of a 2D Dirac crystal as a function of ψ at T = 0 K for different values of w F . By analyzing Fig. 6(a) and comparing the red dashed line where w F = 0, with the blue and orange dashed lines where w F = 0, we observe that at the two regions ψ < 1 and ψ > 1, the static term is dominant, however, near the FSN region, ψ ≈ 1, the higher order terms become more dominant. We observe that the dielectric response function diverges near the FSN region at the point ψ ≈ 1 even for small values of w F = 0 where the speed of sound is small, albeit non-negligible, over the Dirac-electron Fermi velocity. It is only when w F = 0 and we assume that we are in a static regime that the strong variation in the dielectric response function vanishes, and we only have a cuspid. The divergence of the dielectric response function in the dynamic regime and the cuspid in the static regime near the FSN region at τ ≈ 0 is a consequence of the particular enhanced electron-phonon interaction at q ≈ 2k F or ψ ≈ 1. However, in the static regime, we observe that the cuspid vanishes as we increase τ as shown in Fig. 5. This is because, as we increase the temperature more electrons begin to get thermally excited to energy levels close to E F smearing out the "single step" Fermi-Dirac approximation making the particular enhanced electron-phonon interaction more general which smooths out the cuspid observed in Fig. 5(a). It should be noted that due to the strong variations in the dynamic dielectric response function near the FSN region of 2D Dirac crystals as shown in Fig. 6(a), it is important not to simplify the problem to the static case where the energy of acoustic phonons is put equal to zero. Such physical settings can be the study of the thermal and electrical variations of 2D Dirac crystals caused by electron-phonon interactions [20]. Also, when studying the higher order terms of electron-phonon interactions in 2D Dirac crystals with the employment of a Fermion propagator, [21,64] which has a time component, we have to consider the energy of the phonons and cannot work in the static regime. [20] By further analyzing Fig. 6(a) we observe that the difference in the dielectric response function of various 2D Dirac crystals with different values of w F [5,52,53,54,65], is negligible showing the generality of our solution.
We will further discuss the effect of the disorder on the dielectric response function near the FSN region for different values of k F . Disorder affects 2D Dirac systems differently at different values of k F and q introducing δk F and δq both of the order of 1/L dis (with L dis being the order range of the disorder). Studying our expression of the dielectric response function written as a function of ψ, disorder significantly affects our considerations when ∆ψ < δψ. ∆ψ is the full-width height-maximum (FWHM) of the dielectric response function near the FSN region and as can be seen from Fig. 6(a) in the static case (red dashed line) ∆ψ stat ≈ 0.2 and in the dynamic case (blue dashed line) ∆ψ dyn ≈ 0.1. We derive δψ by taking the partial derivative of ψ as follow: We assume δk F = δq ≈ 1/L dis and ψ ≈ 1 near the FSN region. We therefore have: To study the different values of k F where the disorder becomes important we assume the disorder range to be L dis = 10 −5 m. For k F = 10 8 m −1 we have δψ = 0.001. Therefore, δψ is smaller than both ∆ψ stat and ∆ψ dyn resulting in the effect of the disorder to be negligible. Making the value of k F smaller for k F = 10 6 m −1 we have δψ = 0.15. In this scenario for the static dielectric response function, we have δψ < ∆ψ stat = 0.2, and although the disorder is significant our results are still reliable. However, for the dynamic dielectric response function, we have 0.1 = ∆ψ dyn < δψ and our results become unreliable. From Eqs.(33 we observe that as k F → 0, the role of the disorder becomes more significant and in order to counterbalance it, the length of the disorder range should approach infinity, L dis → ∞. We further compare the static dielectric response function, χ 0 (ψ), of a 2D Dirac crystal with that of a Fermion gas at T = 0 K [29]. This is shown in Fig. 6(b). We observe that χ 0 (ψ) of a 2D Dirac crystal is reminiscent of that of a 1D and a 2D Fermion gas. As can be seen from Fig. 6(b) χ 0 (ψ) of a 2D Dirac crystal lies between that of a 1D and 2D Fermion gas not diverging near the FSN region like the former and not approaching it along a constant line like the latter. Hence, when studying χ 0 (ψ) of a 2D Dirac crystal it is important to distinguish it from a 2D Fermi gas. While χ 0 (ψ) of a 2D Fermi gas approaches the FSN along a constant line, for 2D Dirac crystals χ 0 (ψ) increases reaching a maximum at ψ = 1. This makes the study of χ 0 (ψ) for 2D Dirac crystals more intricate in the FSN region. By further analyzing Fig. 6(b) we observe that in the region of ψ >> 1, χ 0 (ψ) of the Fermi gas and 2D Dirac crystals are both proportional to 1/ψ which is equivalent to the result we get when we use the Thomas-Fermi approach. Therefore, the Thomas-Fermi model for calculating the dielectric response function at long wavelengths is equal for 2D Dirac crystals and Fermi gasses.
We are able to solve the phonon line width, Eq.(31) numerically for different 2D Dirac crystals by performing the summation over ψ. To analyze the correlation between the phonon line width and the phonon wave number we derive γ ψ for different values of ψ within the range of 0.2 < ψ < 1. Fig. 7(a) shows that there is a linear correlation between γ ψ and ψ, with the linear fit having the coefficient of determination R 2 = 0.99. We, therefore, conclude that similar to the linear proportionality between the spectral energy and the fermion wave number, E ∝ k, there is also a linear correlation between the phonon line width and the phonon wave number, γ q ∝ q, in 2D Dirac crystals. We further compare the phonon line width approaching the FSN condition of 2D Dirac crystals at different temperatures in the static scenario. By studying Fig. 7(b) we observe a linear increase in the phonon line width or a linear decrease in the lifetime of the phonon mode as we increase the temperature. This suggests that at higher temperatures where the lattice vibrations increase the lifetime of the phonon mode decreases. Our results further match the experimental data collected on germanium [63] which further confirms our results.

Conclusion
In this paper by using the Lindhard model we derive an analytical expression for the dielectric response function of electrons strongly correlated to acoustic phonons in 2D crystals. The dielectric response function derived is a function of the phonon wave number which helps us better understand phenomenons in 2D Dirac crystals such as electron-phonon interactions and phonon self-energy. We see that different from free-electron systems, in 2D Dirac crystals the dielectric response function exhibits a cuspidal point near the FSN region in the static case where the phonon energy has been put equal to zero. We further observe when the phonon energy has not been put equal to zero the dielectric response function of a 2D Dirac crystal varies strongly near the FSN region even when the speed of sound is small, albeit nonnegligible, over the Dirac-electron Fermi velocity. We, therefore, show that the common approach for calculating the dielectric response function of 2D Dirac crystals by putting the energy of the phonons equal to zero is not accurate. We also show the generality of our solution for different 2D Dirac crystals with different Fermi velocities. We further use the dielectric response function obtained for 2D Dirac crystals to derive their phonon line width. We observe that the phonon line width increases as we move towards the FSN and as we go towards higher temperatures which matches the experimental work presented in the literature. *** A Appendix: Detailed calculation of the second addend of the dielectric response function The second addend of the dielectric response function at zero temperature for ψ < 1 is: The same derivations can be made for ψ > 1.