Inelastic Electron Spectroscopy of an H2 molecule placed between 1D Au chains


This tutorial is focused on simulating inelastic electron tunneling spectroscopy (IETS) [Ree08] for a device consisting of two one-dimensional (1D) gold electrodes and an H2 molecule placed in between. After DFT calculations of the dynamical matrix and the Hamiltonian derivatives with respect to vibrations of the H2 molecule, IETS is computed and analyzed based on the lowest-order expansion (LOE) [FPBJ07] and extended LOE (XLOE) [LCF+14] approximations.


It is assumed that you are familiar with Atomistix ToolKit and Virtual NanoLab. If not, go through the basic VNL and ATK Tutorials first.



We consider electron transport through a device consisting of an H2 molecule clamped between 1D gold electrodes. In a spectroscopy experiment, if the applied bias is comparable to the energy of one of the vibrational modes of the H2 molecule, a new inelastic transmission channel opens in which the transmission of electrons is mediated by the emission of a phonon. The coupling process is due to the coupling between electronic and vibrational degrees of freedom, that is, to the electron-phonon coupling. At a voltage where an inelastic channel opens, the derivative of the conductance, i.e., the second-derivative of the current has a peak. In experiments, IETS is defined by

\[{\rm IETS}\equiv \frac{dI^2/dV^2}{dI/dV},\]

and is measured to investigate, e.g., the electron-phonon coupling strength, heating, configuration, and local doping.

In the tutorial, you will first compute the dynamical matrix and the Hamiltonian derivatives with respect to vibrations of the H2 molecule by DFT calculations, and then calculate and analyze the IETS signal based on two different approximations: LOE and XLOE. The rest of this tutorial is hence organized as follows:

  1. Device setup,
  2. Calculation of the dynamical matrix and the Hamiltonian derivatives,
  3. Calculation of IETS based on LOE and XLOE,
  4. Analysis of IETS based on the vibrational modes of the H2 molecule,
  5. Comparison between the LOE and XLOE results.

Device setup

Start up VNL and create a new project. Instead of using builder_icon Builder to create the device configuration, we use the ATK Python script, which provides the Au|H2|Au device. Download and save the script in the project folder, and open it from script_generator_icon Script Generator.

  • Add a calculator_icon New Calculator block and select GGA as the exchange correlation functional in the setting.


The k-point sampling is by default 1x1x93, which should, however, be modified properly if the device is periodic in the A and/or B directions.

  • Set the output-file name “”, and save the python script as “”.


The python script must now be the same as To save time, one may just download the script and use it.

  • Send the script to job_manager_icon Job Manager and run it.

Calculation of IETS

The IETS calculation requires information of the dynamical matrix, which is defined in terms of the force for displacing the hydrogen atoms back and forth in each spatial direction. The computation involves six DFT calculations, which are all independent with each other, and is hence performed efficiently by using six CPU cores.

  • In the script_generator_icon Script Generator, click on analysis_from_file_icon Analysis from File and select “”.
  • Add a analysis_icon Analysis ‣ DynamicalMatrix object and set it as follows:
  • Add also a analysis_icon Analysis ‣ VibrationalModes object and keep the default settings.

The IETS calculation requires information of the Hamiltonian derivatives as well. This calculation also consists of six independent parts, and is hence carried out efficiently using six CPU cores.

  • Add a analysis_icon Analysis ‣ HamiltonianDerivatives block and set it as follows:
  • Add an analysis_icon Analysis ‣ InelasticTransmissionSpectrum, choose LOE as the method, and set the rest as follows:
  • Add another analysis_icon Analysis ‣ TransmissionEigenstates block, and set it in the same way except choosing XLOE as the method.
  • Set the output-file name “”, and save the python script as “”.


The python script must be the same as To save time, one may just download the script and use it.

  • Send the script to job_manager_icon Job Manager and run it.

The calculation will take less than one hour if six CPU cores are used, and save every result in “”.


