A finite electric-field approach to evaluate the vertex correction for the screened Coulomb interaction in the quasiparticle self-consistent GW method
Hirofumi Sakakibara Takao Kotani Masao Obata Tatsuki Oda





Department of Applied Mathematics and Physics, Tottori university, Tottori 680-8552, Japan


Department of Computational Science, Institute of Science and Engineering, Kanazawa University,
Kakuma, Kanazawa 920-1192, Japan



We apply the quasiparticle self-consistent GW method (QSGW) to 
slab models of ionic materials, LiF, KF, NaCl, MgO, and CaO, under electric field.
Then we obtain the optical dielectric constants
ϵ_∞( Slab) from the differences of the slopes of
the electrostatic potential in the bulk and vacuum regions.
Calculated ϵ_∞( Slab) show very good
agreements with experiments. For example,


we have  ϵ_∞( Slab)=2.91 for MgO, in agreement with
the experimental value ϵ_∞( Experiment)=2.96.
This is in contrast to ϵ_∞( RPA)=2.37, which is calculated in 
the random-phase approximation for the bulk MgO in QSGW. 
After we explain the difference between the quasiparticle-based perturbation
theory and the Green's function based perturbation theory,
we interpret the large difference ϵ_∞( Slab)-ϵ_∞( RPA)=2.91-2.37
as the contribution from the vertex correction of the proper
polarization which determines the screened Coulomb interaction W.
Our result encourages the theoretical development of self-consistent
G_0 W approximation along the line of QSGW self-consistency, as
was performed by Shishkin, Marsman and Kresse
[Phys. Rev. Lett. 99, 246403(2007)].
 













71.10.-w, 71.15.-m, 71.15.Dx

A finite electric-field approach to evaluate
the vertex correction for the screened Coulomb interaction
in the quasiparticle self-consistent GW method
    Tatsuki Oda
    December 30, 2023
=======================================================================================================================================================









§ INTRODUCTION


The quasiparticle self-consistent GW (QSGW) is one of the most
reliable method to determine the one-particle effective Hamiltonian
which describes the independent-particle picture,
or the quasiparticle (QP) picture,
for treating electric excitations of materials <cit.>. 
Other competitive methods such as HSE <cit.> and 
Tran-Blaha-09 functional <cit.> 
may work well in many systems, although we may need to use
material-dependent parameters <cit.>.
In contrast, QSGW is virtually parameter free and gives reliable descriptions for
a wide range of materials, not only metals and semiconductors, but also 
transition-metal oxides, type-II superlattice, and 4f systems
<cit.>.
Since heterogeneous mixtures of materials are used in current technologies,
QSGW is worth to be developed more as a tool to treat electronic
structures of such materials, where methods including such material-dependent parameters are hardly
applicable.

However, QSGW as it is has a shortcoming that it gives a systematic
overestimation of the exchange effects. This results in 
a little larger band gaps in QSGW for materials.
In fact, Faleev, van Schilfgaarge, and Kotani <cit.>
made a suggestion that the overestimation is removed if we perform
improved QSGW calculations taking into account the enhancement
of the screening effect due to the electron-hole correlation
in the evaluation of the screened Coulomb interaction W.
This is based on the theoretical consideration combined with the observation
that the calculated optical dielectric constant ϵ_∞ in the
random-phase approximation (RPA) in
QSGW gives ∼ 20 percent smaller ϵ_∞
than experiments for kinds of materials <cit.>.

Such an improved calculation was performed by Shishkin, Marsman and
Kresse, where they include the enhancement of the screening effect <cit.>.
The enhancement is via the vertex correction for the proper polarization P, which
determines W=v/(1-vP), where v denoted the Coulomb interaction.
They approximately include the lowest-order vertex correction 
due to the electron-hole correlation; see Eq.(15) around in Ref.Bruneval05.

