Exciton-phonon coupling and luminescence - Yambo postprocessing

Step 7: Obtain the exciton-phonon coupling

Now, we can finally access our basic building block for exciton-phonon physics. This could be done entirely in python (using Yambopy), or by running Yambo. In this version of the tutorial we present the latter case.

Our objective is obtaining the following quantity:

[math]\displaystyle{ \mathcal{G}^\mu_{\alpha\lambda}(0,q)=\sum_{vv^\prime c k} A^{\alpha, *}_{cv^\prime} (k, q) g_{vv^\prime}^\mu (k,q) A^{\lambda}_{cv}(k,q) - \sum_{cc^\prime vk} A^{\alpha, *}_{c^\prime v} (k+q, q) g_{c^\prime c}^\mu (k+q,q) A^{\lambda}_{cv}(k,q) }[/math]

Here, [math]\displaystyle{ A^{\lambda}_{cv}(k,q) }[/math] are the exciton coefficients extracted from the eigenvector of the two-particles Hamiltonian during the BSE calculation in step 5, while [math]\displaystyle{ g_{nm}^\mu (k,q) }[/math] are the electron-phonon coupling matrix elements obtained in step 6. As you can see, the exciton [math]\displaystyle{ \lambda }[/math] undergoes phonon-mediated scattering to state [math]\displaystyle{ \alpha }[/math] via phonon mode [math]\displaystyle{ \mu }[/math]. The scattering can happen for the hole (valence, [math]\displaystyle{ v }[/math]) or for the electron (conduction, [math]\displaystyle{ c }[/math]).

NB:

(1) This is written in the "backward" momentum transfer convention used by Yambo. The momentum dependence is different in the "forward" transfer convention.

(2) For simplicity, this is written for zero initial exciton momentum. This means that one of the two states involved in the phonon-mediated scattering process will be in the optical limit (and possibly an optically generated exciton), while the other state can have any momentum: this momentum will be the same as the phonon one. This matrix element can be used to describe phonon-assisted absorption and emission spectra.

We have to write a new yambo input, that we can call excph.in, for this. You can copy (and adapt) the one below, or you can generate one by running from the command line:

yambo_ph -excph o

This generates an input to compute luminescence ("o" is for "optics"). The variables that we are interested in are:

excph                            # [R] Exction-phonon
ExcGkkp                          # [R][EXCPH] Exciton-Phonon Matrix Elelements
% ELPhExcStates
 1 | 4 |                             # [EXCPH] Incoming (external) exciton states
%
% ELPhExcSum
 1 | 12 |                             # [EXCPH] Outgoing (virtual) exciton states
%
LoutPath= "./bse_Lfull"                 # [EXCPH] Path of the outgoing L
% ElPhModes
  1 | 9 |                           # [ELPH] Phonon modes included
%

In this input, we have to select the initial exciton states [math]\displaystyle{ \lambda }[/math] with ELPhExcStates, the final exciton states [math]\displaystyle{ \alpha }[/math] with ELPhExcSum and the phonon modes [math]\displaystyle{ \mu }[/math] with ElPhModes. Here we consider the first four states at [math]\displaystyle{ Q=0 }[/math] (corresponding to just two excitons because they are both doubly degenerate -- do not break degeneracies when selecting states!) and the first twelve states at each finite-[math]\displaystyle{ q }[/math] point. We also include all the nine phonon modes of monolayer MoS2.

What about LoutPath? This variable controls the directory where the databases for the final-state excitons [math]\displaystyle{ \alpha }[/math] is located, which can be different from the directory with the initial-state excitons [math]\displaystyle{ \lambda }[/math] read as usual with the -J option when running the code. This makes it possible to compute the coupling between different exciton kinds. However, for our tutorial, we stick with the previously computed Lfull in both cases.

When we are satisfied with the input, we run the code using yambo_ph:

yambo_ph -F excph.in -J excph,bse_Lfull -C excph

If you check the output, you should find the ndb.excph* databases in the excph directory.

Analysis of the couplings

It is a good idea to have a look at what we computed up to now in order to make sure nothing has gone wrong.

