View source for Exciton-phonon coupling and luminescence

[[File:Tdgw-phonon-usc-01-1024x829.jpg| 400px | right]]

In this advanced tutorial, we will calculate exciton-phonon interactions from first principles by interfacing DFPT (for phonon calculations) and BSE (for exciton calculations).

The DFTP calculations are run with Quantum ESPRESSO, while the many-body GW-BSE calculations are run with Yambo. Finally, the exciton-phonon interaction will be obtained by combining and postprocessing the databases computed in the two previous runs. The great advantage of this workflow is that the calculations can be run in the irreducible Brillouin zones both for the electronic momenta ('''k''') and the transfer momenta ('''Q''', '''q''') of excitons and phonons, thus speeding up considerably the jobs while reducing the IO and memory load.

We will first compute the exciton-phonon coupling matrix elements: these are the building blocks needed to construct experimental observables such as phonon-assisted optical spectra (such as luminescence), Raman spectra and exciton lifetimes. We will do this in the case of monolayer MoS<sub>2</sub>, a 2D system with large spin-orbit interaction.

As an example of application, we will consider the case of phonon-assisted luminescence. We will do this in the case of bulk hBN, a layered indirect insulator with strong electron-phonon coupling.

'''Note:''' this tutorial will be updated when new exc-ph tools become available in Yambopy (including full-python postprocessing, Raman spectra, interpolated lifetimes, etc).

== Requirements ==

This is an advanced topic: we assume that you already know something about the theory<ref name="toyozawa2003" /><ref name="antonius2017" /><ref name="cudazzo2020" /><ref name="paleari2019_PhD" /><ref name="paleari2022" /><ref name="lechifflart2023_PhD" /> and applications<ref name="paleari2019" /><ref name="cannuccia2019" /><ref name="chen2020" /><ref name="zanfrognini2023" /><ref name="lechifflart2023" /><ref name="chan2023" /><ref name="marini2024" /><ref name="murali2025" /> of exciton-phonon physics. 

Also, we assume that you already know how to run both a basic '''Yambo''' GW-BSE calculation and a DFPT phonon calculation with '''Quantum ESPRESSO'''.

