Phonon-assisted luminescence by finite atomic displacements

Finite atomic displacements

In recent years, several articles have been published with Yambo, who calculated the phonon-assisted luminescence by means of finite atomic displacements. In this tutorial we will show how to use YamboPy in combination with Yambo to perform these kind of calculations. For a general review of finite displacements calculations of vibrational properties we recommend this article.
Notice that calculation of phonon-assisted luminescence by finite atomic displacements are quite involved and expensive from a computational point of view, for this reason they are limited to relatively simple systems.
In this tutorial we will follow the approach of Ref., and calculate luminescence using the Roosbroeck–Shockley (RS) relation applied to the excitonic case.
If you use this approch to calculate luminescence or phonon-assisted absorption please cite references and .
Important: This tutorial requires that you are familiar with Bethe-Salpeter calculations with Yambo, and phonons with Quantum-Espresso, and YamboPy. If it is not the case please study the corresponding tutorial before continuing.

Locate the lowest indirect exciton (or indirect band gap)

In this tutorial as example we will consider the hexagonal-BN, that is an indirect insulator. The first step to calculate the phonon-assisted luminescence is to locate the lowest indirect exciton. The cleanest strategy to identify the lowest excitation would be to calculate the excitons along the whole Irreducible Brillouin zone(IBZ), then interpolate them as explained in Interpolate exciton dispersion and find the q point that correspond to the minimum.
A simpler approach is to use the GW (or Kohn-Sham) band structure to find the momentum responsible for the indirect band gap, in the major part of cases this correspond to the lowest exciton.
In the h-BN case we have:

The indirect gap is between a point close to K and the M point, see panel (a) of the figure above. We will approximate the momentum responsible for the indirect emission with q=K-M. Notice that then we should build supercells that contain this moment, so it is good to approximate it with fractions of integers that are not too large so as not to have giant supercells. In this case we get q=(1/3,-1/6, 0).
Notice that if you have a direct band gap material, the momentum q corresponding to the lowest transition is q=0 and you do not need to construct special supercells to sample the corresponding vibration, it will be sufficient to displace atoms in the primitive cell.

Calculate phonons at momentum q of the indirect transition

Now we need to calculate all atomic vibrations compatible with the momentum of the lowest indirect transition. In order to do so we first calculate the phonons with QE, here the inputs hBN_phonons.tgz, then we calculate the force constants with q3r.x and now we can interpolate phonons along all the BZ as shown in panel (b) of Figure 1. Then we use matdyn.x to calculate phonons at momentum q=K-M = (1/3, -1/6, 0):

&input
   asr='crystal',
   flfrc='bn.fc',
   flfrq='bn.freq',
   deltaE=1.d0,
   q_in_cryst_coord=.true.
   fleig='bn.eig'
/
2
0 0 0
0.3333333333333333 -0.16666666666666666666 0.0

The obtained phonon modes correspond phonons at the vertical red line in panel(b) of Fig. 1. They are written in the matdyn.modes file with the corresponding eigenvectors and will be later used to displace atoms along the phonon modes. If you open the 'matdyn.modes' file you will find phonon modes for the q=0 and q=(1/3, -1/6, 0) that in cartesian coordinates is (1/3,0,0). We copy the modes corresponding to the second q=(1/3, -1/6, 0) in a separate file, called Q2.modes that will look like:

     diagonalizing the dynamical matrix ...

q =       0.3333      0.0000      0.0000
**************************************************************************
    freq (    1) =       5.194139 [THz] =     173.257815 [cm-1]
(  0.000000  -0.000000    -0.000000  -0.000000    -0.274372  -0.475224   )
( -0.000000   0.000000    -0.000000  -0.000000    -0.445962   0.000001   )
( -0.000000  -0.000000    -0.000000  -0.000000    -0.222981  -0.386214   )
(  0.000000   0.000000    -0.000000   0.000000    -0.548742   0.000000   )
...........................

Generate the supercell

Now we use YamboPy to generate the supercell corresponding to the momentum q=(1/3, -1/6, 0). A tutorial is already present in the file tutorial/supercell/generate_supercells.py. You can edit the file to specify the momentum q, the name of the SCF initial file (uc_filnm) and the name of the phonon mode file (mode_file).

   .....
   uc_filnm = 'hBN.scf.in' #pw input
   Q = [[1,-1,0],[3,6,0]] #qpoint to fold in nondiagonal supercell: Q=(1/3,-1/6,0)
   kpoints = [12,12,4]    #NB: these fractional crystal coordinates must exactly divide the kpoint mesh!
   modes_file = 'Q2.modes' #file with phonon eigenmodes
   .....

run the script to generate the supercell without displace atoms python generate_supercells.py -N. The script will generate a non-diagonal supercell in the file "sc_nondiagonal.bn.scf' that is shown below:

