[Tutorial pdf document]
In this tutorial you will learn how to carry out basic simulations of RAS and HREELS spectra at semiconducting surfaces using the Yambo code. In particular we will demonstrate how reasonable agreement with experimental data can be achieved even at the independent particle (RPA) level of approximation, and illustrate the effect of higher order corrections on the computed spectra. This feature is available only in the yambo_surf executable, so make sure you have compiled it!
The tutorial is split into two parts. In the first part we will deal with RAS (reflectance anisotropy spectroscopy) calculations on the GaAs(001)c(4x4) surface using an asymmetric slab. In the second part we will consider a symmetric Si(001)c(4x2) surface, demonstrating the calculation of RAS including many body effects, and the calculation of HREELS (highresolution electron energy loss spectroscopy) data.
Once the tutorial zip file is unzipped the following folder structure will appear
COPYING README Surface_spectroscopy/
with the latter folder containing a directory for each structural model
>ls Surface_spectroscopy
GaAs001c4x4/ Si001c4x2/
Each system in turn is composed of the subdirectories
PWSCF/ YAMBO/
The reflectance anisotropy signal is defined as the difference between normalized reflectivities for light polarized along perpendicular surface directions x and y. Since the reflectivity from the bulk (in cubic materials) is isotropic, this difference spectra is completely of surface origin. Equivalent expressions for the RAS were derived by Bagchi, Barrera and Rajagopal [5] and by McIntyre and Aspnes [3]. In the language of a threelayer model, the RAS is computed from:
Hence, the main quantity that is calculated in this runlevel is the dielectric tensor of the surface layer. We will mostly work within the independent particle scheme, in the long wavelength limit, and without local fields (i.e. equivalent to yambo o c, with NGBlk=1). The tensor will be calculated for the (perpendicular) polarizations specified by q0x and q0y. Later we will show how to add some higher order effects.
GaAs(001)c(4x4)


Let's begin with the GaAs(001)c(4x4) surface. Enter the YAMBO folder
>cd GaAs001c4x4/YAMBO/
>ls
Experiment INPUTS REFERENCE SAVE
and run yambo_surf in the usual way to generate the setup databases:
yambo_surf F INPUTS/01_init
Note that we use a fairly low number of Gvectors throughout just to keep things running quickly.
The RAS runlevel can be launched using yambo_surf s r.
Fairly complete descriptions of the input file variables are available on the RAS runlevel page.
The preprepared input file INPUTS/02_RAS indicates the runlevel:
sursp # [R] Surface Spectroscopy
ras # [R] Reflectance anisotropy spectroscopy (RAS)
Following some standard variables that you should be familiar with by now are variables specific to the RAS calculation.
% q0x
1.000000  0.000000  0.000000  # [RAS] [cc] Xpolarization direction
%
% q0y
0.000000  1.000000  0.000000  # [RAS] [cc] Ypolarization direction
%
These define the two surface polarization vectors in cartesian coordinates. They must be orthogonal, and form part of a righthanded frame with the surface normal.
There are also several lines allowing the
BulkFile= "Experiment/EPSILON_GaAs_bulk.dat" # [RAS] File containing bulk
BulkForm= "2KK" # [RAS] Format of bulk (`3ri`,`3ir`,`2KK`,`3KK`)
BlkBroad= 0.10000 eV # [RAS] Broadening of bulk eps
Have a quick look first at the contents of the bulk file, defined by BulkFile, and plot the first column with gnuplot.
The form of the data in the bulk file is determined by BulkForm.
Here we are using the experimental bulk epsilon of GaAs for no particular reason  one should really use a computed one as in the case of Si(001) below. As only the imaginary part is given,
the real part is generated through a KramersKronig transformation, using a broadening BlkBroad.
Go ahead and launch the code using the predefined input file:
yambo_surf F INPUTS/02_RAS_full
...
<> [04] Surface spectroscopy module
<> [04.01] Surface geometry and spectral analysis
<> [04.02] Bulk epsilon data import
<> Processed 101 data [ 0.0010.00eV] from 42 records [ 2.03 5.62eV]
<> [04.03] Calculating eps cell.
<> Absorption @ q  x : 0.70711 0.70711 0.00000
<> [XCG] R(p) Tot o/o(of R) : 18592 38704 100
<01s> Xo@q[1] 1101 #################### [100%] (E) (X)
<01s> Absorption @ q  y : 0.70711 0.70711 0.00000
<02s> Xo@q[1] 1101 #################### [100%] (E) (X)
<02s> [04.04] Write eps cell for plotting
<02s> [04.05] Slab symmetries
<02s> [04.06] Surface dielectric function
<02s> [04.07] RAS
<02s> [05] Game Over & Game summary
If you set up the bulk file incorrectly, the code will issue a warning and continue, thus computing anyway the SAVE/ndb.dipoles database.
It's best not to use the J option here, as the ndb.dipoles will be used again.
You should now have some new output files o.ras*, o.epscell*, o.epssurf*.
Check first inside the report file r_sursp_ras* that the geometry is as expected,
[04.01] Surface geometry and spectral analysis
==============================================
Polarization directions:
 X 1.000000 0.000000 0.000000
 Y 0.000000 1.000000 0.000000
