Unexpected phonon frequency shift at Γ when including Quadrupole correction in EPW

[Hust_wty]

I am currently performing electron–phonon interaction calculations for wurtzite GaN using EPW, and I have encountered an issue related to the Quadrupole correction in the phonon self-energy calculation.

As shown in Figure 1

1.png (95.33 KiB) Viewed 290 times

， when the Quadrupole term is not included, the phonon self-energy calculation gives a reasonable result:
at the Γ (0 0 0) point, the phonon frequency ω = 0 meV, and both lambda and gamma___ are 0, which is physically correct.

However, as shown in Figure 2,

2.png (98.03 KiB) Viewed 290 times

when the Quadrupole correction is included, the phonon frequency at Γ starts from 18.3356 meV, which is clearly unphysical.
The Quadrupole tensor itself seems correct — it matches the literature and previously published results (see attached tensor tables).
Could this frequency shift be related to an incorrect phonon renormalization or normalization after including the Quadrupole correction?
the follow is epw input
--
&inputepw
prefix = 'gan'
amass(1) = 69.723
amass(2) = 14.0067
outdir = './output'

iverbosity = 3

elph = .true.
!epbwrite = .true.
!epbread = .false.
epwwrite = .fasle.
epwread = .true.
etf_mem = 1
vme = 'wannier'
!selecqread = .true.


lpolar = .true.
lifc = .false.
asr_typ = 'simple'
use_ws = .true.
lphase = .true.

nbndsub = 14
bands_skipped = 'exclude_bands = 1-12'

wannierize = .false.
num_iter = 1500
iprint = 3
dis_win_max = 24.0660
dis_win_min = 1
dis_froz_max= 11.03265
proj(1) = 'Ga:sp3'
proj(2) = 'N:p'
wdata(1) = 'bands_plot = .true.'
wdata(2) = 'begin kpoint_path'
wdata(3) = 'G 0.0000 0.0000 0.0000 M 0.5000 0.0000 0.0000'
wdata(4) = 'M 0.5000 0.0000 0.0000 K 0.3333 0.3333 0.0000'
wdata(5) = 'K 0.3333 0.3333 0.0000 G 0.0000 0.0000 0.0000'
wdata(6) = 'G 0.0000 0.0000 0.0000 A 0.0000 0.0000 0.5000'
wdata(7) = 'A 0.0000 0.0000 0.5000 L 0.5000 0.0000 0.5000'
wdata(8) = 'L 0.5000 0.0000 0.5000 H 0.3333 0.3333 0.5000'
wdata(9) = 'end kpoint_path'
wdata(10) = 'bands_plot_format = gnuplot'
wdata(11) = 'guiding_centres = .true.'
wdata(12) = 'dis_num_iter = 1000'
wdata(13) = 'num_print_cycles = 10'
wdata(14) = 'dis_mix_ratio = 1.0'
wdata(15) = 'conv_tol = 1E-8'
wdata(16) = 'conv_window = 4'

delta_approx= .false.
elecselfen = .false.
phonselfen = .true.
a2f = .false.

fsthick = 0.2 ! eV
temps = 300 ! K
degaussw = 0.02 ! eV

scissor = 1.5579

dvscf_dir = './save/'

efermi_read = .true
fermi_energy= 10.0328
!band_plot = .true.

!filkf = './eband.check'
filqf = './BTE.qpoints_36'

nq1 = 6
nq2 = 6
nq3 = 6
!nqf1 = 12
!nqf2 = 12
!nqf3 = 12

nk1 = 12
nk2 = 12
nk3 = 12
nkf1 = 60
nkf2 = 60
nkf3 = 60
!mp_mesh_k = .true.
/

28 cartesian
0.000000000000000E+00 0.000000000000000E+00 0.000000000000000E+00
0.000000000000000E+00 0.000000000000000E+00 0.102310586375119E+00
0.000000000000000E+00 0.000000000000000E+00 0.204621172750238E+00
0.000000000000000E+00 0.000000000000000E+00 -0.306931759125357E+00
0.000000000000000E+00 0.192450089729862E+00 0.000000000000000E+00
0.000000000000000E+00 0.192450089729862E+00 0.102310586375119E+00
0.000000000000000E+00 0.192450089729862E+00 0.204621172750238E+00
0.000000000000000E+00 0.192450089729862E+00 -0.306931759125357E+00
0.000000000000000E+00 0.384900179459723E+00 0.000000000000000E+00
0.000000000000000E+00 0.384900179459723E+00 0.102310586375119E+00
0.000000000000000E+00 0.384900179459723E+00 0.204621172750238E+00
0.000000000000000E+00 0.384900179459723E+00 -0.306931759125357E+00
0.000000000000000E+00 -0.577350269189585E+00 0.000000000000000E+00
0.000000000000000E+00 -0.577350269189585E+00 0.102310586375119E+00
0.000000000000000E+00 -0.577350269189585E+00 0.204621172750238E+00
0.000000000000000E+00 -0.577350269189585E+00 -0.306931759125357E+00
0.166666666666667E+00 0.288675134594792E+00 0.000000000000000E+00
0.166666666666667E+00 0.288675134594792E+00 0.102310586375119E+00
0.166666666666667E+00 0.288675134594792E+00 0.204621172750238E+00
0.166666666666667E+00 0.288675134594792E+00 -0.306931759125357E+00
0.166666666666667E+00 0.481125224324654E+00 0.000000000000000E+00
0.166666666666667E+00 0.481125224324654E+00 0.102310586375119E+00
0.166666666666667E+00 0.481125224324654E+00 0.204621172750238E+00
0.166666666666667E+00 0.481125224324654E+00 -0.306931759125357E+00
0.333333333333333E+00 0.577350269189585E+00 0.000000000000000E+00
0.333333333333333E+00 0.577350269189585E+00 0.102310586375119E+00
0.333333333333333E+00 0.577350269189585E+00 0.204621172750238E+00
0.333333333333333E+00 0.577350269189585E+00 -0.306931759125357E+00

To validate my setup, I also tested the official SiC Quadrupole example provided with EPW. the SiC calculation without Quadrupole gives ω = 0 meV at Γ (correct), and the computed scattering rates are reasonable.
When the Quadrupole correction is added, however, ω becomes –10.9810 meV at Γ, again unphysical.

I have checked Gaussian broadening, and all relevant parameters, but the results remain the same.

I have also attached a comparison of my calculated GaN electron–phonon linewidths with and without Quadrupole corrections.

10.png (82.38 KiB) Viewed 290 times

Could you please advise whether this issue could be due to a normalization problem in the phonon renormalization procedure with Quadrupoles, or if there is another known reason for this behavior?