Dielectric Function of silver

This tutorial offers a study of optical dielectric function ε(ω) in silver. The system is relatively simple, and it provides a useful test bed to show what is captured by different levels of theory, and what is left out. We wil use the optics branch of lmf.

Experimentally, some quantities are accurately known: the DC conductivity is σ(ω=0) = 6.3×10⁷ (Ω-m)⁻¹. Also it is well established that a sharp shoulder appears around 4 eV in the imaginary part of ε(ω).

Other quantities are also measured, but with less precision, as we discuss here.

Preliminaries

This tutorial uses several Questaal executables, e.g. blm, lmfa, lmf, pfit, plbnds, and fplot, and if you cover the GW part, the shell script lmgw.sh.

This tutorial assumes you have an 8-processor machine, with the Intel MKL library part of the compilation. We use the following shorthand as a launcher for jobs with 4 processes, and 2 threads.

mpix --omp=2 --mkl=2 -np=4 

you can run in scalar mode, and simply ignore the launcher. Your machine may require details this tutorial cannot anticipate. A typical launch instruction would be

env OMP_NUM_THREADS=2 MKL_NUM_THREADS=2 mpirun -n 4

In this tutorial we will compute the dielectric function within the RPA, mostly with local fields. The optics branch in lmf will calculate Im ε_xx without local fields

Theory

For a local potential and neglecting local field effects, Im ε reads [4] $$\mathrm{Im}\,\epsilon(\omega) = \sum_{vc} \sum_{\mathbf{k}} \left| \int{d^3\mathbf{r} \Psi_c^{\mathbf{k}*}(\mathbf{r}) (\hat{\mathbf{e}}\cdot\mathbf{p}) \Psi_v^{\mathbf{k}}(\mathbf{r})} \right|^2 \delta(E_c^{\mathbf{k}} - E_v^{\mathbf{k}} - \hbar\omega) f_{cv} \quad\quad [1]$$

Here c and v refer to unoccupied and occupied states, respectively, and $\Psi_c$ and $\Psi_v$ the corresponding one-particle wave functions with eigenvalues $E_c$ and $E_v$; $\hat{\mathbf{e}}$ is the unit vector, and $\mathbf{p}$ is the momentum operator.

For an insulator, c and v span the unoccupied and occupied states, respectively. In the metals case we consider here, bands cross the Fermi level so c and v are not cleanly separated. This is taken into account by including the difference in Fermi functions $f_{cv}=f_c-f_v$, and c and v now span all states.

If a nonlocal potential is added to the local one, the velocity operator should replace $\mathbf{p}/m$. However, if the nonlocal potential is diagonal in the band basis of the local potential, it does not change $\Psi_c$ and $\Psi_v$. In such a case it can be shown [4] that an effective momentum operator may be used instead, which is written as $$\left< c\mathbf{k} | \mathbf{p}^\mathrm{eff} | v\mathbf{k}\right> = \left< c\mathbf{k} | \mathbf{p} | v\mathbf{k}\right> \frac{E_c^{\mathbf{k}} - E_v^{\mathbf{k}} }{E_c^{0\mathbf{k}} - E_v^{0\mathbf{k}} } \quad\quad [2]$$

Energies in the numerator and denominator correspond to whether the nonlocal potential is included or not, respectively. You can use this correction factor by setting MEFAC=2 in OPTIONS. We will do this in the present tutorial, when optics are calculated using lmf with the QSGW self-energy.

In the metallic case we must distinguish intraband and interband transitions. The former is a transition with an infinitesimal energy changes from an occupied state to an occupied state of the same band at the Fermi level E_F, while the latter are transitions between different bands. lmf’s optics branch takes into account only interband transitions. When the two are separated the total dielectric function becomes $\epsilon \epsilon^{\mathrm{intra}} + \epsilon^{\mathrm{intra}} - 1$. The GW optics computes $\epsilon$ for a small but finite wave number q, and thus there is no separation.

$\epsilon^{\mathrm{intra}} = 1 - \omega_p^2/\omega^2$. $\omega_p$ is called the “plasma frequency”. It can be interpreted as the physical plasma frequency (frequency of collective oscillations), but the usual construction $\omega_p$ only takes into account the intraband contribution, and corresponds only approximately to the physical plasma frequency.

