In this section, we describe the ultrafast carrier dynamics. In contrast to other PERTURBO calculation modes, here we solve the real-time Boltzman transport equation (rt-BTE):

\[\frac{\partial f_{n\mathbf{k}}(t)}{\partial t} = \frac{e\mathbf{E}}{\hbar}\cdot \nabla _\mathbf{k}f_{n\mathbf{k}}(t)+\mathcal{I}^{e-\mathrm{ph}}[f_{n\mathbf{k}}(t)],\]

where \(f_{n\mathbf{k}}(t)\) is the time-dependent electronic occupation of a Bloch state with band index \(n\) and crystal momentum \(\mathbf{k}\), and \(e\) is the electronic charge. The first term on the right-hand side, called the advection term, is proportional to the external electric field \(\mathbf{E}\) and is responsible for electron drift. The second term is the \(e\)-ph collision integral \(\mathcal{I}^{e-\mathrm{ph}}\) accounting for phonon absorption and emission processes.

First, we will consider the zero field ultrafast dynamics (\(\mathbf{E}=0\)), then the high-field dynamics (full rt-BTE) will be considered.

Zero-field ultrafast dynamics
calc_mode = ‘dynamics-run’

For the ultrafast dynamics, one needs first to perform the 'setup' calculation. So, we assume that the user has already performed this calculation. From the 'setup' calculation, we retain the following files necessary for the dynamics calculations: ‘prefix’.temper and ‘prefix’_tet.h5.

In a typical zero-field dynamics, at \(t=0\), one excites the electronic system with a Gaussian pulse and then follows the electron relaxation due to the electron-phonon scattering. The excited carrier concentration is preserved during the real-time simulation and the integral of the excitation pulse corresponds to the excited carrier density.

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

  • preifx: the same prefix as in ‘prefix’_epwan.h5
  • calc_mode: set to 'dynamics-run'
  • boltz_kdim: \(\mathbf{k}\) grid for electrons, here we use a 80x80x80 grid
  • boltz_qdim: \(\mathbf{q}\) grid for phonons, specify if it is different from the \(\mathbf{k}\) grid, here we use the same \(\mathbf{q}\) grid as \(\mathbf{k}\) grid
  • boltz_emin, boltz_emax: energy window (in eV units), use the same as in the 'setup' calculation (here, 6.4 and 7.4 eV)
  • band_min, band_max: band range
  • ftemper: the filename of a file containing temperature, chemical potential, and carrier concentration values (see the format)
  • time_step: simulation time step \(\Delta t\), set to its typical value, 1 fs
  • boltz_nstep: total number of time steps, set to 50; here we perform a relatively short simulation of 50 fs
  • output_nstep: an optional variable to shorten the output; the output time step \(\Delta t_{out}\) is determined in the following way: \(\Delta t_{out} = \texttt{output_nstep}\times \Delta t\)
  • solver: BTE solver type, set to 'euler', here we use the Euler first order solver of BTE
  • boltz_init_dist: set to 'gaussian', we select the Gaussian initial distribution. To restart the simulation, specify boltz_init_dist='restart', then the distribution of the last step from the previous simulation will be used.
  • boltz_init_e0: in this example, the Gaussian distribution is centered around 7.0 eV
  • boltz_init_smear: we select a 40 meV smearing
  • phfreq_cutoff: we select a 1 meV phonon energy cutoff
  • delta_smear: the broadening to model the Dirac delta function is chosen to 8 meV

Here is the input file (pert.in):

&perturbo
 prefix      = 'si'
 calc_mode   = 'dynamics-run'

 boltz_kdim(1) = 80
 boltz_kdim(2) = 80
 boltz_kdim(3) = 80

 boltz_emin = 6.4
 boltz_emax = 7.4
 band_min = 5
 band_max = 6

 ftemper  = 'si.temper'

 time_step   = 1 !fs
 boltz_nstep = 50
 output_nstep = 2
 solver = 'euler'

 boltz_init_dist = 'gaussian'
 boltz_init_e0 = 7.0 ! eV
 boltz_init_smear = 40 !meV

 tmp_dir = "./tmp"
 phfreq_cutoff = 1 ! meV
 delta_smear = 8 ! meV
/

In this example, we calculate the evolution of the electron distribution. In order to perform the hole dynamics, set the parameter hole to true.

Run perturbo.x (remember to link or copy ‘prefix’_epwan.h5 in the current directory):

$ export OMP_NUM_THREADS=4
$ mpirun -n 8 <perturbo_bin>/perturbo.x -npools 8 -i pert.in > pert.out