The YamboPy script automatically reduce the number of k-point in such a way to have the same sampling of the corresponding primitive cell, however we advice you to check if the new k-point grid is what you expect, and modify it in case of incorrect grid.
In the same way you can generate the input file for non-self-consistent calculation in the supercell.

Calculate the optical spectra in the supercell

Now you have to calculate the optical spectra and exciton in the new supercell. Do not forget to increase the number of conduction bands in both the screening the Bethe-Salpeter equation. For example in the supercell shown above there are 24 atoms, 6 times more than in the original cell, so you will have to multiply the number of bands by 6, if in the base cell you have used 40 bands here you will have to put 240 to have the same convergence parameters. Then print exciton according to their energy with the command ypp -e s.

N. B. : in order to not recalculate W for each atomic displacement you can perform calculations without symmetries in the pristine supercell and then copy the corresponding dielectric constant in the calculations with displaced atoms, this is usually a safe approximation described in Ref..

After solving the BSE in the supecell you will find the excitons that were already present in the primary cell (within the numerical noise) and new lower energy excitons, with zero dipole matrix elements, that correspond to the indirect excitons that have now been mapped at q=0, something like:

#    Maximum Residual Value =  0.38659E+01
#   
#    E [ev]             Strength           Index
#
    5.7165937         0.47290171E-11      1.0000000         <-- Indirect exciton mapped at q=0
    5.7166047         0.58239158E-09      2.0000000         <-- Indirect exciton mapped at q=0
    5.7367945         0.31830746E-09      3.0000000         <-- Indirect exciton mapped at q=0
    5.7368207         0.33674005E-10      4.0000000         <-- Indirect exciton mapped at q=0
    5.80542803        0.100763231E-6      5.00000000        <-- Direct exciton
    5.80598068        0.384604526E-8      6.00000000        <-- Direct exciton
    5.88944197         1.00000000         7.00000000        <-- Direct exciton
    5.88967037        0.102969212E-1      8.00000000        <-- Direct exciton
    ................

This result shows that indirect exciton without the coupling with phonon do not contribute the optical response.

Displace atoms and calculate dipoles derivatives and luminescence

Now you can generate the displaced supercell for all the phonon modes compatible with the indirect band gap q point by doing: python generate_supercells.py -S. The script will produce a set of input files for each of the 12 phonon modes at q=(1/3,-1/6,0). You have to perform BSE calculation for each one of these modes and get exciton dipoles.
For each phonon mode the atoms are displace along the corresponding eigenvector read from the "Q2.modes" file multiplied for 0.1 Borh (the variable Temp=0.1), see the following line in "generate_supercell.py":

 ...  
 sc.displace(eivs,nd_atom_positions,Temp=0.1) # Generate list of displaced supercells as PwIn objects called 'modes_qe'
 ...

this is usually a good value to get correct derivatives.
Once you get dipoles for a give phonon mode you can calculate the second order dipole derivatives using a finite different formula and the fact the first order derivatives are zero by construction.
Now you have all the information necessary to construct the emission spectra.
Using the formula, derived in Ref. you obtain the phonon-assisted luminescence as:

in this formula:

• n_r(ω) is the refractive index, that usually can be approximated with a constant in the luminescence energy range.
• E^S are the excitonic energies in the undistorted supercell, including indirect excitons
• B(E^S,T_exc) is the Bose occupation of excitons, that can be approximated with a Boltzman function, see Ref.~
• T_exc is the excitonic temperature that should be obtained from experimental measurements or used as a parameters, for a discussion see Ref.
• d²|T^S|²/dR²_λq are the derivatives of the exciton dipole matrix elements for a given phonon mode
• Ω_λ,q are the phonon frequencies and the sum on q is on the different supercells obtained from phonon modes responsible for the luminescence, in this tutorial there is just a single q point.
• η is the exciton line-width that is given as input parameter to define the broadening of the spectra

In general due to the Bose (or Boltzman) factor only the lower excitons (often only the lowest one) contribute to the luminescence signal.
Notice that in this formulation we occupied exciton as Boson (or Boltzman particles), in principle it is possible to define a excitonic occupation starting from a quasi-equilibrium distribution of excited electron-hole, for a discussion see Ref..
In order to speed up calculations you can perform all calculation without symmetries and copy the dielectric constant (ndb.em1s* or ndb.pp* files) from the undistorted supercell to the displaced ones. In this case if you do not need to calculate GW correction, for example using a scissor, you can also reduce the number of bands in the the non-self-consistent calculation to the one required by the BSE.
If you apply above formula considering only the lowest two excitons plus an appropriate scissor operator in the hBN case you get:

Notice that the dipole matrix elements that will enter in optical absorption and luminescence depend form the direction of the electric field. In this example the luminescence has been averaged on x and y directions because we were interested in the in-place luminescence. In general is it possible to select specific directions to reproduce particular experimental setups.