See InelasticTransmissionSpectrum in ATK Manual, which provides a brief description about IETS and how it is computed from the dynamical matrix and the Hamiltonian derivatives. See also, e.g., Ref. [FPBJ07] for more details.


Vibrational modes of the H2 molecule

In the Labfloor, look at the vibrational modes of the H2 molecule by selecting the labfloor_vibrationalmode_icon Vibrational Mode object in the “” file with the help of Vibration Visualizer plugin.


There are, in total, six modes: Four transverse modes 0, 1, 2, 3; two longitudinal modes 4 and 5. The two pairs of modes (0 and 1) and (2 and 3) are, respectively, doubly degenerated. Note that the temperature is set T = 1000 K in the movies below to see the vibrations clearly.


Mode 0: 25.44 meV; Mode 1: 25.63 meV (the small energy difference is due to numerical error)


Mode 2: 60.54 meV; Mode 3: 60.56 meV (the small energy difference is due to numerical error)


Mode 4: 127.46 meV


Mode 5: 264.06 meV


Then look at IETS. Select one of the labfloor_inelastictransmissionspectrum_icon InelasticTransmissionSpectrum objects \(T({\bf k},{\bf q})\) in LabFloor and click on the Inelastic Transmission Spectrum Analyzer plugin, a new plugin available from ATK2017.


Set the panel as follows and click Show.



Because the density of states (DOS) of the gold electrodes is almost flat within this bias window, there only exists a minor difference between the LOE and XLOE results. See the discussion below for more details.


There appear three significant inelastic features (steps in conductance and peaks in IETS). Note that the inelastic conductance is increased by the coupling to the low-frequency transverse modes 0-3 and decreased by the coupling to the high-frequency longitudinal mode 5. Analyzing the scattering states clarifies the reason of the increases.

  • Drag the labfloor_transmissioneigenstate_icon TransmissionEigenstates object \(\psi_{i,{\bf k},\epsilon}\) in LabFloor into viewer_icon Viewer. Then drag the twoprobe_icon DeviceConfiguration into viewer_icon Viewer.
  • Go to properties and, under the Isosurface flag, set the Isovalue to 0.5.

The scattering eigenstate has a clear d-symmetry on the gold chain. Hence, it couples only very weakly to the s-like states of the hydrogen atoms. The elastic transmission is thus small. The coupling to the transverse modes 0-3 can, however, open new inelastic channels and leads to positive steps in conductance and steps in IETS.

Detailed analysis of transmission: Comparison between LOE and XLOE

In the same window of the Inelastic Transmission Spectrum Analyzer plugin, now choose Transmissions in Plot type, set the options as shown below, and click Show. In this case, one will look at the positive symmetric and asymmetric parts of the inelastic transmission:


To understand some details of the transmission, note that, in either LOE or XLOE, the current as a function of the bias voltage \(V\) and the temperature \(T\) reads

\[I(V,T) = \sum_{\lambda} \Big[ \mathcal{I}^{sym}_{\lambda}(V,T) \mathcal{T}^{sym}_{\lambda}(\epsilon_F) + \mathcal{I}^{asym}_{\lambda}(V,T) \mathcal{T}^{asym}_{\lambda}(\epsilon_F) \Big],\]

where \(\epsilon_F\) is the Fermi energy, and the rest of the new symbols are defined as follows: In either approximation, the universal functions appear in the same form,

