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*):

- prefix: the same prefix as in
*‘prefix’_epr.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, we use the boltz_nstep parameter to control the number of time steps (here, set to 50). Alternatively, the following parameters can be used to control the number of time steps:

- boltz_acc_thr: drift acceleration threshold in units of cm s
^{-2}; if specified, the real-time simulation will stop when at least 10 last iterations had a drift acceleration lower than the threshold - boltz_nstep: only applied if boltz_acc_thr is specified; minimum number of time steps taken before real-time simulation stops (even if the drift acceleration threshold criteria has already been met)

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’_epr.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:

We also obtain an outputted YAML file, here called *si_dynamics-run.yml*, containing the inputs and outputs of the `'dynamics-run'`

calculation. The sections that are common to all Perturbo calculations are described here.

An example of the fields unique to the `'dynamics-run'`

calculation is given below.

```
grid:
number of tetrahedra selected: 380160
number of reducible k-points: 72165
number of irreducible k-points: 1797
scatter parallelization:
MPI task index:
1:
number of q-points: 254451
number of (k,q) pairs: 48597938
2:
number of q-points: 254484
number of (k,q) pairs: 48606145
......
dynamics-run:
```

Because most data is outputted to the HDF5 files, very little data is outputted to the *si_dynamics-run.yml* file. The *grid* field gives information on the number of tetrahedra, reducible \(\mathbf{k}\) points, and irreducible \(\mathbf{k}\) points in the energy window of interest. The specific \(\mathbf{k}\) points and tetrahedra are given in the *‘prefix’_tet.h5* file. The *scatter parallelization* field gives information on the number of \(\mathbf{k}\) points and their distribution across the MPI tasks. The *dynamics-run* field is empty, as the results are outputted to the *‘prefix’_cdyna.h5* HDF5 file.

We may next use Perturbopy to export the data to Python for postprocessing. For more details, see the Perturbopy website, which explains how Perturbopy processes the HDF5 files and can be used to quickly generate plots. Alternatively, users can process HDF5 files using their own Python scripts with 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 *dynamics-pp* section.

We may also want to run the calculation including only certain phonon modes to determine the effects of different modes. To do this, we can use the optional parameter ph_mode_exclude_ranges, which excludes phonon modes in the rt-BTE by setting their `|g`

factors to 0. For example:^{2}|

- to exclude modes from 2 to 5, write in the input file:
`ph_mode_exclude_ranges(1,:)=2,5`

. The index 1 reflects that this is the first range to exclude. - to also exclude modes from 7 to 10, add another line in the input file:
`ph_mode_exclude_ranges(2,:)=7, 10`

. The index 2 reflects that this is the second range to exclude. - to also exclude just one mode (say, mode 12), write:
`ph_mode_exclude_ranges(3,:)=12, 12 .`

. The index 3 reflects that this is the third range to exclude.

Note that the maximum number of ranges is 30, and the indices of the phonon modes to exclude must be smaller than the total number of modes in a system.

## 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’_epr.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:

We also obtain an outputted YAML file, here called *si_dynamics-pp.yml*, containing the inputs and outputs of the `'dynamics-pp'`

calculation. The sections that are common to all Perturbo calculations are described here.

An example of the fields unique to the `'dynamics-pp'`

calculation is given below.

```
grid:
number of tetrahedra selected: 380160
number of reducible k-points: 72165
number of irreducible k-points: 1797
dynamics-pp:
time units: fs
concentation units: num. of carriers/u.c.
number of snaps: 25
time:
- 0.0000000000000000E+00
......
- 5.0000000000000000E+01
concentration:
- 1.1419980691444462E-03
......
- 1.1419968066740211E-03
```

Again, the *grid* field gives information on the number of tetrahedra, reducible \(\mathbf{k}\) points, and irreducible \(\mathbf{k}\) points in the energy window of interest. The specific \(\mathbf{k}\) points and tetrahedra are given in the *‘prefix’_tet.h5* file. The *dynamics-pp* field gives, under the *time* fields and *concentration* fields respectively, the time in fs and carrier concentration per unit cell (integral of \(\bar{f}(E,t)\) over the energy) at each of the time steps (here, there are 25).

We may next use Perturbopy to export the data to Python for postprocessing. For more details, see the Perturbopy website, which explains how Perturbopy processes the output and can be used to quickly generate plots. Alternatively, users can process HDF5 files using their own scripts (i.e. Python, Julia, or other high-level languages). For example, we provide here a simplistic Python script showing an example how to manipulate this HDF5 file. 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")
```

## High-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 nonzero external fields.

The steps for running high-field ultrafast dynamics are similar to those for zero-field ultrafast dynamics, with the following modifications to the input script.

To specify the external electric field in V cm^{-1}, we modify the following boltz_efield parameter by adding a few lines to the previous input file. This parameter is three-dimensional. Below, we modify only the *x* component of the external field, and leave the other two components as zero.

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_efield(1) = 100.0
boltz_efield(2) = 0.0
boltz_efield(3) = 0.0
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
/
```

We obtain the *‘prefix’_cdyna.h5* HDF5 output file, which is organized almost the same as the file shown above. However, it has the following additional field under the `dynamics_run_[i]`

group:

`dynamics_run_[i]`

`efield`

: the*x*,*y*, and*z*components of the external electric field in V cm^{-1}

We also obtain the *‘prefix’_dynamics-run.yml* output file, which is identical to the example above.

We may next use Perturbopy to export the data to Python for postprocessing. For more details, see the Perturbopy website

## High-field ultrafast dynamics

*calc_mode = ‘dynamics-pp’*

The steps for running high-field ultrafast dynamics post-processing are identical to those for zero-field ultrafast dynamics. Please refer to the section above for details.

After running the calculation, we obtain the similar *‘prefix’_popu.h5*, described above. The *prefix_dynamics-pp.yml* file described above is also analogous, but has additional fields to specify the drift velocities.

```
dynamics-pp:
time units: fs
concentation units: num. of carriers/u.c.
number of snaps: 25
electric field: [ 100.00000, 0.00000, 0.00000, ]
velocity units: cm/s
time:
- 0.0000000000000000E+00
......
- 5.0000000000000000E+01
concentration:
- 1.1419980691444460E-03
......
- 1.1419975237966593E-03
velocity:
- 3.4859250185988308E-03
......
- -1.5315116904290144E+04
```

Like in the zero-field case, the *dynamics-pp* field gives, under the *time* fields and *concentration* fields respectively, the time in fs and carrier concentration per unit cell at each of the time steps. There is now an additional field *velocity*, which gives the drift velocity in cm s^{-1} unitsat each time. It also gives the electric field, though please note that this echos the *boltz_efield* input in the `'dynamics-pp'`

input file and may not be meaningful if this value differs from that of the `'dynamics-run'`

input file.

We may next use Perturbopy to export the data to Python for postprocessing. For more details, see the Perturbopy website