It is not easy to know what to expect (apart from symmetry and gauge compliance of the matrix elements), but one can work out the exciton-phonon selection rules in advance, check that the magnitude is reasonable, etc.

It is also not easy to meaningfully plot this quantity. We have to make sure that we are not breaking degenerate states, otherwise the plots will not be invariant.

First of all, we have to know our system: in monolayer MoS₂, the first four excitons are all doubly degenerate. The first exciton responsible for a bright peak in the absorption spectrum (the A peak), is the second state, corresponding to state indices (3,4) in fortran indexing or (2,3) in python indexing.

All these information can be obtained by analyzing the BSE results (this stuff is explained in the BSE tutorials) and by knowledge of the system or class of systems from the literature.

Thus, a good quantity to plot may be the norm of the matrix elements, summed over the degenerate subspace of exciton A, for a certain number of scattered final states mediated by certain phonon modes:

[math]\displaystyle{ F_A(q)= \sqrt{ \sum_{\alpha \in A,\lambda,\mu} |\mathcal{G}_{\alpha\lambda}^\mu (0,q)|^2 } }[/math]

In order to do this, we create a python script analyse_excph.py in which we first load the excph dabatases using the YamboExcitonPhononDB from yambopy. You can find this script in the yambopy directory, in tutorials/exciton-phonon. First, we select the exciton and phonon states to be included in F_A, together with the path of databases and plot details:

# Exciton in states
exc_in  = [2,3]     # A: 2,3 -- B: 6,7
exc_out = [0,1,2,3] # first 4 states (dispersion of triplet state and A)
ph_in  = 'all'

# Paths of databases
ns_db1 =f'{path}/SAVE/ns.db1'
ndb_exc=f'{path}/excph'

...

Then, we load the data:

# Read lattice and k-space info
ylat = YamboLatticeDB.from_db_file(filename=ns_db1),Expand=True)
print(ylat)

# Read exc-ph databases
X = YamboExcitonPhononDB(ylat,save_excph=ndb_exc)
print(X)

The quantity [math]\displaystyle{ F_A(q) }[/math] is defined inside the plotting function as

G_squared = excph.excph_sq
G2plt = np.zeros(len(G_squared))

if exc_in  == 'all': exc_in  = range(G_squared.shape[2])
if exc_out == 'all': exc_out = range(G_squared.shape[3])
if ph_in   == 'all': ph_in   = range(G_squared.shape[1])

G_squared = G_squared[:, ph_in, :, :].sum(axis=(1))
G_squared = G_squared[:, exc_in, :].sum(axis=(1))
G_squared = G_squared[:, exc_out].sum(axis=(1))

F_q = np.sqrt( G_squared )*ha2ev # Switch from Ha to eV

And finally, we have to make a plotting function. For this tutorial we will use the default scatterplot provided by YamboExcitonPhononDB

excph.plot_excph(F_q,plt_cbar=plt_cbar,**kwargs)
...

You can get more experience on using Yambopy for these kinds of visualization by following the Yambopy tutorials. In fact, remember that this scripts and all the other yambopy tutorial scripts are just suggestions, not source code written in stone: if you know numpy and matplotlib you can do your own analysis and your own plots, you just need to import the required Yambopy modules to load the data.

In our case, the resulting plot is the following.

This can be checked against Fig. 2(d) of reference , although you have to keep in mind that our results are badly undersampled in terms of the reciprocal-space grid, as can be easily seen, and the quantity plotted is not exactly the same. However, the main features are already there since they are dictated mostly by crystal symmetries.

Now that we have the exciton-phonon matrix elements, we can use them to build several kinds of observables. Below, we give an example related to phonon-assisted luminescence, but we may update this tutorial in the future to include more cases.

Step 8: Compute phonon-assisted luminescence

We want to compute the experimental optical signature due to the phonon-assisted recombination of an exciton (as sketched in the figure).

The signal from the phonon replicas can be modeled as a second-order scattering process involving one phonon and one photon. In addition, we can calculate the renormalization of the "direct" exciton recombination, due to the appearance of the satellites. However, in order to model both phonon replicas (satellites) and main peak renormalization, we pay the price of approximating the replicas oscillator strength. In the yambopy postprocessing, currently we do not compute peak renormalization but rather we calculate the full expression for the satellites.