Their results are theoretically quite satisfactory in the sense
that both band gaps and ϵ_∞,
which are calculated simultaneously and self-consistently without
parameters as in HSE, are in agreement with experiments.
For example, calculated values ϵ_∞=2.96
and band gap E_ G=8.12eV for MgO are in agreement with the experiments, 2.95, and 7.83
eV, respectively. See scGW(e-h) in Table I in Ref.shishkin_accurate_2007.
Furthermore, based on these theoretical analyses, we can introduce QSGW80
to avoid the very expensive computational costs of the method by
Shishkin et al.. QSGW80 is just a simple hybridization, 80 % QSGW+ 20 % GGA
to include such enhancement of the screening effectively.
The hybridization is very different from the hybridization in HSE,
where the mixing ratio α between GGA and the Hartree-Fock strongly affects to final results.
QSGW80 works well to describe experimental band gaps <cit.>.
The performance of QSGW80 is systematically examined in
Ref.deguchi_accurate_2016 by Deguchi et al.,
where we see both the calculated band gaps and effective masses are in good
agreements with experiments. QSGW80 is successfully used for practical applications, for example,
to the type-II superlattice of InAs/GaSb
<cit.>.

In this paper, we evaluate ϵ_∞ not in bulk calculations
with such approximations used in Ref.shishkin_accurate_2007,
but by slab models with finite electric bias voltage.
We treat five ionic materials, LiF, KF, NaCl, MgO and CaO.
We put a slab in the middle of vacuum region in a supercell.
The electric field is applied by the effective screening medium
(ESM) method given by Otani and Sugino <cit.>.

We obtain ϵ_∞ from the ratio of slopes of the electrostatic
fields at slab region and at vacuum region.
Our approach is based on the self-consistent method, thus we do not need to
utilize approximations as was used in Ref.shishkin_accurate_2007.
Since we explicitly treat the response to the bias,
our method includes higher-order effects in a self-consistent manner.



Our findings are that the calculated ϵ_∞ in QSGW for the
slab models are very close to experimental values.
This is in contrast to the fact that ϵ_∞ in RPA of QSGW
are generally ∼ 20 percent smaller than experimental values.
This indicates that the vertex correction at the level of derivative of the QSGW self-energy
should make W be in agreements with experiments.
Our results is consistent with Table II in Ref.shishkin_accurate_2007.


We can interpret the enhancement of screening, represented by the enlargement of ϵ_∞,
as the size of the vertex correction for the proper polarization P.









Note that the vertex correction we evaluate is not what is defined in
the Hedin's equation <cit.>.
In the equation, we see P= -i G G Γ, that is,
the vertex function Γ is for the correction to P= -i G G, 
where G denotes the Green's function.
Instead, we rather evaluate Γ for P=-i G_0 G_0 Γ,
where G_0 is the bare Green's function. 
To clarify the above theoretical point on Γ, we give an extensive discussion in Sec.<ref>.
We explain role of Γ in the two kinds of perturbation theories.
In Sec.<ref>, we explain QSGW+ESM, an implementation of QSGW
combined with ESM. The QSGW+ESM for slab models should be very useful
not only for our purpose here, but also for others where usual GGA+ESM have difficulties.
In Sec.<ref>, we show our results of ϵ_∞.
Then they are interpreted as the vertex correction. In
Sec. <ref>, we give a
rationale of QSGW80, followed by a summary.













§ QP-BASED PERTURBATION VS. G-BASED PERTURBATION


To make our motivation in this paper clarified, we have to clarify
the difference between the quasiparticle-based perturbation (QbP) and
the Green's function-based perturbation (GbP). 
QbP is based on the Landau-Silin's QP theory,
while GbP is on the Hedin's one. 
To illustrate the difference between QbP and GbP,
we give a narrow-band model as follows. The model represent situations where we have good
QP picture (= independent-particle picture).