We obtain the ‘prefix’_cdyna.h5 HDF5 output file (this file can be also found in the “References” directory). This file contains all the necessary output information about the performed simulation. This file is organized as follows:

  • band_structure_ryd: electronic bandstructure in Ry; each column corresponds to the band index \(n\) \((~\texttt{band_min}\leq n \leq \texttt{band_max})\)
  • dynamics_run_[i]: an HDF5 group that contains information about the i th simulation.
    If the simulation was restarted (boltz_init_dist='restart') one or more times, one will have several dynamics_run_[i] groups, otherwise, only dynamics_run_1 will be present. A group dynamics_run_[i] is structured as follows:
    • num_steps: the number of output time steps (taking into account output_nstep), can be different for different dynamics_run_[i]

    • snap_t_0: \(f_{n\mathbf{k}}(t_0)\)
      \(\vdots\)
    • snap_t_[j]: the distribution function \(f_{n\mathbf{k}}\) for time \(t_j\): \(t_j = t_0 + j\Delta t_{out}\), where \(\Delta t_{out}\) is the output time step. Each column of the array corresponds to the band index \(n\)
      \(\vdots\)
    • snap_t_[num_steps]: \(f_{n\mathbf{k}}(t_{\texttt{num_steps}})\)

    • time_step_fs: the output time step \(\Delta t_{out}\) (can be different for different dynamics_run_[i])
  • num_runs: total number of performed simulations (corresponds to the number of dynamics_run_[i] groups).


The ‘prefix’_cdyna.h5 file structure can be schematically represented as follows:

diagram_hdf5_dynamics-run

The HDF5 files can be easily processed by Python package h5py. As an example, we present here a simple Python script that visualizes the distribution function for the time \(t_5\) of the simulation and for the first band (in the selected band range):

#!/usr/bin/env python3
import h5py
import matplotlib.pyplot as plt

prefix='si'
snap_number=5
band_index=0

# load the HDF5 file
h5file = h5py.File(prefix+'_cdyna.h5', 'r')

# get the data
ryd2ev = h5file['band_structure_ryd'].attrs['ryd2ev']
energy_ev = h5file['band_structure_ryd'][:,band_index] * ryd2ev
dist_func = h5file['dynamics_run_1']['snap_t_'+str(snap_number)][:,band_index]
h5file.close()

# plot the data
plt.plot(energy_ev,dist_func,marker='o',linestyle='')
plt.xlabel('Energy (eV)')
plt.ylabel('Distribution function')
plt.show()

In order to postprocess this file using perturbo.x, see the next section.

Dynamics post-processing
calc_mode = ‘dynamics-pp’

In this section we aim to calculate the Brillouin zone-averaged energy-dependent carrier population \(\bar{f}(E,t)\). Having calculated the distribution function \(f_{n\mathbf{k}}(t)\), one can find \(\bar{f}(E,t)\) in the following way:

\[\bar{f}(E,t) = \sum_{n\mathbf{k}} f_{n\mathbf{k}}(t) \delta(\epsilon_{n\mathbf{k}}-E).\]

The integral of \(\bar{f}(E,t)\) over the energy gives the number of carriers per unit cell as a function of time.

In order to calculate the \(\bar{f}(E,t)\) quantity, one needs to have all the files required for the calc_mode='dynamics-run' calculation (previous section) and the HDF5 output file ‘prefix’_cdyna.h5 from the dynamics-run calculation. To perform the postprocessing, use a similar to the previous section input file, but change the calculation mode to calc_mode='dynamics-pp'. Run perturbo.x (remember to link or copy ‘prefix’_epwan.h5 in the current directory):

$ mpirun -n 1 <perturbo_bin>/perturbo.x -npools 1 -i pert.in > pert.out

On the output, we obtain the following files:

  • si_popu.h5: an HDF5 file that contains all the necessary information for \(\bar{f}(E,t)\)
  • si_cdyna.dat: an ASCII file containing the number of carriers per unit cell as a function of time

The si_popu.h5 HDF5 file is organized as follows:

  • energy_distribution: a group that contains the populations for all the time instants of the dynamics-run simulation

    • popu_t1: \(\bar{f}(E,t_1)\)
      \(\vdots\)
    • popu_t[j]: the carrier population \(\bar{f}(E,t_j)\) at time \(t_j\)
      \(\vdots\)
    • popu_t[num_steps+1]: \(\bar{f}(E,t_{\texttt{num_steps+1}})\)

  • energy_grid_ev: the grid of energies in eV; the number of energy grid points is given by \(\frac{ \texttt{emax} - \texttt{emin} }{ \texttt{boltz_de} }+\texttt{3}\)
  • times_fs: the array of time instants in fs


The si_popu.h5 HDF5 file can be schematically represented as follows:

diagram_hdf5_dynamics-pp

Similarly to the previous section, we provide here a simplistic Python script showing an example how to manipulate this HDF5 file. For example, to plot the electron population for the time \(t_{25}\), run:

#!/usr/bin/env python3
import h5py
import matplotlib.pyplot as plt

prefix='si'
snap_number=25

# load the HDF5 file
h5file = h5py.File(prefix+'_popu.h5', 'r')

# get the data
energy_ev = h5file['energy_grid_ev'][()]
population = h5file['energy_distribution']['popu_t'+str(snap_number)][()]
h5file.close()

# plot the data
plt.plot(energy_ev,population,marker='o',linestyle='')
plt.xlabel('Energy (eV)')
plt.ylabel('Electron population')
plt.show()

It is also convenient to postprocess and visualize the data in HDF5 file using other high level languages, such as Julia. For example, the following Julia script does the same thing as the above Python script:

using HDF5, Plots

prefix = "si"
fname = prefix * "_popu.h5"
snap_number = 25

# read the data
energy_ev = h5read(fname, "energy_grid_ev")
population = h5read(fname, "energy_distribution/popu_t"*string(snap_number))

# plot
plot(energy_ev, population, xlabel="Energy (eV)", ylabel="Electron population")