The luminescence expression implemented in Yambo comes from the first-order exciton-phonon response function adapted using the Van Roosbroeck-Shockley relation. It employs the diagonal exciton-phonon self-energy only, so that we can call this "exciton-quasiparticle approximation".

[math]\displaystyle{ I^{PL}(\omega)=I^{(0)}(\omega)+I^{Sat}(\omega) }[/math]

The first-order renormalised direct part is

[math]\displaystyle{ I^{(0)}(\omega)=-\omega \sum_\alpha \omega^2 \frac{|D_\alpha|^2 [1-R_\alpha^{em}(T)-R^{abs}_\alpha (T)]N_\alpha^{exc}(T_{exc})}{\omega - E_\alpha +\mathrm{i}\eta} }[/math]

The satellite parts are (phonon absorption | photon emission):

[math]\displaystyle{ I^{Sat,abs}(\omega)=-\omega\sum_{\beta q \alpha \mu}\left\{ \frac{1}{N_q}\frac{(\omega-2\Omega_{\mu q})^2 |D_\alpha |^2 T^{em}_{\mu \beta \alpha q} N^{exc}_{\beta q}(T_{exc})n_{\mu q}(T)}{\omega - [E_{\beta q}+\Omega_{\mu q}]+\mathrm{i}\eta} \right\} }[/math]

and (phonon emission | photon emission):

[math]\displaystyle{ I^{Sat,em}(\omega)=-\omega\sum_{\beta q \alpha \mu}\left\{\frac{1}{N_q} \frac{(\omega+2\Omega_{\mu q})^2 |D_\alpha |^2 T^{abs}_{\mu \beta \alpha q} N^{exc}_{\beta q}(T_{exc})[1+n_{\mu q}(T)]}{\omega - [E_{\beta q}-\Omega_{\mu q}]+\mathrm{i}\eta} \right\} }[/math]

Here, the [math]\displaystyle{ T }[/math]-factors contain the satellite weights

[math]\displaystyle{ T^{em}_{\mu\beta\alpha q}=\frac{|\mathcal{G}_{\beta\alpha}^\mu(q)|^2}{|E_{\beta q}-E_{\alpha}+\Omega_{\mu q}|^2}, \quad T^{abs}_{\mu\beta\alpha q}=\frac{|\mathcal{G}_{\beta\alpha}^\mu(q)|^2}{|E_{\beta q}-E_{\alpha}-\Omega_{\mu q}|^2} }[/math]

Which, once summed, form the renormalization factors (equal to minus the energy-derivative of the exc-ph self-energy evaluated at the "bare" exciton energy, analogous to GW)

[math]\displaystyle{ R^{abs}_\alpha = \sum_{\beta \mu q} \frac{1}{N_q} T^{em}_{\mu\beta\alpha q}[1+n_{\mu q}(T)], \quad R^{em}_\alpha = \sum_{\beta \mu q} \frac{1}{N_q} T^{abs}_{\mu\beta\alpha q}n_{\mu q}(T) }[/math]

In these expressions, the oscillator strength of the peak is given by the exciton-phonon coupling matrix elements [math]\displaystyle{ \mathcal{G} }[/math] multiplied by the exciton dipoles [math]\displaystyle{ D }[/math] (they are called "residuals" in Yambo). Here [math]\displaystyle{ E_\lambda }[/math] and [math]\displaystyle{ E_{\alpha q} }[/math] are the energies of the optical and finite-momentum excitons, respectively, while [math]\displaystyle{ \Omega_{\mu q} }[/math] are the phonon energies.

The phonon temperature-dependent Bose-Einstein occupation function is [math]\displaystyle{ n_{\mu q}(T) }[/math]. We include the processes of phonon emission ([math]\displaystyle{ \propto n(T)+1 }[/math]), as well as of phonon absorption ([math]\displaystyle{ \propto n(T) }[/math]). Therefore, [math]\displaystyle{ I^{(1)}_{PL}(\omega;T) }[/math] describes light emission by recombining excitons mediated by either phonon absorption or emission.