In advance, remind that we mainly have two kinds
of excitations in the paramagnetic electronic systems.
That is, the multi-particle excitations, 
and the collective excitation such as plasmons.
The former is described by QPs interacting each other.
Note that plasmons is located at high energy
because of the long-range Coulomb interaction <cit.>.


These excitations can be hybridized. For example, we know pseudo plasmon
in Silver, where one-particle excitations of 3d electrons are hybridized with plasmons.




 §.§ narrow-band model to explain the QP-based perturbation


The QPs based on the Landau-Silin's Fermi liquid theory is originally for
metals <cit.>. However, the idea of QPs are
rather easily applicable to insulators. We can consider QbP based on the QPs.
To illustrate this, let us consider a narrow-band model, a paramagnetic case given by a
Hamiltonian H, which has an one-body term represented by
finite numbers of Wannier functions in the primitive cell,
and the Coulomb-like interaction e^2/ϵ' |-'|,
where ϵ' is a constant.
We consider a case that it gives the QPs shown in Fig.<ref>
given by the one-particle Hamiltonian H_0.



Based on the perturbation, we expect that the lifetimes of all the electrons are infinite
because the band gap E_ G is large enough to forbid all the
electrons decaying into lower energy electrons accompanying electron-hole pairs, holes as well.
In other words, all bands are within the threshold of impact ionization <cit.>.
Thus QPs described by G_0=1/(ω-H_0) should give well-defined
one-particle excitations of the narrow-band model.

This conclusion should be essentially kept even when we fully turn on the
interaction as long as the following conditions are well satisfied.
First, excitons, binding states of electron-hole pairs,
should be only slightly at lower energies than E_ G.
Second, plasmons should be located at high enough energies 
so that the plasmon are hardly hybridized with the one-particle excitations.
Under these assumptions, we have well-defined QPs.
We can consider a path of adiabatic connection with keeping the
well-defined QPs given by H_0 since the QP spectrum
is clearly separated from the other excitations.
It can be written as H_λ =H_0+ (H_λ-H_0) for λ is from
zero to unity, where H=H_λ=1. Note that (H_λ-H_0) need to contain
λ-dependent one-body term.



In this narrow-band model, QbP should be suitable;
QPs are interacting each other by (H_λ=1-H_0).
As for the proper polarization function, we can take all the
non-interacting two-body QP excitations correctly by P_0= -i G_0 G_0.
Based of this QbP, G_0W_0 approximation, Σ=i G_0  W_0 where W_0=v/(1-vP_0),
is physically justified. That is, it describes how the motion of QPs in G_0 is perturbed by the
dynamical self-interaction given by W_0 in RPA.
Most of all first-principles calculations in literatures are implicitly
based on this QbP, while some latest literatures <cit.>
are based on GbP explained in Sec.<ref>.
QSGW is a method to determine H_0 self-consistently in QbP.




 §.§ vertex function  Γ_ QbP in the QP-based perturbation

To improve G_0W_0, we may include electron-hole correlation.
The corrections replace P_0= -i G_0 G_0 with P where
we include the correlation via the Bethe-Salpeter equation (ladder
diagrams). That is, we include two-body spectrum in the proper polarization
accurately. This lets us to use W=1/(1-vP) instead of W_0, resulting the G_0 W
approximation as Σ=i G_0 W.
In this paper, we concentrate on the G_0 W approximation
as in the case of scGW(e-h) in Ref.shishkin_accurate_2007.
In QbP, we thus define the vertex function Γ_ QbP for W
as P= -i G_0 G_0 Γ_ QbP. Roughly speaking, P=P_0 Γ_ QbP.
We give the numerical evaluation for this Γ_ QbP via the
evaluation of ϵ_∞ as shown in Sec.<ref>.


To go beyond Σ=i G_0 W approximation here, 
we need to take into account three-particle intermediate states.
However, it is not theoretically straightforward because of
a double counting problem that Σ=i G_0 W
already partially takes into account such states.
We will need to construct theories of
three-particle problem without double counting
along the line of first-principles calculations.
We do not treat this problem in this paper.














 §.§ G-based perturbation and vertex function Γ_ GbP