Besides the QE executables <code>pw.x</code> and <code>ph.x</code>, we also use the yambo phonon-specific executable <code>yambo_ph</code> and the python utility '''Yambopy'''. The auxiliary code '''LetzElPhC''' (executable <code>lelphc</code>) will be used to obtain the electron-phonon matrix elements by reading the same electronic wavefunctions used by Yambo (and stored in the <code>SAVE</code> directory), while also making full use of crystal symmetries. [https://gitlab.com/lumen-code/LetzElPhC LetzElPhC] will be run by Yambopy, but it must nonetheless be installed. Finally, the exciton-phonon properties can be computed either using <code>yambo_ph</code> or using Yambopy itself.

[[File:Workflow scheme.png|800px|center]]

== Step 0: Pseudopotentials, equilibrium structure and convergence ==

In a real calculation, it is important to ensure that both the pseudopotential and the lattice parameters that we are using are compatible and perform well for the electronic excited states and for the lattice vibrations simultaneously. Furthermore, you have to make sure that the wave function cutoff <code>ecutwfc</code> is converged with respect to the DFPT step and not just to the DFT one. This is in addition to the other customary convergence tests for DFT, DFPT, GW and BSE calculations.

This is often the most time-demanding step when starting on a new system.

For the sake of this tutorial, we assume that we have already done all these tests and we are starting the final workflow to get the exciton-phonon properties.

== Step 1: scf calculation ==

First of all, we run a standard scf calculation with <code>pw.x</code> for Yambo. We stick with non-symmorphic symmetries. At the end, we will have the QE <code>save</code> directory.

This is the input <code>mos2.scf</code>:

 &control
           wf_collect = .true.,
          calculation = "scf",
            verbosity = 'high',
           pseudo_dir = '$PSEUDO_DIR',
               prefix = "mos2",
               outdir = '.',
  /&end
  &system
              ecutwfc = 100.0,
          occupations = 'fixed',
                ibrav = 4,
            celldm(1) = 5.9000811881,
            celldm(3) = 6.7795677253,
                  nat = 3,
                 ntyp = 2,
             lspinorb = .true.
             noncolin = .true.
             assume_isolated = '2D'
         force_symmorphic = .true.
  /&end
  &electrons
     electron_maxstep = 200,
          mixing_beta = 0.7,
             conv_thr = 1.d-08,
  /&end
   ATOMIC_SPECIES
    Mo  95.940      Mo_ONCV_PBE_FR-1.0.upf
    S    32.065     S_ONCV_PBE_FR-1.1.upf
   ATOMIC_POSITIONS { crystal }
 Mo       0.333333333   0.666666667   0.000000000
 S        0.666666667   0.333333333   0.073413577
 S        0.666666667   0.333333333  -0.073413577
  K_POINTS { automatic }
 6 6 1 0 0 0

Here we are using full relativistic pseudopotentials from the SG-15 database.

We can run it on our machine (for example using 4 MPI tasks) as:

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

== Step 2: nscf calculation for Yambo ==

Copy the QE <code>save</code> directory from the scf calculation and run the nscf calculation for any number of empty states, with the correct <code>k</code>-grid we want to use in Yambo. Here we are using a badly underconverged grid of 6x6x1.

This reciprocal-space grid will also match the momentum transfer <code>q</code> grid on which excitons and phonons will be defined!

The electronic wavefunctions computed at this step and stored in the new nscf <code>save</code> directory will be used both by Yambo and by the electron-phonon code: this is important because using different sets of wavefunctions would lead to a phase mismatch issue in the exciton-phonon matrix elements.

The nscf input <code>mos2.nscf</code> is
 &control
           wf_collect = .true.,
          calculation = "nscf",
            verbosity = 'high',
           pseudo_dir = '$PSEUDO_DIR',
               prefix = "mos2",
               outdir = '.',
  /&end
  &system
              ecutwfc = 100.0,
          occupations = 'fixed',
                ibrav = 4,
            celldm(1) = 5.9000811881,
            celldm(3) = 6.7795677253,
                  nat = 3,
                 ntyp = 2,
             lspinorb = .true.
             noncolin = .true.
                nbnd  = 250
             assume_isolated = '2D'
         force_symmorphic = .true.
  /&end
  &electrons
     electron_maxstep = 200,
          mixing_beta = 0.7,
             conv_thr = 1.d-08,
  /&end
   ATOMIC_SPECIES
    Mo  95.940      Mo_ONCV_PBE_FR-1.0.upf
    S    32.065     S_ONCV_PBE_FR-1.1.upf
   ATOMIC_POSITIONS { crystal }
 Mo       0.333333333   0.666666667   0.000000000
 S        0.666666667   0.333333333   0.073413577
 S        0.666666667   0.333333333  -0.073413577
  K_POINTS { automatic }
 6 6 1 0 0 0


Again, we run the calculation

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

== Step 3: dvscf phonon calculation ==

Now we run the phonon calculation.

Copy the <code>save</code> directory from the '''scf''' calculation and run <code>ph.x</code> for a dvscf calculation with a standard <code>q</code>-grid matching the <code>k</code>-grid we wanna use in Yambo.

At the end, we will have the <code>_ph0</code> directory containing the variation of the self-consistent potential, <math>\Delta V_{SCF}(q)</math>, and the <code>*.dyn</code> files with the phonon energies and eigenvectors.

NB: one could further refine the phonon energies by enforcing the acoustic sum rule, including non-analytic long-range contributions, interpolating to finer grids... all of this can be done within Quantum ESPRESSO and will not be covered in this version of the tutorial.

The input is <code>mos2.dvscf</code>:
 
 mos2_dvscf
 &inputph
   tr2_ph=1.0d-12,
   verbosity='high'
   prefix='mos2',
   fildvscf = 'mos2-dvscf',
   electron_phonon = 'dvscf',
   fildyn='mos2.dyn',
   epsil=.false.,
   ldisp=.true.,
   recover=.true.,
   nq1=6,
   nq2=6,
   nq3=1
 /

And now we run as

 nohup mpirun -np 8 ph.x -inp mos2.dvscf > dvscf.out &

This time we use nohup and more processes because this calculation may take some time. It is a good idea to set <code>recover=.true.</code> as in a real calculation you will easily breach walltime, and in this way you can safely restart.

== Step 4: create Yambo <code>SAVE</code> directory ==

This is just the standard Yambo initialization: run 
 
 p2y 

and then 

 yambo 

in the '''nscf''' <code>save</code> folder and then move the newly generated <code>SAVE</code> directory to a convenient place.

== Step 5: run a BSE calculation ==

Now we switch from QE to Yambo. Here, we forgo the GW step for simplicity (we can use a scissor operator to open the band gap).

This calculation has a couple of differences with respect to a standard BSE calculation for optical absorption. We can look at the input file <code>bse.in</code>:

 <span style="color:blue"># Runlevels</span>
 optics                       # [R OPT] Optics
 rim_cut                      # [R RIM CUT] Coulomb potential
 bss                          # [R BSS] Bethe Salpeter Equation solver
 em1s                         # [R Xs] Static Inverse Dielectric Matrix
 bse                          # [R BSE] Bethe Salpeter Equation.
 bsk                          # [R BSK] Bethe Salpeter Equation kernel
 <span style="color:blue"># RIM and cutoff settings</span>
 RandQpts=1000000             # [RIM] Number of random q-points in the BZ
 RandGvec= 100            RL    # [RIM] Coulomb interaction RS components
 CUTGeo= "slab z"               # [CUT] Coulomb Cutoff geometry: box/cylinder/sphere X/Y/Z/XY..
 <span style="color:blue"># Static screening</span>
 X_and_IO_CPU= "1 1 4 2 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)
 Chimod= "hartree"            # [X] IP/Hartree/ALDA/LRC/BSfxc
 % BndsRnXs
    1 |  200 |                 # [Xs] Polarization function bands
 %
 % LongDrXs
  1.000000 | 1.000000 | 1.000000 |        # [Xs] [cc] Electric Field
 %
 NGsBlkXs= 8000            mRy    # [Xs] Response block size
 <span style="color:blue"># BSE</span>
 BS_CPU= "4.1.2"                   # [PARALLEL] CPUs for each role
 BS_ROLEs= "k.eh.t"                 # [PARALLEL] CPUs roles (k,eh,t)
 BS_nCPU_diago=4              # [PARALLEL] CPUs for matrix diagonalization
 BSEmod= "causal"             # [BSE] resonant/causal/coupling
 BSKmod= "SEX"                # [BSE] IP/Hartree/HF/ALDA/SEX/BSfxc
 BSSmod= "d"                  # [BSS] (h)aydock/(d)iagonalization/(s)lepc/(i)nversion/(t)ddft`
 BSENGexx=  40000      mRy    # [BSK] Exchange components
 ALLGexx                      # [BSS] Force the use use all RL vectors for the exchange part
 BSENGBlk=  8000       mRy    # [BSK] Screened interaction block size
 <span style="color:red">Lkind="full"                  #[BSE,X] bar(default)/full/tilde</span>
 % KfnQP_E
  1.000000 | 1.000000 | 1.000000 |        # [EXTQP BSK BSS] E parameters  (c/v) eV|adim|adim
 %
 % BEnRange
   0.00000 |  4.00000 | eV    # [BSS] Energy range
 %
 % BDmRange
   0.05000 |  0.05000 | eV    # [BSS] Damping range
 %
 BEnSteps= 2000               # [BSS] Energy steps
 % BLongDir
  1.000000 | 0.000000 | 0.000000 |        # [BSS] [cc] Electric Field
 %
 <span style="color:red">% BSEQptR</span>
 <span style="color:red"> 1 | 7 |                     # [BSK] Transferred momenta range</span>
 <span style="color:red">%</span>
 % BSEBands
    25 |  28 |                 # [BSK] Bands range
 %
 WRbsWF                      # [BSS] Write to disk excitonic the FWs
 
This file was generated using the command: <code> yambo -X s -o b -k sex -y d -r</code><br>

First of all, we compute the excitons for all the momenta in the irreducible Brillouin zone for our discrete grid via the <code>BSEQptR</code> variable. This will be a '''finite-momentum''' BSE calculation, analogous to the phonon one.

Second, we change the variable <code>Lkind</code> from <code>bar</code> to <code>full</code>. In Yambo, <code>Lkind="bar"</code>, which is the default for optical absorption, means that we are computing the excitonic response function without the long-range component of the exchange interaction. This cannot be used when computing the exciton momentum dependence, where the long-range exchange interaction can play a role, therefore we have to include it with <code>Lkind="full"</code>. This allows for the calculation of the excitonic longitudinal-transverse splitting (in 3D systems) as well.

We can now run the code:

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

At the end of the calculation, we have obtained the <code>ndb.BS_diago_Q*</code> databases inside the directory <code>bse_Lfull</code>. They contain information on the exciton energies and wavefunctions at each momentum. Do not forget to check the report and logs of your calculation in the same directory to make sure that the code is doing what you want.

== Step 6: obtain the electron-phonon matrix elements ==

We have finished the heavy simulations. Now it's time for the postprocessing. The first order of business is the reconstruction of the electron-phonon coupling matrix elements from the dvscf results and the electronic wavefunctions.

In order to do this, we will run the <code>lelphc</code> executable of the '''LetzElPhC''' code. We will run via command line using yambopy, although it will be instructive to have look at the <code>lelphc</code> input files later.

We run in the same directory where the Yambo <code>SAVE</code> is (remember than you can also virtually move it with a symbolic link).

Type:

 yambopy l2y

to see the help for the calculation. For example, if we want to do a serial run of LetzElPhC for bands from <math>n_i</math> to <math>n_f</math>, we should type:

 yambopy l2y -ph path/of/ph_input.in -b n_i n_f

Here <math>n_i</math> and <math>n_f</math> are integers representing the initial and final band indices. 

These should coincide with those used for the Bethe-Salpeter kernel, i.e. those specified in the <code>BSEBands</code> variable of the BSE input file (this is not strictly necessary, but certainly efficient since these calculations use a lot of disk space). 

For our system, we want to do a parallel calculation with 4 qpools and 2 kpools. In addition, we want to explicitly specify the path of the <code>lelphc</code> executable and avoid automatically deleting the LetzElPhC data. So we type:

 yambopy l2y -ph path/of/dvscf/mos2.dvscf -b 25 28 -par 4 2 -lelphc path/to/lelphc_exe --debug

At the end, check the <code>SAVE</code>

 ls SAVE/ndb.elph_gkkp*

to see that it has created the Yambo-compatible electron-phonon databases.

If you saved the <code>lelphc.in</code> input file, you can inspect it:

 # LetzElPhC input for yambo generated by yambopy
 nqpool      = 2
 nkpool      = 4
 start_bnd   = 25
 end_bnd     = 28
 save_dir    = ./SAVE
 kernel      = dfpt
 ph_save_dir = dvscf/ph_save
 convention = yambo

You can also run it as it is, but the code will generate the database <code>ndb.elph</code>. In order to convert it to the <code>ndb.elph_gkkp*</code> databases of Yambo, you still need a couple of lines of python using the Yambopy class <code>ConvertElectronPhononDB</code> in <code>yambopy/letzelph_interface/lelph2y.py</code>.

Notice the variable <code>convention=yambo</code>: what does it mean? At variance with QE and many other codes, Yambo uses the "backward" momentum transfer convention for electronic scatterings. That is, an electronic transition goes from band <math>n</math> and momentum <math>k-q</math> to band <math>m</math> and momentum <math>k</math>. In the "forward" momentum transfer convention (the more standard one), the transitions go from <math>nk</math> to <math>mk+q</math>. Therefore, this variable ensures that the electron-phonon coupling matrix elements are computed as <math>\langle mk|dV|nk-q\rangle</math>. This will have consequences also in the formulation of the ''exciton''-phonon coupling matrix element.

== 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>\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>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>g_{nm}^\mu (k,q)</math> are the electron-phonon coupling matrix elements obtained in step 6. As you can see, the exciton <math>\lambda</math> undergoes phonon-mediated scattering to state <math>\alpha</math> via phonon mode <math>\mu</math>. The scattering can happen for the hole (valence, <math>v</math>) or for the electron (conduction, <math>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 <code>excph.in</code>, 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>\lambda</math> with <code>ELPhExcStates</code>, the final exciton states <math>\alpha</math> with <code>ELPhExcSum</code> and the phonon modes <math>\mu</math> with <code>ElPhModes</code>. Here we consider the first four states at <math>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>q</math> point. We also include all the nine phonon modes of monolayer MoS2.

What about <code>LoutPath</code>? This variable controls the directory where the databases for the final-state excitons <math>\alpha</math> is located, which can be different from the directory with the initial-state excitons <math>\lambda</math> read as usual with the <code>-J</code> 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 <code>Lfull</code> in both cases.

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

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

If you check the output, you should find the <code>ndb.excph*</code> databases in the <code>excph</code> 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<sub>2</sub>, 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 <code>(3,4)</code> in fortran indexing or <code>(2,3)</code> 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>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 <code>analyse_excph.py</code> in which we first load the excph dabatases using the <code>YamboExcitonPhononDB</code> from yambopy. 
You can find this script in the yambopy directory, in <code>tutorials/exciton-phonon</code>. 
First, we select the exciton and phonon states to be included in <code>F_A</code>, 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>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 <code>YamboExcitonPhononDB</code>

 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 [https://wiki.yambo-code.eu/wiki/index.php?title=First_steps_in_Yambopy 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 <code>numpy</code> and <code>matplotlib</code> 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.

[[File:1L MoS2 ExcPh.png|500px|center]]

This can be checked against Fig. 2(d) of reference <ref name="chan2023" />, 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 ==

[[File:Luminescence scheme.png|250px|right]]

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:

<math>I^{(1)}_{PL}(\omega;T) \propto \frac{1}{N_q}\sum_{ s \mu \beta q} \left|\sum_\lambda\frac{D_\lambda \mathcal{G}^{\mu q}_{\alpha q,\lambda}}{E_{\alpha q}-E_\lambda -s\Omega_{\mu q}} \right|^2 N_{\alpha q}(T_{exc}) F^s_{\mu q}(T)\delta(\omega - \left[E_{\alpha q}-s\Omega_{\mu q}\right])</math>

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

The occupation function <math>F</math> is <math>F^s_{\mu q}(T)=n_{\mu q}(T)+\frac{1+s}{2}</math>. Here, <math>n(T)</math> is the temperature-dependent phonon Bose-Einstein occupation function. As it can be seen, <math>s=1</math> corresponds to processes of phonon ''emission'' (<math>\propto n(T)+1</math>), while <math>s=-1</math> corresponds to processes of phonon ''absorption'' (<math>\propto n(T)</math>). Therefore, <math>I^{(1)}_{PL}(\omega;T)</math> describes ''light'' emission by recombining excitons mediated by either ''phonon'' absorption or emission. 

The quantity <math> 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>T_{exc}</math> tends to coincide with the lattice temperature <math>T</math> more or less above 100 K, while at very low temperature (< 10 K), <math>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 <code>hbn.scf</code>:

 &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 <code>hbn.nscf</code>:

 &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 <code>hbn.dvscf</code>:

 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 <code>bse.in</code> (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
 <span style="color:red">Lkind= "Lfull"                    # [BSE] Lbar (default) / full</span>
 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
 %
 <span style="color:red">% BSEQptR</span>
 <span style="color:red"> 1 | 14 |                     # [BSK] Transferred momenta range</span>
 <span style="color:red">%</span>
 <span style="color:red">% BSEBands</span>
 <span style="color:red">   7 | 10 |                         # [BSK] Bands range</span>
 <span style="color:red">%</span>
 % 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 <code>bse_Lbar.in</code> with the following changes: 

 <span style="color:red">Lkind= "Lbar"                    # [BSE] Lbar (default) / full</span>
 <span style="color:red">% BSEQptR</span>
 <span style="color:red"> 1 | 1 |                     # [BSK] Transferred momenta range</span>
 <span style="color:red">%</span>

4b. So now we make a second BSE run in a different directory specified by <code>-J</code>. Here, we also pass to yambo the directory of the previous run as it includes the important screening databases <code>ndb.em1s*</code> 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 <code>lelphc</code> 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 <code>excph.in</code> using <code>yambo_ph -excph o</code>. 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
 <span style="color:red">ExcPhOptics                      # [R][EXCPH] Exciton-Phonon Optics</span>
 <span style="color:red">BoseTemp=10.000000         Kn    # Bosonic Temperature</span> 
 <span style="color:red">EXCTemp= 10.000000          Kn    # [EXCPH] Excitonic Temperature (for luminescence spectra)</span>
 % 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
 <span style="color:blue"># These are the plot parameters, same as in other parts of yambo</span> 
 % 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 <code>BoseTemp</code> and <code>EXCTemp</code> are the lattice and excitonic temperature, respectively.

Now we run the calculation with <code>yambo_ph</code>. Here, we read as <math>Q=0</math> excitons with <code>-J</code> 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 <code>excph</code> 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 <code>o-excph.pl_bse_ph_ass</code> 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

[[File:Luminescence plot.png|500px|center]]

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. 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 <ref name="zanfrognini2023" />.

== References ==

<references>
<ref name="toyozawa2003" >Toyozawa, Yutaka, and Chris Oxlade, ''Optical processes in solids'', [https://m.booksee.org/book/1121964?force_lang=en Cambridge University Press, (2003)]. </ref>
<ref name='lechifflart2023'>P. Lechifflart, F. Paleari, D. Sangalli, C. Attaccalite, ''First-principles study of luminescence in hexagonal boron nitride single layer: Exciton-phonon coupling and the role of substrate'', 
[https://doi.org/10.1103/PhysRevMaterials.7.024006 Phys. Rev. M, '''7''' (2), 024006 (2023)]; [https://arxiv.org/abs/2212.10407 arXiv2212.1047]</ref>
<ref name='cannuccia2019'>E. Cannuccia, B. Monserrat and C. Attaccalite, ''Theory of phonon-assisted luminescence in solids: Application to hexagonal boron nitride'', [https://doi.org/10.1103/PhysRevB.99.081109 Phys. Rev. B '''99''', 081109(R) (2019)]; [https://arxiv.org/abs/1807.11797 arXiv1807.11797]</ref>
<ref name='paleari2019'>F. Paleari et al., ''Exciton-Phonon Coupling in the Ultraviolet Absorption and Emission Spectra of Bulk Hexagonal Boron Nitride'', [https://doi.org/10.1103/PhysRevLett.122.187401 Phys. Rev. Lett. '''122''', 187401 (2019)]; [https://arxiv.org/abs/1810.08976 arXiv1810.089776] </ref>
<ref name='chen2020'>''Exciton-Phonon Interaction and Relaxation Times from First Principles'',
Hsiao-Yi Chen, Davide Sangalli, and Marco Bernardi, [https://doi.org/10.1103/PhysRevLett.125.107401  Phys. Rev. Lett. '''125''', 107401 (2020)]; [https://arxiv.org/abs/2002.08913 arXiv 2002.08913 (2020)]</ref>
<ref name="lechifflart2023_PhD">P. Lechifflart, ''Exciton-phonon coupling and phonon-assisted luminescence in hexagonal Boron Nitride nanostructures'', [https://hal.science/tel-04266805v1 PhD Thesis, University of Marseille (2023)]; [https://www.yambo-code.eu/wiki/images/5/54/These_final.pdf From the yambo website]</ref>
<ref name='paleari2019_PhD'>F. Paleari, ''First-principles approaches to the description of indirect absorption and luminescence spectroscopy: exciton-phonon coupling in hexagonal boron nitride'', [https://wwwen.uni.lu/research/fstm/dphyms/people/fulvio_paleari PhD thesis, University of Luxembourg (2019)]</ref>
<ref name='zanfrognini2023'>M. Zanfrognini et al., ''Distinguishing different stackings in layered materials via luminescence spectroscopy'', [https://doi.org/10.1103/PhysRevLett.131.206902 Phys. Rev. Lett. '''131''', 206902 (2023)]; [https://arxiv.org/abs/2305.17554 arXiv 2305.17554] </ref>
<ref name='marini2024'>G. Marini, M. Calandra, P. Cudazzo, ''Optical absorption and photoluminescence of single layer boron nitride from a first principles cumulant approach'', [https://doi.org/10.1021/acs.nanolett.4c00669 Nano Lett., '''24''', 20, 6017 (2024)]; [https://arxiv.org/abs/2402.03826 arXiv 2402.03826 (2024)]</ref>
<ref name='antonius2017'>G. Antonius, S. G. Louie, ''Theory of exciton-phonon coupling'', [https://doi.org/10.1103/PhysRevB.105.085111 Phys. Rev. B, '''105''', 085111 (2022)]; [https://arxiv.org/abs/1705.04245 arXiv1705.04245 (2017)]</ref>
<ref name='paleari2022'> F. Paleari, and A. Marini, ''Exciton-phonon interaction calls for a revision of the “exciton” concept'', [https://doi.org/10.1103/PhysRevB.106.125403 Phys. Rev. B, '''106''', 125403 (2022)]; [https://arxiv.org/abs/2205.02783 arXiv 2205.02783]</ref>
<ref name='cudazzo2020'> P. Cudazzo, ''First-principles description of the exciton-phonon interaction: A cumulant approach'', [https://doi.org/10.1103/PhysRevB.102.045136 Phys. Rev. B, '''102''', 045136 (2020)]; [https://orbilu.uni.lu/bitstream/10993/44769/1/main.pdf Open access pdf from Luxembourg University]</ref>
<ref name='chan2023'>Y-h Chan, J. B. Haber, M. H. Naik, J. B. Neaton, D. Y. Qiu, F. H. da Jornada, S. G. Louie, ''Exciton Lifetime and Optical Line Width Profile via Exciton–Phonon Interactions: Theory and First-Principles Calculations for Monolayer MoS2'', [https://doi.org/10.1021/acs.nanolett.3c00732 Nano Lett., '''23''', 9 (2023)]; [https://arxiv.org/abs/2212.08451 arXiv 2212.08451 (2023)]</ref>
<ref name='murali2025'>M. Nalabothula, S. Reichardt, L. Wirtz, ''Origin of Interlayer Exciton–Phonon Coupling in 2D Heterostructures'', [https://doi.org/10.1021/acs.nanolett.5c00355 Nano Lett., '''25''', 15 (2025)], [https://arxiv.org/abs/2407.16111 arXiv 2407.16111 (2025)]</ref>
</references>