Surface normal : 0.000000 0.000000 1.000000
Surface vector indices (ix,iy,iz) : 1 2 3
Supercell height az (a.u.) : 42.40000
No cut off function in matrix elements.
and now plot the contents of the RAS output file o.rasrpa.
Not so useful. Look at the diagram of the surface structure model. Do you see why?
Now change the polarization directions to correspond to meaningful directions in the surface plane: put q0x perpendicular to the dimer axes, q0z parallel to the dimer axes (the X and Y in the figure are the Cartesian axes, as used in the DFT input file) and run the code again.
The new output should look something like this:
[04.01] Surface geometry and spectral analysis
==============================================
Polarization directions:
 X 0.707107 0.707107 0.000000
 Y0.707107 0.707107 0.000000
and the output file o.rasrpa_01:
# Version 3.3.1 Revision 2087
# http://www.yambocode.org
#
# 1 2 3 4 5 6 7 8
#hw (eV) RAS e1(xy) e2(xy) A B A x e2 B x e1
0.000000 0.000000 1.341996 0.000000 0.123842 0.000000 0.000000 0.000000
0.1000E1 0.5830E7 1.342 0.2266E2 0.1238 0.1810E4 0.2807E3 0.2429E4
0.2000E1 0.2333E6 1.342 0.4534E2 0.1238 0.3613E4 0.5615E3 0.4850E4
The different columns help in the interpretation of the computed RAS spectra, by splitting the RAS into surface and bulk components
(note the bulk dielectric function in the denominator).
Therefore, it is useful to rewrite the RAS as
You should now be able to identify which peaks in the computed RAS derive from the absorption, and thus can be easily analysed.
The "zero" position and extent of the cutoff function are determined by the variables CutZero and CutStep .
To use this, edit the relevant lines in the input, e.g. INPUTS/04_RAS_front (uncomment Cutoff):
%
Cutoff # [RAS] Cutoff mode (Uncomment to use)
CutZero= 0.250000 # [RAS] Zero position of cutoff fn (Frac)
CutStep= 0.500000 # [RAS] Width of cutoff function (Frac)
Note that the values here are in fractional units of the cell axis, i.e. [0..1] range.
Run the code now with
yambo_ras F INPUTS/04_RAS_front J 04_RAS_front
and note the creation of a new database in 04_RAS_front/ndb.dipoles_cut. This run is relatively slow because of a double sum on G.
This will remove the contribution from the back surface. Let's go one more step and isolate the contribution from the topmost layers containing the dimers:
yambo_ras F INPUTS/05_RAS_dimer J 05_RAS_dimer
and have a look at the results all together.
We can understand three important things:
Last, let's see how our calculated spectrum compared with experiment. GaAs(001)c(4x4) can be prepared in different ways. During MBE (molecular beam epitaxy) growth, a different reconstruction appears depending on whether As2 molecules (beta phase) or As4 molecules (alpha phase) are used. The measured RAS data appears in the Experiment folder: make a comparison with our "best" calculated spectrum. Can you now identify which phase corresponds to the GaAs heterodimer structure we have assumed?
You're right  its the alpha phase.
The beta phase in fact has pure unbuckled AsAs dimers. Hence the surface richer in As (beta) is grown using As2 molecules, opposite to what you might expect!
The agreement is already very good, in spite of the approximations we have used in the calculation.
What do you think is the most important parameter to converge/change in order to improve
the agreement with experiment across the whole spectral range?
(2) Si(001)c(4x2)


