Cumulant spectral function

In this section of the tutorial, we will compute the many-body electron/hole spectral function from cumulant expansion for the first-principles e-ph calculations. These calculations include the following sequential three calculation modes: calc_mode = 'spectral-re', calc_mode = 'spectral-se', calc_mode = 'spectral-cum'. The retarded cumulant approach is based on an exponential ansatz for the retarded Green’s function of an electronic state \(|nk\rangle\) (with band index n and momentum \(\mathbf{k}\)) in the time domain:

\[G^{R}_{nk}(t) = -i\theta(t) e^{-i\epsilon_{nk}t} e^{C_{nk}(t)}\]

where t is time, \(\theta(t)\) is the Heaviside step function, and \(\epsilon_{nk}\) is the noninteracting electron energy from DFT. The e-ph interactions are included in the cumulant \(C_{nk}(t)\), at the leading order is,

\[C_{nk}(t) = \int_{-\infty}^{\infty} d\omega \; \frac{\beta_{nk}(\omega)}{\omega^{2}} \left( e^{-i \omega t} + i \omega t - 1 \right),\]

where \(\beta_{nk}(\omega) \equiv \frac{\lvert \mathrm{Im}\,\Sigma_{nk}(\omega + \epsilon_{nk}) \rvert}{\pi}\) is the frequency-dependent imaginary part of the e-ph self-energy. Evaluating \(\beta_{nk}(\omega)\) on a dense frequency grid is the most time-consuming part of the calculation. The spectral function can be obtained for each state from the retarded Green’s function in the frequency domain, using \(A_{nk}(\omega) = -\frac{\mathrm{Im}\,G_{nk}^R(\omega)}{\pi}.\)

The methodology of first-principles electron-phonon cumulant spectral function is described in this paper: Phys. Rev. Research, 1, 033138, (2019)

Real-part on-shell Fan-Migdal self-energy
calc_mode = ‘spectral-re’

Directory: example11-sto-cumulant-spectral/perturbo/pert-spectral-re/ link

Computes: The on-shell real-part of the Fan-Migdal self-energy.

Note: There are soft phonon modes in SrTiO\(_3\) from DFPT calculations. Only for the tutorial purpose, we set a fictiously large phonon frequency cutoff to avoid the soft modes by setting phfreq_cutoff = 10.d0 meV. The best way is to include anharmonic effects, which can stablize the soft modes at finite temperatures. For full inclusion of anharmonic effects using TDEP, please follow this tutorial.

For the 'spectral-re' calculation, specify the following variables in the input file (pert.in):

prefix: the same prefix as in ‘prefix’_epr.h5
calc_mode: set to 'spectral-re'
band_min, band_max: band range
fklist: the filename of a file containing the list of momentum coordinates of electrons \(\mathbf{k}\) (see the format)
ftemper: the filename of a file containing temperature, chemical potential, and carrier concentration values (see the format)
phfreq_cutoff: cutoff below which a phonon will not be included in the calculation, here set to 10 meV phonon energy.
delta_smear: the broadening to model the Dirac delta function over energy difference, here chosen to be 10 meV.
sampling: sampling method for random \(\mathbf{q}\) points used in e-ph self-energy calculation. The default value is 'uniform', indicates sampling random \(\mathbf{q}\) points in the first BZ following uniform distribution. Another option is 'cauchy', sampling random \(\mathbf{q}\) points following Cauchy distribution, which is useful for polar materials. Note that random \(\mathbf{q}\) points from other importance sampling methods or \(\mathbf{q}\) points on regular MP grid is also possible, one just needs to pre-generate the \(\mathbf{q}\) points list to a file, and pass the file to perturbo.x via fqlist.
nsamples: number of random \(\mathbf{q}\) points sampled to compute the imaginary part of the e-ph self-energy for each \(\mathbf{k}\) point

Here is the input file (pert.in):

&perturbo
 prefix='sto'
 calc_mode = 'spectral-re'

 band_min = 1
 band_max = 3

 fklist = 'band.kpt'
 ftemper = 'sto.temper'

 phfreq_cutoff = 10.d0 
 delta_smear = 10.d0

 sampling = 'cauchy'
 nsamples = 100000