Traditionally $\omega_p$ is taken to be a scalar, but the intraband contribution was recently rederived [5] specifically to discuss its tensorial nature. It can be written as $$D^{ij} = \frac{e^2}{2\pi^2} \sum_{n} \int_\mathrm{BZ} {d^3\mathbf{k} v^i(\mathbf{k}) v^j(\mathbf{k})} \delta(E_F - E_n^{\mathbf{k}}) \quad\quad [3]$$

with the usual plasma frequency defined as $$\omega_p^2 = \frac{\mathrm{Tr} D/3}{4\pi}. \quad\quad [4]$$

lmf has the ability to compute $\omega_p^2$ in its tensorial form, as we show below. Since Ag is cubic, $\omega_p^2$ is a scalar, and the tensorial nature doesn’t matter here.

Scattering mechanisms are another important omission. Scattering can occur in part from electron-electron scattering, but also from phonons (phonons are usually the dominant ones). Scattering mostly modifies the intraband contribution. The intraband contribution has the same form as the Drude model. When this model includes a relaxation time $\tau$, $\omega^2\rightarrow{\omega(\omega + i/\tau)}$ [6]. We use the Drude model to modify the intraband contribution as $$\epsilon^{\mathrm{intra}} - 1 = - \frac{\omega_p^2}{\omega(\omega + i/\tau)}\quad\mathrm{which\ has\ an\ imaginary\ part}\quad \frac{\omega_p^2\,\tau}{\omega{\,\left(\omega^2\,\tau^2+1\right)}}. \quad\quad [5]$$

Calculation of $\tau$ is highly nontrivial. Here we will use an experimental value for it instead.