Let's move onto the Si(001)c(4x2) surface. Enter the YAMBO folder
>cd ../Si001c4x2/YAMBO/
>ls
Experiment INPUTS REFERENCE SAVE
and again run yambo_surf to generate the setup databases:
yambo_surf F INPUTS/01_init
Here we have a centrosymmetric slab, and don't necessarily need to use the cutoff function (do you see why a symmetric slab is not suitable for studying RAS of GaAs(001)?)
You can quickly confirm this using:
yambo_surf F INPUTS/02_RAS_full J 02_RAS_full
yambo_surf F INPUTS/03_RAS_half J 03_RAS_half
Note that here we are using a properly calculated Si bulk epsilon datafile (DFTRPA), and have rigidly shifted it also by 0.5eV.
BlkShift= 0.500000 eV # [RAS] GW shift of bulk eps
BulkFile= "Experiment/EPSILON_Si_bulk_RPA.dat" # [RAS] File containing bulk
Plot the two output files, and compare with experimental data in the Experiment folder.
Again, it already looks pretty good at the RPA level. This ab initio modeling is a bit of a joke!
A centrosymmetric slab allows us to move up a level, and add in some many body effects (at least until a real space cutoff is included...hi, Bernardo!).
We do this through separate yambo calculations for the same two polarizations as before, and post process the results.
Let's assume you've already done a GW calculation for the slab, which is no mean feat, and obtained an opening of the gap of +0.6eV.
From this we can do a rough BSE calculation, just looking at the low energy part of the spectrum to save time.
Remember the steps from previous lessons: first calculate the screening yambo b, then construct and solve BSE kernel yambo o b k sex y h.
yambo F INPUTS/03_W
yambo F INPUTS/04_BSE_X
yambo F INPUTS/04_BSE_Y
Have a look inside these input files: you should be familiar with most of the variables by now. Compare the components of epsilon for the cell,
computed within BSE and within RPA, from our earlier calculation.
Finally a trivial postprocessing step ypp_surf a ras, where the two output files are read and the RAS computed in the normal way:
ypp_surf F INPUTS/05_BSE_RAS
Just take care to define the parameter d correctly: it should be half the cell height in this case.
The result? Not a huge amount of difference between the RPA and the BSE results...It depends mostly on the value of the scissors operator correction. Can you understand why this is, by looking at the formula for the RAS?
This expression consists of two parts: the kinematic term, which contains all the information about the experimental geometry, and the loss function, which contains all the interesting juicy physics of electron scattering. Within the dipole regime, we can use a true threelayer model to derive the loss function, similar to what was done for RAS. The most flexible and precise version of this model is that of Del Sole/Selloni[6]:
where
So you see, also in the case of HREELS, the problem reduces to one of computing the dielectric function of a thin surface layer. The dependence on transferred momentum q is partially included, as we take the optical dielectric functions for convenience. This works quite well for typical experiments. All the tricks we discussed before regarding the cutoff function and so on apply also here, but be careful: the d parameter is actually used here, and must be determined correctly.
The modeling is a little complicated due to the kinematic term and different possible scattering geometries. Here we just take you through the basic steps of the calculation.
A typical input file is generated using yambo_surf s e (extra options available with V qp or V gen):
yambo_surf s e F INPUTS/06_REELS_full
Go ahead and run the code
yambo_surf F INPUTS/06_REELS_full J 06_REELS_full
and plot the output. You should be able to guess what the peaks at 1.8, 3.5, 5.1eV correspond to! (Hint: remember the RAS spectra).
One major approximation being made here is that the electron does not penetrate the crystal surface. We can estimate its penetration depth from the universal scattering curve:
At 40eV incident energy, it is about 6A, or 12 au. We can use this information then to augment the surface loss function with a term that represents excitation of bulk modes by the electron inside the material. This is done by defining a nonzero penetration depth:
yambo_surf s e F INPUTS/06_REELS_bulk J 06_REELS_bulk
Run the code again, with this modified input file,
yambo_surf F INPUTS/06_REELS_bulk J 06_REELS_bulk
and now compare with the experimental data in the Experiment folder.
It should be clear now that the higher energy features are excitations from the E1 and E2 critical points of the bulk.