/

Before running perturbo.x, remember to put the following files in the run directory:

‘prefix’_epr.h5: here sto_epr.h5
ftemper: here sto.temper obtained in the 'setup' calculation, specifying the temperature, Fermi level and carrier concentration.
fklist: here band.kpt, specifying the \(\mathbf{k}\) points on which the spectral function is calculated.

The input files for the 'setup' calculation can be found in the following directory: “example11-sto-cumulant-spectral/perturbo/pert-setup”.

Run perturbo.x:

export OMP_NUM_THREADS=12
mpirun -n 1 perturbo.x -npools 1 -i pert.in > pert.out

This calculation usually takes a few minutes. We obtain the following two output files:

‘prefix’.resigma includes the on-shell real-part of self-energy of each electron states.
‘prefix’.resigma_mode includes the on-shell real-part of self-energy of each electron states, resolved by phonon mode.
‘prefix’_spectral-re.yml is the YAML output file containing information about the input and output parameters.

Note: One should always check the convergence of the e-ph self-energy with respect to the number of \(\mathbf{q}\) points (nsamples) and the smearing parameter (delta_smear). Check this paper for more detail.

Frequency-dependent imaginary part of the e-ph self-energy
calc_mode = ‘spectral-se’

Directory: example11-sto-cumulant-spectral/perturbo/pert-spectral-se/ link

Computes: The frequency-dependent imaginary part of the Fan-Migdal self-energy.

For the 'spectral-se' calculation, there are some extra variables in the input file than the 'spectral-re' calculation:

cum_inner_emin: Lower boundary for inner window for the frequency grid of the imaginary part of the self-energy. Here set to -0.2 eV.
cum_inner_emax: Upper boundary for inner window for the frequency grid of the imaginary part of the self-energy. Here set to 0.2 eV.
cum_outer_emin: Lower boundary for outer window for the frequency grid of the imaginary part of the self-energy. Here set to -0.6 eV.
cum_outer_emax: Upper boundary for outer window for the frequency grid of the imaginary part of the self-energy. Here set to 0.6 eV.
cum_de: The energy step in inner energy window for the frequency grid of the imaginary part of the self-energy. Here set to 0.4 meV.
cum_outer_np: The outer window energy step is cum_de * cum_outer_np. Here set to 4.

Here is the input file (pert.in):

&perturbo
prefix='sto'
calc_mode = 'spectral-se'

!--- cumu converges para------  
   cum_inner_emin = -0.2d0 !eV
   cum_inner_emax =  0.2d0 !eV
   cum_outer_emin = -0.6d0 !eV
   cum_outer_emax =  0.6d0 !eV 
   cum_de   =  0.4d0      !meV  
   cum_outer_np = 4
!------------------------------

 band_min = 1
 band_max = 3
 fklist = 'band.kpt'
 ftemper = 'sto.temper'
  
 phfreq_cutoff = 10.d0 
 delta_smear = 10.d0
 sampling = 'cauchy'
 nsamples = 100000
/

Before running perturbo.x, remember to put the following files in the run directory:

‘prefix’_epr.h5: here sto_epr.h5
ftemper: here sto.temper obtained in the 'setup' calculation, specifying the temperature, Fermi level and carrier concentration.
fklist: here band.kpt, specifying the \(\mathbf{k}\) points on which the spectral function is calculated.

Note: The ftemper and fklist files are the same as the ones used in the 'spectral-re' calculation. To generate these files, run the 'setup' calculation using input files found here.

Run perturbo.x:

export OMP_NUM_THREADS=12
mpirun -n 1 perturbo.x -npools 1 -i pert.in > pert.out

This calculation usually takes a few minutes for a single \(\mathbf{k}\) point. We obtain the following two output files:

‘prefix’_selfenergy.h5 , an hdf5 file, includes the frequency-dependent imaginary part of self-energy of each electron states.
‘prefix’_spectral-se.yml is the YAML output file containing information about the input and output parameters.