Simulating Pump-Probe Spectroscopy

by Ronaldo Rodrigues Pela for exciting oxygen

Purpose: In this tutorial, you will learn how to employ real-time time-dependent density-functional theory (RT-TDDFT) to simulate a simple pump-probe experiment. We consider two-dimensional MoS₂ to illustrate.

---

0. Define relevant environment variables

Read the following paragraphs before starting with the rest of this tutorial!

Before starting, be sure that relevant environment variables are already defined as specified in How to set environment variables for tutorials scripts. Here is a list of the scripts which are relevant for this tutorial with a short description.

• EXECUTE-single.sh: (Bash) shell script for running a single exciting calculation.
• PLOT-multitask.py: Python script for making plots from output files of RT-TDDFT.

From now on the symbol $ will indicate the shell prompt.

Requirements: Bash shell. Python numpy, lxml, matplotlib.pyplot, and sys libraries.

Important note: All input parameters that will appear will be given in atomic units (a.u.)! For this tutorial, it can be useful to remind the conversion between a.u. and the SI for times: 1 a.u. = 0.02418884254 fs.

---

1. Theoretical background

In a pump-probe experiment, there are basically two laser pulses that shine a sample with different goals. The first pulse (pump) excites the sample during a specific time, generating a non-equilibrium state. The second pulse, usually much weaker than the first, probes specific properties, such as optical constants. With this technique, it is possible to analyze the impact of excitations on the optical properties, and how the systems relax back to the equilibrium.

In exciting, it is possible to simulate this kind of experiment, employing RT-TDDFT. In one possible scenario, we first consider a pump field described by a cosine modulated by a sin squared. This field mimics a gaussian envelope, which is quite common in experiments. The expression for the vector potential $\mathbf{A}_\textrm{pump}(t)$ is

(1)

\begin{align} \mathbf{A}_\textrm{pump}(t) = \mathbf{A}_0 f(t) \cos ( \omega t + \phi ), \end{align}

where $f(t)$, the envelope function, is given as

(2)

\begin{align} f(t) = \left\{ \begin{array}{ll} 0 , & t \le t_0 \quad\mbox{ or }\quad t \ge t_0 + w \\ \sin^2( \pi(t-t_0)/w)\quad & t_0 \le t \le t_0 + w \\ \end{array} \right. \end{align}

The following figure illustrates the considered field

After some time $t$ when the pump is not acting any longer, we can consider a delta-kick as the probe

(3)

\begin{align} \mathbf{A}_\textrm{probe}(t) = -c\,\mathbf{E}'_0\, u( t-t_\textrm{probe} ), \end{align}

where $u( t-t_\textrm{probe} )$ is a step-function applied at time $t=t_\textrm{probe}$. The electric field of the probe is

(4)

\begin{align} \mathbf{E}_\textrm{probe}(t)=-\frac{1}{c}\frac{d}{dt}\mathbf{A}_\textrm{probe}(t) = \mathbf{E}'_0\, \delta( t-t_\textrm{probe} )\:. \end{align}

In this tutorial, we will analyze how the imaginary part of the dielectric function (related to the absorption spectrum) of the sample excited by the pump differs from its ordinary counterpart (evaluated starting with the sample on the ground state). With this purpose, we consider the general relationship between the Fourier-transforms of the current density $J$ and the electric field $E$:

(5)

\begin{align} \sigma_{\alpha\beta}(\omega) = \frac{J_\alpha(\omega)}{E_\beta(\omega)}, \quad \quad\varepsilon_{\alpha\beta}(\omega) = \delta_ {\alpha\beta}+ \frac{4\pi \mathrm{i}\,\sigma_{\alpha\beta}(\omega)}{\omega}, \end{align}

in which $\sigma$ and $\varepsilon$ are the optical conductivity and the dielectric tensors, respectively; the indexes $\alpha$, $\beta$ denote the cartesian directions $x$, $y$, or $z$ [1].

To simulate the pump-probe experiment, we can evaluate the dielectric function as indicated in Eq. (5), but taking the current density as the difference between the values after the pump and the probe ($J_\textrm{pump-probe}$) and the pump pulse $J_\textrm{pump}$, as obtained from two separate calculations. And as comparison, we consider a third calculation just with the probe pulse, from which we can extract the dielectric function from the system starting on the groundstate.

Please refer to Ref.[2] for one different example of simulation of pump-probe spectroscopy with exciting.

---

2. Ground-state calculation

i) Preparation of the input file

The first step is to create a directory for the system that we want to investigate, MoS2-pump_probe and to enter it.

$ mkdir MoS2-pump_probe
$ cd MoS2-pump_probe

As initial condition for the time evolution done in RT-TDDFT, we need the electron density and potential from a ground-state calculation. Therefore, we create the file input.xml that could look like the following:

