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 implemented in PERTURBO can be found in these papers:

*Nano Lett.*

**17**, 5012-5019, (2017),

*Comput. Phys. Commun.*

**264**, 107970, (2021)

**high-field**ultrafast dynamics implemented in PERTURBO can be found in this paper:

*Phys. Rev. B*

**104**, L100303, (2021)

## Zero-field ultrafast dynamics

*calc_mode = ‘dynamics-run’*

**Computes:**Ultrafast hot carrier dynamics via the time-dependent Boltzmann transport equation: set an initial carrier distribution and calculate its evolution in time at zero external fields.

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*.

**Note:**For the ultrafast dynamics we will need a larger energy window. So, the

`'setup'`

calculation should be performed with `boltz_emin = 6.4`

and `boltz_emax = 7.4`

.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.

**Note:**For the ultrafast dynamics, the carrier density is determined by the initial condition and not by the

*‘prefix’.temper*file, in contrast to other calculation modes.

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*simulation.^{ th}

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:

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’*

**Computes:**Postprocessing of the ultrafast dynamics calculations: carrier population as a function of energy and time.

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:

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:

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")
```