lmf computes Im ε only. However the real part can be derived from it via a Kramers Kronig transformation $$\chi(\omega)={\frac {1}{i\pi }}{\mathcal {P}}\!\!\int_{-\infty }^{\infty }{\frac {\chi (\omega ')}{\omega '-\omega }}\,d\omega'. \quad\quad [6]$$

The KK transformation connects the real and imaginary parts of any complex function that is analytic in the upper half-plane; the dielectric function and conductivity are two instances of this. Questaal has a code, kkt, which obtain the real part from the imaginary part by calculating the Hilbert transform shown in the preceding equation.

Finally, to compare to experiment we will use the well known relation between dielectric function and conductivity $$\epsilon = 1 + 4\pi{i}\sigma/\omega . \quad\quad [7]$$

Thus the real part of the intra-atomic (Drude) contribution to $\sigma$ is $$\mathrm{Re}\,\sigma^\mathrm{Drude} = \frac{\omega_p^2\,\tau}{4\pi\left(\omega^2\,\tau^2+1\right)} \quad\quad [8]$$

Tutorial

Ag assumes an fcc sublattice. Copy the box below to init.ag, which takes the room temperature lattice constant.

LATTICE
   ALAT=7.720 PLAT= -0.0 0.5 0.5  0.5 -0.0 0.5  0.5 0.5 -0.0
SITE
   ATOM=Ag  X=0 0 0

Create a ctrl file and a symmetry lines file

blm --brief --rsham~el=-.4,-1.2 --upcase --nfile --optics  --gw~emax=3~nk=8 --nk=12  ag
sed -i.bak s:'NPTS=1001 WINDOW=0,1':'NPTS=3001 WINDOW=0,5 DW=0.002,0.02':g ctrl.ag

The first line makes a ctrl file. The second line modifies tokens in the OPTICS category to span a wider window than the default values blm supplies.

Get the free atom density and make Ag self-consistent at the LDA level

lmfa ag 
mpix --omp=2 --mkl=2 -np=4 lmf ag 
cp rst.ag rst.lda

Next, calculate the optics with lmf. Note in the ctrl file, the MODE token is controlled by variable loptic, which makes it possible to turn on the optics part through the command line

mpix --omp=2 --mkl=2 -np=4 lmf ag -vnkabc=24 -vloptic=1 --quit=rho

This should invoke lmf’s dielectric function maker. You should see this block in the output

 OPTINT:  eps2(w) using occ bands (1,30) and unocc (1,30)
          3001 points in energy window (0,5)
          Ef=-0.117188  Treating bands as metal
          Integration with tetrahedron method, mode 1

and opt.ag should be generated. It has four columns: frequency (in Ry), followed by Im ε_xx, Im ε_yy, and Im ε_zz. Ag is cubic so they are all equivalent.
Plot Im ε_xx with fplot -x 0,12 -ab 'x1*13.6' opt.ag, or use your favorite graphics package. The fplot instruction converts the x axis to eV. and generates the figure shown.

Experimental data can be found in Refs [1,2,3]. We will use Ref. [3] as the benchmark because it is the most recent, and the samples are thought to be better.

Several observations are noteworthy here:

• Im ε vanishes identically below 2  eV, contrary to experiment; see Ref. [3], Fig. 4. Since the DC conductivity σ(ω=0) is related to Im ε(0), it implies σ(ω=0) vanishes. This is because the scattering mechanisms (especially electron-phonon scattering) have been left out. Remarkably, scattering is needed if the D.C. conductivity is not to vanish.
• A sharp shoulder begins a little below 3 eV with the midpoint of the rise close to 3 eV. This shoulder is observed experimentally, with the midpoint at ∼4 eV. The shoulder originates from transitions between the filled Ag d bands and the sp like bands around and above E_F. The rise starts too soon in the calculation, because the LDA puts the Ag d states too shallow. This will be revisited in the QSGW section.
• The peak value of Im ε is about 6, and it falls off to 2 at about 7 eV. In Ref. [3] its peak value is ∼4.

kkt will yield Re ε from Im ε.

kkt opt.ag

It writes to kk.out by default, which contains Re ε and Im ε side by side. Inspect the file and observed that ε_∞=Re ε(ω=0)≈4.1. A precise experimental determination of this quantity is difficult because it is buried in the (divergent) Drude term. the scattering time $\tau$ and plasma frequency $\omega_p$ have been variously reported, with similar uncertainties.
For experimental benchmark we will use the following estimates reported in Ref. [3]: $\omega_p{=}8.9{\pm}0.2\,\mathrm{eV},\enspace\tau=17{\pm}3{\times}10^{-15}\,\mathrm{sec},\enspace\epsilon_\infty=5{\pm}2$.

Re ε(ω) diverges as $-\omega_p^2/\omega^2$, owing to the intraband term (see theory section). $\hbar\omega_p$ is also difficult to pin down precisely experimentally, but values between 8.5 and 9.5 have been published (see Table 1 of Ref.[5]). The expression for $\omega_p$ described in the theory section has been implemented Questaal; to calculate it, do the following

mpix --omp=2 --mkl=2 -np=4 lmf ag -vnkabc=72 --quit=rho
lmdos ag -vnkabc=72 --dos~wp

(A conservative k mesh was used to avoid the need to monitor convergence.)

You should see the following printed to stdout:

 Plasma tensor at Ef=-0.117834 Ry  (Tr(wp^2)/3)^1/2 = 9.098 eV
              Dij (Ry^2)                           Dij (eV^2)
   0.447120   0.000000   0.000000        82.76848    0.00000    0.00000
   0.000000   0.447120   0.000000         0.00000   82.76848    0.00000
   0.000000   0.000000   0.447120         0.00000    0.00000   82.76848

Thus lmf predicts a very reasonable value for $\hbar\omega_p$. It is slightly less than that reported in Ref. [3], presumably because of the pseudopotential approximation they use.

kkt also implements the intraband term modified by $\tau$ as describe in the theory section. Invoke it as follows

kkt -outfile=eps.lda -units:out=ev -wp=9.1/13.6,235 -cnvtmid -skip0=2 opt.ag

-wp=9.1/13.6,235 tells kkt to add a Drude term with $\hbar\omega_p{=}9.1/13.6$ Ry (you must supply this quantity in Ry). The second number is a dimensionless quantity $\omega_p\tau{=}$: note that $\omega_p\tau{=}9.1\,\mathrm{eV}(6.58{\times}10^{-16}\,\mathrm{eV{-}sec})^{-1}17{\times}10^{-15}\,\mathrm{sec}=235$. -cnvtmid tells kkt to treat the frequencies as a histogram with the response function the value at the center of each histogram. It doesn’t matter much, but for small $\omega$ it makes a little difference. Finally, since the KK algorithm has difficulty stabilizing the transformation very near the origin, the first two energies are omitted (-skip0=2).

Inspect eps.lda. Comparing to Fig. 6 of Ref. [3], the agreement is seen to be very good. Both fall off as $\omega^{-3}$ (in this regime $\omega\tau$ is already much larger than 1). At $\omega=0.1\,\mathrm{eV}$ and $\omega=1.0\,\mathrm{eV}$, $\mathrm{Im}\,\epsilon{\sim}2800$ and 3, respectively. These are essentially the same values in Fig. 6 of Ref. [3], confirming that the Drude model works well. (To draw a figure that compare with Fig. 6, try something like fplot -frme:lxy 0,1,0,1 -x .05,1 -y 1,10000 -ord x3 eps.lda.)

To compare to the DC σ we can read off $\sigma(\omega=0)={\omega_p^2\,\tau}/4\pi$ from Eq. [8] in the theory section. In atomic units, $\sigma_0 = {\hbar\omega_p\,\omega_p\tau}/{4\pi} = 9.1/13.6{\times}235/4\pi\,\mathrm{Ry} = 12.5\,\mathrm{Ry}.$ $\sigma_0$ is usually expressed in SI units, ohm^-1 m^{^-1}. To convert, use the conductance quantum $G_0 = 2e^2/h = 7.748{\times}10^{-5}$ ohm^-1. Also, $G_0/a_0 = 2e^2/ha_0 = 2e^2/a_0/(2\pi\hbar) = 0.6366\,\mathrm{sec}^{-1}.$ This implies $1\,\mathrm{sec}^{-1} = 0.6366^{-1}{\times}7.748{\times}10^{-5}/0.5292{\times}10^{-10}$ (ohm-m)^-1 = $2.300{\times}10^{6}$ (ohm-m)^-1.
We finally obtain $\sigma_0 = 12.5{\times}2.300{\times}10^{6} = 2.9{\times}10^{7}$ (ohm-m)^-1. This is a factor of two smaller than the experimental value, $6.3{\times}10^{7}$ (ohm-m)^-1. The discrepancy means that at very low frequency some modification of the Drude model is needed. Ref. [3] folded the discrepancy into an energy-dependent relaxation time, approximately frequency independent for ${\hbar\omega}$ constant above 0.1 eV, but increasing by a factor of ∼2 approaching 0 from the flat region (see Fig. 10 of Ref. [3].)

How QSGW modifies the conclusions

We noted that the shoulder rises at too low a frequency because the LDA does a poor job at placing the Ag d states. In this section we will make the QSGW self-energy and revisit the preceding section

First, make the LDA bands

lmchk ctrl.ag --syml~n=50~lblq:G=0,0,0,L=1/2,1/2,1/2,X=1,0,0,W=1,1/2,0,K=3/4,3/4,0~lbl=LGXWGK
mpix --omp=2 --mkl=2 -np=4 lmf ag -vnkabc=24 --quit=rho
mpix --omp=2 --mkl=2 -np=4 lmf ag -vnkabc=24 --band:fn=syml
plbnds -range=-8,6 -fplot~scl=.7~sh~ts=1 -ef=0 -scl=13.6 -dat=lda bnds.ag

Plot the bands with fplot -f plot.plbnds, or use your favorite graphics package. The Ag d bands lie in the interval (−6,−3) eV.

Note that files bnds[12345].dat were created. We will use these files later, to compare against QSGW.

To run QSGW, we need to create the GWinput file. While we are at it, we can make the pqmap file, so the code runs faster.

lmfgwd -vsig=0 ctrl.ag --job=-1
lmfgwd -vsig=0 ctrl.ag --job=0 --lmqp~rdgwin --batch~np=8~nplm=16~npgwd=16~nodes=1~pqmap@fill=.8@plot~vanilla16

Now we are ready to run lmgw.sh to create the QSGW self-energy.

lmgw.sh --incount --sym --tol=1e-5 --split-w0 --mpirun-1 "env OMP_NUM_THREADS=2 MKL_NUM_THREADS=2 mpirun -n 1" --mpirun-n "env OMP_NUM_THREADS=2 MKL_NUM_THREADS=2 mpirun -n 4" ctrl.ag

Note:  you may need to modify --mpirun-1 and --mpirun-n depending on your architecture.

First, make and plot the QSGW bands

mpix --omp=2 --mkl=2 -np=4 lmf ag -vnkabc=24 --quit=rho
mpix --omp=2 --mkl=2 -np=4 lmf ag -vnkabc=24 --band:fn=syml
plbnds -range=-8,6 -fplot~scl=.7~sh~ts=1 -ef=0 -scl=13.6 -dat=qsgw bnds.ag

Plot the bands with fplot -f plot.plbnds, or use your favorite graphics package. You can confirm that the Ag d bands are deeper by ∼1 eV than the LDA, in the interval (−3,−7) eV. The figure on the left draws both the LDA bands (red) and the QSGW bands (blue), together with some ARPES data (circles).

We can repeat the other parts also. Excluding the Drude term

mpix --omp=2 --mkl=2 -np=4 lmf ag -vnkabc=24 -vloptic=1 --quit=rho
kkt opt.ag

shows that QSGW causes ε_∞ to drop from 4.0 to 3.3. Both agree with observed data within the experimental uncertainty.

The plasma frequency is also slightly modified.

mpix --omp=2 --mkl=2 -np=4 lmf ag -vnkabc=72 --quit=rho
lmdos ag -vnkabc=72 --dos~wp

$\omega_p$ increases from 9.1 to 9.5, slightly larger than the LDA result.

kkt -outfile=eps.qsgw -units:out=ev -wp=9.5/13.6,245 -cnvtmid -skip0=2 opt.ag

Not surprisingly, the low-frequency behavior of $\mathrm{Im}\,\epsilon$ is very similar, since the same Drude term is used for both and $\omega_p$ and $\tau$ are both similar. Over a broader window, the QSGW response function compares very well with experiment. The following plots $\mathrm{Im}\,\epsilon$ for both LDA (red) and QSGW (blue), and was used to generate the image shown:

fplot -x .5,10  -y 0,8 -lt 2,col=1,0,0 -ord x3 eps.lda -ord x3 -lt 1,col=0,0,1 eps.qsgw

You can also generate the reflectivity.

kkt -outfile=refl.qsgw -units:out=ev -wp=9.5/13.6,245 -cnvtmid -skip0=2 -out=refl opt.ag
fplot -x 0,5 refl.qsgw 

You should see a figure quite similar to that of Fig. 8 in Ref. [3], falling sharply from 1 slightly above 3 eV., reaching a minimum at just under 4 eV.

Comparing dielectric functions calculated by the GW code and by lmf

The QSGW code calculates the dielectric function for finite wave number q, with and without local fields. We seek the q→0 limit, which we accomplish approximately with a small value of q along the x axis. q and the frequency mesh are supplied in the OPTICS category of the ctrl file, with the tags EPSQ=0,0,0.015 DELRE= 0.002 0.02. Because q is nonzero, the optics will pick up the intraband term, but with no scattering, and also frquencies below q·v_F are not picked up; thus below some cutoff $\epsilon(\omega)$ should be ignored.

To obtain the optics, do something like the following (details will depend on your architecture)

touch meta bz.h5
rm -r meta *.h5
m1="env OMP_NUM_THREADS=8 MKL_NUM_THREADS=8 mpirun -n 1"
mn="env OMP_NUM_THREADS=2 MKL_NUM_THREADS=2 mpirun -n 4"
ml="env OMP_NUM_THREADS=1 MKL_NUM_THREADS=1 mpirun -n 8"
md="env OMP_NUM_THREADS=1 MKL_NUM_THREADS=1 mpirun -n 8"
mc="env OMP_NUM_THREADS=4 MKL_NUM_THREADS=4 mpirun -n 2"
lmgw.sh -vnkgw=24 --epsq --mpirun-1 "$m1" --mpirun-n "$mn" --cmd-lm "$ml" --cmd-gd "$md" --cmd-vc "$mc" --split-w0 --tol=2e-5 ctrl.ag

This script should generate eps.h5. Inspect it with h5ls eps.h5 and you should see this:

eps                      Dataset {1, 298}
eps_nlfc                 Dataset {1, 298}
freq                     Dataset {298}
nfreq                    Dataset {1}
nq                       Dataset {1}
q                        Dataset {1, 3}

To retrieve $\epsilon(\omega)$ in the same format as above, do the following

mcx '-ff12.6,100(f20.12)' -r:h5~id=freq eps.h5 -s13.6 -r:h5~z~id=eps eps.h5 -real -r:h5~z~id=eps eps.h5 -s0,-1 -real -ccat -ccat > dat-24

Here we used a q mesh of 24 divisions. You can check convergence in the mesh by changing -vnkabc=24.

Note that below 0.25 eV, $\mathrm{Im}\,\epsilon(\omega)$ becomes sharply positive, and $\mathrm{Re}\,\epsilon(\omega)$ turns around and becomes positive. Thus 0.25 eV is the cutoff just referred to. To remove this part, rework the mcx command:

mcx '-ff12.6,100(f20.12)' -r:h5~id=freq eps.h5 -s13.6 -r:h5~z~id=eps eps.h5 -real -r:h5~z~id=eps eps.h5 -s0,-1 -real -ccat -ccat -inc 'x1>0.25' > dat-24

Compare $\mathrm{Im}\,\epsilon(\omega)$ as generated by lmf to the QSGW result. For an apples-to-apples comparison we need to remove the scattering from the Drude term. We also skip the first few points where the QSGW result is inapplicable.

kkt -outfile=eps.qsgw0 -units:out=ev -wp=9.5/13 -cnvtmid -skip0=10 opt.ag

Draw a figure of $\mathrm{Im}\,\epsilon(\omega)$ calculated by : lmf with both $\omega_p$ and $\tau$ (red), lmf with $\omega_p$ only (blue), and by QSGW (green)

fplot -x .25,10  -y 0,8 -ord x3 -lt 1,col=1,0,0 eps.qsgw -ord x3 -lt 1,col=0,0,1 eps.qsgw0 -ord x3 -lt 1,col=0,1,0 dat-24

You should see a figure very similar to the one shown on the left. It is also useful to compare $\mathrm{Re}\,\epsilon(\omega)$:

fplot -x .25,10  -y -150,10 -ord x2 -lt 1,col=1,0,0 eps.qsgw -ord x2 -lt 1,col=0,0,1 eps.qsgw0 -ord x2 -lt 1,col=0,1,0 dat-24

You should see a figure very similar to the figure on the right.

Both $\mathrm{Re}\,\epsilon(\omega)$ and $\mathrm{Im}\,\epsilon(\omega)$ are similar between all three calculations, but there are differences, which can be summarized as follows:

• lmf uses the momentum operator for matrix elements, scaled by the correction Eq. [2] in the theory section, while the QSGW optics evalues matrix elements of $\exp(i \mathbf{q}{\cdot}\mathbf{r})$, with a small q.
• lmf does not pick up the intra-atomic term, but it can be added by a separate determination of $\omega_p$, via kkt.
• GW optics picks up both intra and interatomic terms, but without scattering and valid only above a cutoff frequency.
• The intra-atomic term calculated by two different ways are nearly identical.
• GW optics includes nonlocal field effects, while lmf does not.
• The scattering has a marked effect on $\mathrm{Im}\,\epsilon(\omega)$, but differences are much less noticable for $\mathrm{Re}\,\epsilon(\omega)$. This is because $\mathrm{Re}\,\epsilon$ is independent of $\tau$ once $\omega\tau{\gg}1$, while $\mathrm{Im}\,\epsilon$ is proportional to $1/\tau$.
• lmf neglects local fields.

References

[1] P. B. Johnson and R. W. Christy. Optical Constants of the Noble Metals. Phys. Rev. B 6, 4370 (1972)

[2] E.D. Palik and G. Ghosh. Handbook of optical constants of solids. Academic Press, 1998.

[3] Honghua U. Yang, Jeffrey D’Archangel, Michael L. Sundheimer, Eric Tucker, Glenn D. Boreman, and Markus B. Raschke, Optical dielectric function of silver. Phys. Rev. B 91, 235137 (2015)

[4] R. Del Sole and Raffaello Girlanda, Optical properties of semiconductors within the independent-quasiparticle approximation, Phys. Rev. B 48, 11789 (1993).

[5] Bernardo S. Mendoza and W. Luis Mochán, Ab initio theory of the Drude plasma frequency. J. Opt. Soc. Am. B 38, 1918-1926 (2021).

[6] Michael Marder, Condensed Matter Physics. John Wiley & Sons. DOI 10.1002/9780470949955. See Eqn 23.25.