The quantity [math]\displaystyle{ N_{\alpha q}(T_{exc}) }[/math] is the exciton occupation function. Luminescence is technically an out-of-equilibrium process, but we can assume that for very low density of excitations and in steady-state conditions, the exciton population can be approximately described by an equilibrium distribution evaluated at an effective temperature. Here, we use the Boltzmann distribution. Experimentally, [math]\displaystyle{ T_{exc} }[/math] tends to coincide with the lattice temperature [math]\displaystyle{ T }[/math] more or less above 100 K, while at very low temperature (< 10 K), [math]\displaystyle{ T_{exc} }[/math] may vary between 10-50 K. It goes without saying that this needs to carefully be checked in your realistic calculations.

Running the jobs

In order to study luminescence in a paradigmatic system, we switch to bulk hexagonal boron nitride and we repeat the workflow. As you can easily see, one can think about automatizing the execution of all these calculations via scripting or more advanced tools. However, in the case of very large simulations (memory-limited or disk-space limited) or for systems whose electronic and lattice properties are fragile with respect to tiny calculation details, one must be very careful and run many basic tests.

Fortunately, we are running a fast underconverged example. We use LDA pseudopotentials from the pseudo-dojo library and the following are the calculations steps.

1. Input hbn.scf:

&control
    calculation='scf',
    prefix='hBN',
    restart_mode='from_scratch'
    pseudo_dir = '$PSEUDO_DIR',
    outdir = './tmp'
    verbosity = 'high'
    wf_collect=.true.
/
&system
    ibrav = 4,
    celldm(1) = 4.703675849
    celldm(3) = 2.603711434
    nat= 4,
    ntyp= 2,
    force_symmorphic=.true.
    ecutwfc = 100,
/
&electrons
  diago_david_ndim = 2
  diago_full_acc=.true.
  diago_thr_init=5.0e-6
  mixing_mode = 'plain'
  mixing_beta = 0.7
  conv_thr =  1.0d-16
/
ATOMIC_SPECIES
 B 10.81100  B_LDA_dojo.UPF
 N 14.00674  N_LDA_dojo.UPF
ATOMIC_POSITIONS {crystal}
N             0.6666666670        0.3333333330        -0.250000000000
B             0.3333333330        0.6666666670        -0.250000000000
B             0.6666666670        0.3333333330        0.25000000000
N             0.3333333330        0.6666666670        0.25000000000
K_POINTS {automatic}
6 6 2 0 0 0

mpirun -np 4 pw.x -inp hbn.scf > scf.out

2. Input hbn.nscf:

&control
    calculation='nscf',
    prefix='hBN',
    restart_mode='from_scratch'
    pseudo_dir = '$PSEUDO_DIR'
    outdir = './'
    verbosity = 'high'
    wf_collect=.true.
/
&system
    ibrav = 4,
    celldm(1) = 4.703675849
    celldm(3) = 2.603711434
    nat= 4,
    ntyp= 2,
    force_symmorphic=.true.
    ecutwfc = 100,
	nbnd = 120
/
&electrons
  diago_david_ndim = 2
  diago_full_acc=.true.
  diago_thr_init=5.0e-6
  mixing_mode = 'plain'
  mixing_beta = 0.7
  conv_thr =  1.0d-16
/
ATOMIC_SPECIES
 B 10.81100  B_LDA_dojo.UPF
 N 14.00674  N_LDA_dojo.UPF
ATOMIC_POSITIONS {crystal}
N             0.6666666670        0.3333333330        -0.250000000000
B             0.3333333330        0.6666666670        -0.250000000000
B             0.6666666670        0.3333333330        0.25000000000
N             0.3333333330        0.6666666670        0.25000000000
K_POINTS {automatic}
6 6 2 0 0 0

mpirun -np 4 pw.x -inp hbn.nscf > nscf.out

3. Input hbn.dvscf:

hbn_dvscf
&inputph
  tr2_ph=1.0d-12,
  verbosity='high'
  prefix='hBN',
  fildvscf = 'hBN-dvscf',
  electron_phonon = 'dvscf',
  fildyn='hBN.dyn',
  epsil=.false.,
  ldisp=.true.,
  recover=.true.,
  nq1=6,
  nq2=6,
  nq3=2
/

nohup mpirun -np 8 pw.x -inp hbn.dvscf > dvscf.out &

4. Input bse.in (we include 2 valence and 2 conduction bands):

optics                           # [R] Linear Response optical properties
bss                              # [R] BSE solver
bse                              # [R][BSE] Bethe Salpeter Equation.
dipoles                          # [R] Oscillator strenghts (or dipoles)
em1s
DIP_CPU= "1 8 1"                      # [PARALLEL] CPUs for each role
DIP_ROLEs= "k c v"                    # [PARALLEL] CPUs roles (k,c,v)
DIP_Threads=0                    # [OPENMP/X] Number of threads for dipoles
X_and_IO_CPU= "1 1 1 8 1"                 # [PARALLEL] CPUs for each role
X_and_IO_ROLEs= "q g k c v"               # [PARALLEL] CPUs roles (q,g,k,c,v)
X_and_IO_nCPU_LinAlg_INV=-1      # [PARALLEL] CPUs for Linear Algebra (if -1 it is automatically set)
X_Threads=0                      # [OPENMP/X] Number of threads for response functions
BS_CPU= "8 1 1"                       # [PARALLEL] CPUs for each role
BS_ROLEs= "k eh t"                     # [PARALLEL] CPUs roles (k,eh,t)
BS_nCPU_LinAlg_INV=-1            # [PARALLEL] CPUs for Linear Algebra (if -1 it is automatically set)
BS_nCPU_LinAlg_DIAGO=-1          # [PARALLEL] CPUs for Linear Algebra (if -1 it is automatically set)
X_Threads=0                      # [OPENMP/X] Number of threads for response functions
K_Threads=0                      # [OPENMP/BSK] Number of threads for response functions
% QpntsRXs
   1 | 14 |                         # [Xs] Transferred momenta
%
% BndsRnXs
   1 | 120 |                         # [Xs] Polarization function bands
%
NGsBlkXs= 10                Ry    # [Xs] Response block size
% LongDrXs
 1.000000 | 1.000000 | 1.000000 |        # [Xs] [cc] Electric Field
%
BSEmod= "resonant"               # [BSE] resonant/retarded/coupling
BSKmod= "SEX"                    # [BSE] IP/Hartree/HF/ALDA/SEX/BSfxc
Lkind= "Lfull"                    # [BSE] Lbar (default) / full
BSSmod= "d"                      # [BSS] (h)aydock/(d)iagonalization/(s)lepc/(i)nversion/(t)ddft`
% DipBands
   1 | 120 |                         # [DIP] Bands range for dipoles
%
DipApproach= "G-space v"         # [DIP] [G-space v/R-space x/Covariant/Shifted grids]
DipComputed= "R V P"             # [DIP] [default R P V; extra P2 Spin Orb]
BSENGexx= 30000            Ry    # [BSK] Exchange components
#ALLGexx                       # [BSS] Force the use use all RL vectors for the exchange part
BSENGBlk=  9000            Ry    # [BSK] Screened interaction block size [if -1 uses all the G-vectors of W(q,G,Gp)]
% KfnQP_E
 1.25997 | 1.08816 | 1.12683 |        # [EXTQP BSK BSS] E parameters  (c/v) eV|adim|adim
%
% BSEQptR
 1 | 14 |                     # [BSK] Transferred momenta range
%
% BSEBands
   7 | 10 |                         # [BSK] Bands range
%
% BEnRange
  0.50000 | 8.00000 |         eV    # [BSS] Energy range
%
% BDmRange
 0.050000 | 0.050000 |         eV    # [BSS] Damping range
%
BEnSteps= 1000                    # [BSS] Energy steps
% BLongDir
 1.000000 | 1.000000 | 0.000000 |        # [BSS] [cc] Electric Field
%
WRbsWF                        # [BSS] Write to disk excitonic the WFs