\[\begin{split}\mathcal{I}^{sym}_{\lambda}(V,T) &=& \frac{G_0}{2e}\sum_{\sigma=\pm}\sigma(\hbar\omega_{\lambda}+\sigma eV) \left( \coth\frac{\hbar\omega_{\lambda}}{2k_BT} - \coth\frac{\hbar\omega_{\lambda}+\sigma eV}{2k_BT} \right), \\ \mathcal{I}^{asym}_{\lambda}(V,T) &=& \frac{G_0}{2e}\int_{-\infty}^{\infty}d\epsilon \big[ n_F(\epsilon-eV)-n_F(\epsilon) \big] \cdot \mathcal{H}_{\epsilon'}\big[ n_F(\epsilon'-\omega_{\lambda})-n_F(\epsilon'+\omega_{\lambda}) \big](\epsilon),\end{split}\]

where \(G_0=2e^2/h\) is the conductance quantum, \(n_F(\cdot)\) is the Fermi-Dirac distribution function, and \(\mathcal{H}_{\epsilon'}[\cdot](\epsilon)\) means the Hilbert transform. In XLOE, the transmission functions are given by

\[\begin{split}\mathcal{T}^{sym}_{\lambda}(\epsilon) &=& {\rm Tr}\big[ M_{\lambda}\tilde{A}_L(\mu_L)M_{\lambda}A_R(\mu_R) \big] + {\rm Im}\:B_{\lambda}(\epsilon), \\ \mathcal{T}^{asym}_{\lambda}(\epsilon) &=& 2{\rm Re}\:B_{\lambda}(\epsilon),\end{split}\]


\[\begin{split}\mu_L&=&\epsilon, \\ \mu_R&=&\epsilon\pm \hbar\omega_{\lambda},\end{split}\]


\[B_{\lambda}(\epsilon) = {\rm Tr}\big[ M_{\lambda}A_R(\mu_L)\Gamma_L(\mu_L)G^r(\mu_L)M_{\lambda}A_R(\mu_R) -M_{\lambda}G^a(\mu_R)\Gamma_L(\mu_R)A_R(\mu_R)M_{\lambda}A_L(\mu_L) \big].\]

This XLOE expression of the transmission functions reduces to the LOE expression by setting \(\mu_L=\mu_R=\epsilon\). Importantly, in the expression of the current, note that the Green’s function, the self-energies, and the spectral functions in the transmission functions are evaluated at two different energies \(\epsilon_F\pm \hbar\omega_{\lambda}\) in XLOE, whereas they are all evaluated at the same energy \(\epsilon_F\) in LOE. See Ref. [LCF+14] for more details.

For direct comparison between the LOE and XLOE results, see here below the inelastic transmissions for each phonon mode. The lines of positive/negative (a)symmetric transmissions indicate the plots of \(\mathcal{T}^{(a)sym}_{\lambda}(\epsilon)\) as functions \(\epsilon\) with setting \(\mu_R=\epsilon\pm \hbar\omega_{\lambda}\). The positive and negative (a)symmetric parts are exactly the same in LOE by definition, but different in XLOE.


The LOE and XLOE results are in good agreement in areas where the transmission varies slowly as a function of \(\epsilon\). Around -0.11 eV, however, the wide-band approximation in LOE, i.e., the approximation \(\mu_L=\mu_R=\epsilon\), breaks down due to the presence of a resonance in the DOS of the 1D gold chain. As a result, there appears a difference between the LOE and XLOE results. The symmetric and asymmetric inelastic transmissions for mode 5 show the difference very clearly. The single peak at -0.11 eV in LOE splits into two peaks in XLOE; the separation is about 264 meV, which corresponds to the phonon energy of mode 5.


[FPBJ07](1, 2) T. Frederiksen, M. Paulsson, M. Brandbyge, and A.-P. Jauho. Inelastic transport theory from first principles: Methodology and application to nanoscale devices. Phys. Rev. B, 75:205413, May 2007. doi:10.1103/PhysRevB.75.205413.
[LCF+14](1, 2) J.-T. Lü, R. B. Christensen, G. Foti, T. Frederiksen, T. Gunst, and M. Brandbyge. Efficient calculation of inelastic vibration signals in electron transport: Beyond the wide-band approximation. Phys. Rev. B, 89:081405, Feb 2014. doi:10.1103/PhysRevB.89.081405.
[Ree08]M. A. Reed. Inelastic electron tunnelling spectroscopy. Materials Today, 2008. doi:10.1016/S1369-7021(08)70238-4.