<input>
 
   <title>MoS2</title>
 
   <structure speciespath="$EXCITINGROOT/species/">
      <crystal scale="6.0274">
         <basevect>  0.5  0.86602540378   0.0 </basevect>
         <basevect> -0.5  0.86602540378   0.0 </basevect>
         <basevect>  0.0  0.00000000000  10.0 </basevect>
      </crystal>
      <species speciesfile="Mo.xml" rmt="2.3">
         <atom coord="0.666667  0.666667  0.0"/>
      </species>
      <species speciesfile="S.xml"  rmt="2.05">
         <atom coord="0.333333  0.333333 -0.04938"/>
         <atom coord="0.333333  0.333333  0.04938"/>
      </species>
   </structure>
 
   <groundstate
      do="fromscratch"
      xctype="LDA_PW"
      ngridk="12 12 1"
      rgkmax="5.0d0"
      nempty="5"
      epsengy="1.d-8">
   </groundstate>
 
</input>

N. B.: Do not forget to replace in the input.xml the string "$EXCITINGROOT" by the actual value of the environment variable $EXCITINGROOT. You can do this by directly editing the input.xml file or by using the following command:

$ SETUP-excitingroot.sh

ii) Running the ground-state calculation

You can start the calculation by invoking the script EXECUTE-single.sh.

$ EXECUTE-single.sh GS

After the calculation is finished, move to the subdirectory GS.

$ cd GS

You can check the output files, especially the main output file INFO.OUT, for general information. In the case of a successfully finished calculation, the last lines of the INFO.OUT will contain the message

...
================================================================================
| EXCITING OXYGEN stopped                                                      =
================================================================================

Check if this is indeed the case and if the EFERMI.OUT and STATE.OUT files are present. They contain the Fermi level and the converged electron density and potential, respectively, and are the starting point of the RT-TDDFT calculation.

Please note: To obtain reliable ground-state density and potential, it is necessary to perform convergence tests with respect to the k-point mesh (parameter ngridk) and the size of the basis set (parameter rgkmax). For details see Simple convergence tests.

---

3. Performing the pump-probe calculations

i) First case: pump-probe

Move to the previous directory, create a new one named pump-probe, and enter it:

$ cd ../
$ mkdir pump-probe
$ cd pump-probe

You can prepare the input file (input.xml) of a RT-TDDFT calculation starting from the existing ground-state input file. Copy the following files to the current directory. You can perform all of this with the following commands:

$ cp ../GS/{input.xml,EFERMI.OUT,STATE.OUT} .

Edit the file input.xml and avoid the explicit ground-state calculation with the attribute do = "skip" inside the element groundstate.

...
   <groundstate 
      do="skip"
      ... >
      ...
   </groundstate>
...

Also, to perform the "RT-TDDFT" calculation, add the element xs (as shown below) the input file inside the input element.

...
   <xs
      xstype="RT-TDDFT"
      ngridk="6 6 1"
      vkloff="0.01 0.02 0.004"
      nempty="5"
      nosym="true"
      reducek="false">
 
      <realTimeTDDFT
         propagator="AETRS"
         timeStep="0.5d0"
         endTime="500.d0"
         vectorPotentialSolver="improvedeuler">
         <laser fieldType="total">
            <sinSq amplitude="20.d0" omega="0.5d0" phase="0.d0"
               t0="0.1d0" pulseLength="80.d0" direction="x" />
            <kick amplitude="0.1d0" t0="100.d0" width="0.d0" direction="x" />
         </laser>
      </realTimeTDDFT>
   </xs>
...

Please check all relevant parameters for a "RT-TDDFT" calculation in Input Reference.

Having modified the input file as described above, you can launch the RT-TDDFT calculation inside the subdirectory pump-probe using

$ time exciting_smp

Be aware that, depending on the computer you are using, this calculation may take several minutes up to some hours.

ii) Second case: pump

Now, we can study the effect of the pump alone on MoS₂. Move to the previous directory, create a new one named pump, and enter it:

$ cd ../
$ mkdir pump
$ cd pump

Copy the files input.xml, EFERMI.OUT and STATE.OUT from the groundstate calculation to the current directory:

$ cp ../GS/{input.xml,EFERMI.OUT,STATE.OUT} .

As done before, edit the file input.xml to avoid the explicit ground-state calculation by setting do = "skip" inside the element groundstate.

...
   <groundstate 
      do="skip"
      ... >
      ...
   </groundstate>
...

Now, you need to add the element xs to input.xml as shown below:

...
   <xs
      xstype="RT-TDDFT"
      ngridk="6 6 1"
      vkloff="0.01 0.02 0.004"
      nempty="5"
      nosym="true"
      reducek="false">
 
      <realTimeTDDFT
         propagator="AETRS"
         timeStep="0.5d0"
         endTime="500.d0"
         vectorPotentialSolver="improvedeuler">
         <laser fieldType="total">
            <sinSq amplitude="20.d0" omega="0.5d0" phase="0.d0"
               t0="0.1d0" pulseLength="80.d0" direction="x" />
         </laser>
      </realTimeTDDFT>
   </xs>
...

Now you are ready to execute the RT-TDDFT calculation inside the subdirectory pump using

$ time exciting_smp

Have in mind that depending on your computational resources, this calculation may take some minutes or even hours.

iii) Third case: probe

Finally, we can carry out a calculation to obtain the dielectric function of MoS₂ starting from the groundstate. Similar as you did before, you need to go to the previous directory, create a new one named probe, and enter it:

$ cd ../
$ mkdir probe
$ cd probe

You can again copy the files input.xml, EFERMI.OUT and STATE.OUT from the groundstate calculation to the current directory:

$ cp ../GS/{input.xml,EFERMI.OUT,STATE.OUT} .

As done before, edit the file input.xml to skip the ground-state calculation:

...
   <groundstate 
      do="skip"
      ... >
      ...
   </groundstate>
...

Now, add the element xs to input.xml as shown below:

...
   <xs
      xstype="RT-TDDFT"
      ngridk="6 6 1"
      vkloff="0.01 0.02 0.004"
      nempty="5"
      nosym="true"
      reducek="false">
 
      <realTimeTDDFT
         propagator="AETRS"
         timeStep="0.5d0"
         endTime="500.d0"
         vectorPotentialSolver="improvedeuler">
         <laser fieldType="total">
            <kick amplitude="0.1d0" t0="100.d0" width="0.d0" direction="x" />
         </laser>
      </realTimeTDDFT>
   </xs>
...

Then, perform the RT-TDDFT calculation (inside the subdirectory probe) using

$ time exciting_smp

As warned before, this calculation may take more some minutes or hours, depending on the computational power that you have.

iv) Analyzing the dielectric function

First, it is convenient to obtain the difference $J_\textrm{pump-probe}-J_\textrm{pump}$. We first move to the parent directory, MoS2-pump_probe/, and then preprocess the files JIND.OUT in the folders pump-probe/ and pump/ with the command:

$ cd ../
$ PLOT-multitask.py --preprocess --sub -f pump-probe/JIND.OUT pump/JIND.OUT -o DIFF.OUT

With the previous command, the difference $J_\textrm{pump-probe}-J_\textrm{pump}$ is stored in the file DIFF.OUT. It assumes that we take the x components, otherwise, we would need to add --y or --z.

Now, you can plot the current densities $J_\textrm{pump}$, $J_\textrm{probe}$, $J_\textrm{pump-probe}$ and the difference $J_\textrm{pump-probe}-J_\textrm{pump}$ with [3]:

$ PLOT-multitask.py --jind -f pump/JIND.OUT probe/JIND.OUT pump-probe/JIND.OUT DIFF.OUT --legend 'lower right'

Your figure should be similar to the one given bellow:

Now, we can compare the imaginary part dielectric function in the pump-probe case, and compare it with the probe case. First, we determine the dielectric function for the probe case:

$ cd probe/
$ PLOT-multitask.py --preprocess --get_eps --xx -f AVEC.OUT -f JIND.OUT --wcut 0.05 -o EPSILON.OUT

Now, we move back to parent directory, and obtain the dielectric function for the pump-probe:

$ cd ../
$ PLOT-multitask.py --preprocess --get_eps --xx -f probe/AVEC.OUT -f DIFF.OUT --wcut 0.05 -o EPSILON-DIFF.OUT

And, finally, we can plot the imaginary part of both dielectric functions in the same Figure using:

$ PLOT-multitask.py --imag_eps -f probe/EPSILON.OUT EPSILON-DIFF.OUT --xlim 5 30 --ylim 0 1

Your figure should be similar to the one bellow:

Be aware that the calculations shown here employ parameters which are not converged, and should be considered only for introductory purposes. To improve the quality of the calculation, special attention needs to be paid to the following parameters:

• ngridk
• rgkmax
• timeStep and endTime

---

Bibliography

1. Note that to obtain $\sigma_{\beta\alpha}$ it is convenient to take electric field as a delta-kick $E_0\,\delta(t)$ in a specific direction $\alpha$ such that $E_\alpha(\omega) = E_0$.

2. Ronaldo Rodrigues Pela, Claudia Draxl. "All-electron full-potential implementation of real-time TDDFT in exciting" (link).

3. More information about the script PLOT-multitask.py can be found here.