Let us consider how we can apply GbP to the narrow-band model.
In contrast to G_0 in QbP, the one-body Green function G has
complex meanings. We have imaginary part of G at high energies,
representing QPs hybridized with plasmons (plasmarons).
Because of sum rule for Im[G], the QP parts are suppressed by Z-factor.
Thus there is a problem that P= -i G G do not contain the two-body non-interacting
excitations with the correct weight, in contrast to the case P_0=-i G_0 G_0.

In principle, this problem is corrected by including the vertex function
Γ_ GbP in the Hedin's equation 
to determine the one-particle Green's function G(1,2) <cit.>.
Because Hedin's equation is theoretically rigorous, we expect
-i G G Γ_ GbP≈ P_0=-i G_0 G_0 in the model, under the
discussion of Sec.<ref> that P_0 gives good approximation for the model.
That is, contributions related to the collective excitations and
renormalization factors Z in -i G G should be virtually taken away
by the factor Γ_ GbP.
However, such numerical calculations should be computationally very
demanding <cit.>.


Similar discussion of Z-factor cancellation is also seen when we
multiply G_0 to W. That is,
we should have G_0 W ≈ G W Γ_ GbP since QbP correctly treats the model. 
Our analysis here is consistent with Takada's analysis
based on the Ward identity <cit.>.

QbP should be generally superior to GbP even in real materials.
In contrast to GbP, QbP is quite simple and physically convincing.
We should not be confused with the similarity of QbP and GbP.







In this paper, we evaluate Γ_ QbP.
In the following, we calculate the enhancement of the screening effect.
Then we evaluate the size of the ratio P/P_0
from the comparison between calculated ϵ_∞ in RPA and ϵ_∞
in the slab models. This ratio gives the size of Γ_ QbP.















§ QSGW COMBINED WITH EFFECTIVE SCREENING MEDIUM METHOD


To calculate ϵ_∞, we put a slab
in the middle of vacuum region in a supercell.
ϵ_∞ is calculated from the difference of slopes
in the vacuum region and in the slab region under small bias voltage.
The supercell we use are detailed at the beginning of Sec.<ref>.
In such calculations, we can obtain ϵ_∞ beyond the bulk
calculation in RPA as we explain in the next paragraph.
That is, we can obtain ϵ_∞ including the effect of the vertex correction.