nohup mpirun -np 8 yambo -F bse.in -J bse_Lfull -C bse_Lfull &

Importantly, since we want to describe the phonon-assisted recombination process of an *optical* exciton (i.e., emitting a transverse photon), this time we also run an additional calculation at `Q=0` omitting the nonanalytic long-range Coulomb exchange. Make a second input bse_Lbar.in with the following changes:

Lkind= "Lbar"                    # [BSE] Lbar (default) / full
% BSEQptR
 1 | 1 |                     # [BSK] Transferred momenta range
%

4b. So now we make a second BSE run in a different directory specified by -J. Here, we also pass to yambo the directory of the previous run as it includes the important screening databases ndb.em1s* that we do not want to recompute from scratch.

mpirun -np 8 yambo -F bse_Lbar.in -J bse_Lbar,bse_Lfull -C bse_Lbar

5. Now we run lelphc with yambopy to get the el-ph matrix elements:

yambopy l2y -ph path/of/dvscf/hbn.dvscf -b 7 10 -par 4 2 

6. And finally we generate the exciton-phonon input excph.in using yambo_ph -excph o. Now, we take a look at all the additional variables that we didn't check before, specifying also the details for luminescence the spectrum calculation.

excph                            # [R] Exction-phonon
ExcGkkp                          # [R][EXCPH] Exciton-Phonon Matrix Elelements
ExcPhOptics                      # [R][EXCPH] Exciton-Phonon Optics
BoseTemp=10.000000         Kn    # Bosonic Temperature 
EXCTemp= 10.000000          Kn    # [EXCPH] Excitonic Temperature (for luminescence spectra)
% ELPhExcStates
 1 | 4 |                             # [EXCPH] Incoming (external) exciton states
%
% ELPhExcSum
 1 | 12 |                             # [EXCPH] Outgoing (virtual) exciton states
%
% ElPhModes
 1 | 12 |                           # [ELPH] Phonon modes included
%
LoutPath= "./bse-L_full"                 # [EXCPH] Path of the outgoing L
FANdEtresh= 0.100000E-5    eV    # [ELPH] Energy treshold for Fan denominator
EXCPHdEtresh= 0.100000E-5  eV    # [ELPH] Energy treshold for exc-ph denominator
LDamping= 0.500000E-3      eV    # [EXCPH] Damping of exc-ph self-energy
# These are the plot parameters, same as in other parts of yambo 
% EnRngeXd
  4.00000 | 5.00000 |         eV    # [Xd] Energy range
%
% DmRngeXd
 0.00500000 | 0.00500000 |         eV    # [Xd] Damping range
%
ETStpsXd= 4000                    # [Xd] Total Energy steps

Here BoseTemp and EXCTemp are the lattice and excitonic temperature, respectively.

Now we run the calculation with yambo_ph. Here, we read as [math]\displaystyle{ Q=0 }[/math] excitons with -J the ones without the long-range Coulomb exchange:

yambo_ph -F excph.in -J excph,bse_Lbar -C excph

and we find our outputs in the excph directory.

NB: Step 6 could have been equivalently run in yambopy, limiting the use of the yambo code to just step 4. This latter option is more flexible, as it allows for a greater degree of control by the user. We are in the last stages of the development and it will be available soon.

Results

If we check the output directory from the step 6 calculation, we find the o-excph.pl_bse_ph_ass output files containing the luminescence spectra. We can plot them with gnuplot or any other tool:

gnuplot
> set xrange[4.2:5]
> p 'o-excph.pl_bse_ph_ass' u 1:2 w l lc rgb 'red' lw 3

Here, the signal corresponds to a finite-momentum exciton that recombines with the help of several different phonon modes, both optical and acoustic. Each phonon mode whose coupling with the exciton is allowed can generate a peak, and the energy shifts of these peaks with respect to the initial exciton energy correspond to the phonon energies. The direct recombination peak is completely suppressed. This result is underconverged, but the main features are all there. In the plot, we show a more converged example using a 12x12x4 grid (all the other parameters being equal). These plots can be compared with Fig. 4(a) of reference .