To explain how the effect is included, let us consider slab calculations in the case of GGA at first.
We first perform self-consistent calculation under zero bias.
Then we perform self-consistent calculation under the finite bias
(theoretically, it should be infinitesimally small).
Then we have the difference of the electron density
δ n() between the two calculations.
Simultaneously, we have a corresponding response of the one-particle
potential given as
δ V() = ∫ d^3r' v(-') δ n(') +
∂ V^ GGA_ xc/∂ n()δ n().
The last term is the difference in the exchange-correlation (xc) potential
caused by δ n() self-consistently.
That is, the derivative δ V()/δ n(') contains
the contribution of the xc kernel f_ xc=∂ V^ GGA_ xc/∂ n().

Under the bias, we can obtain ϵ_∞, from the ratio of
slopes of the electrostatic potential in the vacuum region and in the slab region.

It should contain the contribution from f_ xc, which is identified as the
vertex correction in GGA.



This is essentially the same in QSGW.



Recall that the self-energy in QSGW denoted as
V^ QSGW_ xc(,') is a static non-local potential,
replacing V^ GGA_ xc.
The derivative of the one-particle potential is given as
δ V(,') = ∫ d^3r” v(-”) δ n(”) +
 δ V^ QSGW_ xc(,'), where the last term play
a role of ∂ V^ GGA_ xc/∂ n()δ n().
Note that δ V^ QSGW_ xc(,') is determined
self-consistently, although it is not so simply given as
∂ V^ QSGW_ xc/∂ n()δ n().




Our calculations include the contribution of δ V^ QSGW_ xc(,')
self-consistently as in the case of GGA.
Our method is on the same spirit of solving
the Bethe-Salpeter equation in Ref.PhysRevLett.122.237402.

Our QSGW+ESM is implemented in a first-principles package ecalj
<cit.> which is based on a mixed-basis
method, the augmented plane wave (APW) and Muffin-tin (MT) orbital method (the PMT method)
<cit.>.






The PMT method is an all-electron full potential method which uses not
only the APW basis in the LAPW method, but also the MT orbitals in
the LMTO method simultaneously in the expansion of eigenfunctions.
It also use the local orbital basis <cit.>.
On top of the PMT method, we had implemented the QSGW method <cit.>.
In PMT, we use very localized untuned MTOs which contains damping factor
exp(- κ r), where κ are fixed
to be 1/bohr and/or 2/bohr, together with low-cutoff APWs (≤3 Ry).
We do not need empty spheres since the APWs can handle vacuum regions of slab models.
The charge density is represented in the three component representation,
'smooth part','true part within MT', and 'counter part within MT' as in
the case of PAW method <cit.>.
In contrast to the other GW methods which requires the
Wannier-interpolation technique to make band plots in
the whole Brillouin zone, we can make band plots easily without
resorting to the technique <cit.>.
In the following, we show how to implement ESM in
the PMT method, after an explanation of general theory of ESM.



 §.§ The electrostatic potential in the Effective screening medium method



We apply the ESM method <cit.> to slab
models under an external electric field.
We treat a supercell with periodic boundary condition
where we have a slab with periodicity in the xy-plane. The slab is at the middle of supercell.
Position in the cell is specified by =(_//,z). 
Planes at z=-z_0 and at z=z_0 are the left and right ends of the
supercells. The electrostatic potential is calculated from the charge
density in the supercell assuming two electrodes are at z=± z_0
(we set z_1=z_0 in Fig.1 of Ref.otani_first-principles_2006.)
for applying voltage to the supercell.
As we summarize as follows, the ESM in DFT is formulated from the total energy minimization,
however, it is not true in QSGW since QSGW itself is not formulated from
the total energy minimization. After we obtain the following key equation
eq:vone to determine electrostatic potential.
We use it even in QSGW.







Let us start from the energy functional of DFT in the ESM. It is written as

    E[n] = E^ kin[n] +E^ xc[n] + E^ es[n] + E^ app[n].

Here, we have kinetic energy E^ kin[n], xc
energy E^ xc[n], and the electrostatic energy E^ es[n] terms. In addition,
the last term is the applied electrostatic
term E^ app[n]=∫ d^3r V^ app() (n()+n_ N()),
where n() and n_ N() are the electron density
and the the charge density of nuclei, respectively; 
V^ app() is a linear function of z, representing the external field.


In ESM, we enforce the periodicity in the supercell for the electrostatic potential.
Thus we use V^ app()s() instead of 
V^ app(), where we introduce a support function s() 
which is unity for most of all regions, but is going to be zero at
z=-z_0 and z=z_0. It is different from unity only near the boundaries, 
z≈ -z_0 or z≈ z_0.
Thus the potential V^ app()s() recover the periodicity
of the supercell. A constant can be added to V^ app()
so that it keep smooth periodicity over z=± z_0.
As long as we use large enough vacuum region, we
have little electrons near the boundaries.
Thus the choice of s() is irrelevant.

A key in ESM is that we use the Green function v̅(,') for the
electrostatic energy E^ es[n] instead of the Coulomb interaction
v(-') in usual the GGA calculations.
As in Ref.otani_first-principles_2006,
v̅(,') contains not only the Coulomb interaction v(-') but
also the effects due to the polarization of virtual electrodes,
which are at z=± z_1 (we use z_0=z_1 in our calculations here).

Polarization of the slab occurs with keeping the electrostatic potential being constants at electrodes.
Corresponding to V^ app()s(),
we use s()v̅(,')s(') instead of v̅(,') in practice. 
Then we have well-defined Kohn-Sham total energy with keeping the
periodic boundary condition for given V^ app().

The minimization of E[n] with respect to n() gives the
Kohn-Sham potential V() as

    V()= ∫ d^3r' v̅(,') (n(')+n_ N(')) +  V^ app() + V^ xc(),

Hereafter, we skip s() for simplicity.


In QSGW <cit.>, we cannot derive its
fundamental equation from the energy minimization.
Thus the formulation of QSGW+ESM is not exactly along the line above.
However, we can use the one-particle potential of eq:vone in the
self-consistent cycle, where V^ xc() is replaced by a
static version of the self-energy <cit.>.
Thus, in principle, it is straightforward to perform QSGW+ESM.





 §.§ ESM in the PMT method

In the PMT method, electron density (and also the charge density) 
is represented by the three component formalism described in
Ref.kotani_formulation_2015,
originally introduced by Soler and Williams <cit.>.

At first, space is divided into MT regions and interstitial regions.
Then electron density is represented by three components as
n={n_0(),{n_1,a()},{n_2,a()}} where a is the
index of atomic sites in the primitive cell. 
Following Ref.kotani_formulation_2015,
this is simply expressed as n=n_0 ⊕ n_1 ⊖ n_2.
The 0th component n_0() is the spatially smooth functions,
expanded in analytic functions, that is, planewaves, Gaussians, and smooth Hankel functions
<cit.>.
The 1st components n_1,a()
is the true electron density within MT at _a.
The 2nd components n_2,a() is the counter part,
that is, the projection of n_0() into the MT at _a.
n_0() and n_2,a()  are identical within MT at _a up
to given angular momentum cutoff in their spherical harmonics expansion.

We can get all charge density 
n^ Zcv=n^ Zcv_0 ⊕ n^ Zcv_1 ⊖ n^ Zcv_2
by adding the ion-core contribution to n.
Then we apply the multipole transformation clearly defined in Ref.
kotani_formulation_2015, resulting 
n̅^ Zcv_0 ⊕n̅^ Zcv_1 ⊖n̅^ Zcv_2 as shown in Eq.(28-30) in Ref.kotani_formulation_2015.
The transformation makes
n̅^ Zcv_0(), n̅^ Zcv_1,a(), and n̅^ Zcv_2,a() 
have the same multipole in each MT site at _a, although physically observable density unchanged.




The 1st components n̅^ Zcv_1,a(), unchanged by the transformation, are  
the sum of ion-core charge density and n_1,a().


From the smooth density n̅^ Zcv_0(), we can calculate
electrostatic potential as
V^ es_0()= ∫ d^3r' v̅(,') n̅^
Zcv_0(') +  V^ app(). 
This gives correct interstitial part of the potential V^ es_0()
calculated from all charge density. The values of V^ es_0()
at MT boundaries are used to determine the electrostatic potential within MTs.



We can use usual procedure to determine the electrostatic potential within MTs.
In each MT, we have 1st and 2nd components
n̅^ Zcv_1,a() and n̅^ Zcv_2,a(), which have the same multipole. 
With the condition that the electrostatic potential is zero at the MT
boundary, we can calculate the potential generated by the difference of
the 1st and 2nd components.

Thus we finally have the electrostatic potential V^ es() represented in the three component formalism.
With this potential, we can perform self-consistent calculations for slab models.






















§ RESULTS




 §.§ Optical dielectric constants via the slab model




In Fig.<ref>, we illustrate our treatments 
in the slab models for five NaCl-structure ionic materials, 
where we use ± z_0=± 30 a.u.
We use slabs made of nine layers, 18 atoms in the supercell.
We use experimental lattice constants of bulk materials, without relaxation of atomic positions.
The electrostatic potential V_z^ es(z,E) are
the average of V_0^ es() in the xy plane under the bias voltage E.
We plot the cases of E=0.2 Ry and of E=0.0 Ry.
We show Δ V_z^ es(z)=V_z^ es(z,0.2  Ry)-V_z^ es(z,0.0  Ry)
in the bottom panel in Fig.<ref>. 
Δ V_z^ es(z,E) changes linearly as a function of z in the vacuum
regions and in the slab regions.
From the ratio of two slopes of Δ V_z^ es 
in the slab region and in the vacuum region, we obtain
ϵ_∞( Slab).


Our main results are ϵ_∞ calculated from slab models in QSGW,
ϵ_∞( QSGW,Slab), in Table <ref>.
Note that ϵ_∞( Slab) contains the effect of vertex
corrections based on QbP (See Sec.<ref>), because changes of the self-energy caused by the bias E
are self-consistently taken into account.
Numerical reliability of our calculations


are estimated to be ≲ 1 percent.
See supplemental materials for computational details <cit.>.
In Table <ref>, we also show bulk values
ϵ_∞( RPA). To obtain them, we first perform
self-consistent calculations in QSGW for bulk materials.
Then we calculate ϵ_∞ in the random-phase
approximation (RPA) with/without local field correction (LFC).
We also show ϵ_∞ in GGA together.












The QSGW values are in good agreements with experiments.
For example, ϵ_∞( QSGW,Slab)=1.94 for LiF gives surprisingly
good agreement with ϵ_∞( Experiment)=1.96. 
In contrast, ϵ_∞( QSGW,RPA)=1.67 is very smaller than ϵ_∞( QSGW,Slab)=1.94.
These are generally true in all other materials.
We see that ratios η=ϵ_∞( QSGW,RPA)/ϵ_∞( QSGW,Slab)
in Table <ref> are ∼0.8. 
This is consistent with Ref.kotani07a
where ϵ_∞ for ZnO, Cu_2O, MnO, and NiO are presented.
From a point of view to estimate the enhancement factors
(≈ vertex Γ ) of the proper polarization, we may consider ratios
γ=ϵ_∞( Slab)-1/ϵ_∞( RPA)-1.
As shown in Table <ref>, γ∼1.4. 
Since ϵ_∞( QSGW,Slab) gives very good agreements with ϵ_∞( experiment),
we can say that the vertex correction for bulk should give the
difference between ϵ_∞( RPA) and ϵ_∞( experiment) very well, where
the vertex correction is calculated at the level of the functional
derivative of the self-energy in QSGW; See Sec.<ref>.



This is in contrast to the case of GGA. For example, look into the case
of LiF. The difference ϵ_∞( GGA,Slab)-ϵ_∞( GGA,RPA)=2.01-1.95=0.06 is
very small. The difference is originated from the xc kernel f_ xc in the density
functional perturbation theory. This is consistent with results in Ref.botti_long-range_2004 
where they explicitly evaluate f_ xc in GGA for bulk materials.
Note that ϵ_∞( GGA,Slab)=2.01 is a little larger than
ϵ_∞( experiment)=1.96; this is true for all other materials.
We see that the contributions of vertex corrections f_ xc do not necessarily improve
agreements; ϵ_∞( GGA,Slab) give poorer
agreement with ϵ_∞( experiment) than ϵ_∞( GGA,RPA).































 §.§ Rationale for QSGW80



Our result in Sec.<ref> shows that the vertex correction
should be included in the proper polarization P to obtain
ϵ_∞ in agreement with experiments.
We have to use such P in the QSGW self-consistent cycle.
Such improved QSGW self-consistency can be identified as a
self-consistent method in the G_0W approximation on the basis of QbP.
Ref.shishkin_accurate_2007 by Shishkin et al. gives a method on this idea.

However,


their method is too expensive for computational efforts to apply wide-range of materials.
In fact, although their method was applied to calculate ionization potentials
in Ref.PhysRevLett.112.096401, it was not
so satisfactory because calculations are performed in the combination of
simple materials (bulk calculations) with the supercell calculations in GGA.
Furthermore, no papers available to treat transition-metal oxides such as LaMnO_3 in their method.

We have to develop such an improved QSGW method applicable to wide range of materials.
Two requirements are the computational efficiency and the theoretical validity.

As a possibility to respect the efficiency, 


we can consider a hybridization method between QSGW and the density functional xc <cit.>.
In QSGW80, a simple hybridization, 80 % QSGW+ 20 % GGA, we can see that it
works well for wide range of materials. 
Our present results support the method of QSGW80, which takes only the
80 percent of QSGW self-energy. We can identify QSGW80 as a simplification
of the method in Ref.shishkin_accurate_2007.
Ref.PhysRevB.92.041115 also presents one another approximation at the level
of QSGW80 for the vertex correction in QSGW, resulting in similar good agreement with experiments.













Let us examine how QSGW80 is justified for materials calculated here.
This is by the fact that
η=ϵ_∞( QSGW,RPA)/ϵ_∞( QSGW,Slab) in Table <ref>
are approximately 80 %, and show little material-dependency.
Thus we expect that QSGW80 can mimic QSGW with the vertex corrections;
too large screened-exchange effect is reduced by the factor 0.8, with
adding 0.2 GGA term so as to keep the total size of the xc term.
In Table <ref>, we show band gaps in QSGW and QSGW80 for materials
treated here. The band gaps are systematically too large in QSGW
in comparison with experimental values <cit.>,
while QSGW80 gives rather better agreements with experimental values.
In Ref.deguchi_accurate_2016, we checked the performance of the QSGW80 for ranges of materials.
As in the case of Ref.shishkin_accurate_2007, QSGW80 is
theoretically reasonable in the sense that the band gaps are improved by using the corrected W.
To go beyond QSGW80, we have to develop methods to take the vertex
correction into W as was done in Ref.shishkin_accurate_2007 in a simple manner.
Considering the fact that QSGW80 works well as shown in Ref.deguchi_accurate_2016,
we may expect simple methods to represent the vertex correction
by a scalar factor or by limited number of parameters.
As long as we know, the vertex correction can be relatively insensitive
to materials, thus we expect some simple method might be available. 




§ SUMMARY

To clarify the importance of quasiparticle self-consistency in QSGW,
we have explained the quasiparticle based perturbation in Sec.<ref>.
Then we emphasize the importance of the self-consistency in the G_0 W
approximation. Then we obtain the quasiparticles (independent-particle) given by
H_0, and the interaction between the quasiparticles given by W.
The vertex correction in QbP is introduced.

We have performed QSGW calculations for slab models under electric field by means
of the ESM method. The calculated ϵ_∞ are in good agreements with experimental values.
Compared with ϵ_∞ in bulk calculation in RPA,
we evaluated the size of vertex corrections
as the functional derivative of the static self-energy in QSGW.
Our results on ϵ_∞ give a support to the method
by Shishkin, Marsman and Kresse <cit.>.
As a simplified substitution of their method, we examined the performance of
QSGW80 <cit.> for materials treated here.
The method QSGW+ESM developed for the calculations should be useful
even for other purposes such as bias-dependent
spin susceptibility in material theory, as well as practical device
applications and materials designs.

To go beyond usual QSGW, we should develop improved QSGW method
in the G_0 W approximation, whereas we should use accurate W by including the vertex correction.
Then we have virtually best division of H=H_0+(H-H_0) where H_0
gives the optimum independent particle picture.



T. K. thanks to supporting by JSPS KAKENHI Grant Number 17K05499.
We also thank the computing time provided by Research Institute for Information Technology (Kyushu University).
H. S. thanks to the computing resources provided by
 the supercomputer system in RIKEN (HOKUSAI) and the supercomputer system in ISSP (sekirei).




