Inelastic current in a silicon pn junction¶
Version: 2017.2
Inelastic scattering effects due to electronphonon interactions can play a fundamental role in determining the carrier transport through a device. In this tutorial, you will investigate the impact of such effects on the current through a silicon \(p\)\(n\) junction in reverse bias. In this system, carrier transport is dominated by bandtoband tunnelling between the valence band of the \(p\) region and the conduction band of the \(n\) region. This process can be enhanced considerably by inelastic scattering effects.
Specifically, you will use the InelasticTransmissionSpectrum module implemented in ATK to calculate the transmission spectrum of the device in the presence of electronphonon interactions. You will also use the Inelastic Transmission Spectrum Analyzer plugin in VNL to perform a detailed analysis of the current and determine which phonon mode is the main responsible for the current enhancement.
The methodology implemented in ATK is based on the Lowest Order Expansion (LOE) [FPBJ07] and the eXtended LOE (XLOE) [LCF+14] methods. In this tutorial, you will use the XLOE method, which treats the finite energy difference between initial and final states in the transition \(\mathbf{k} \to \mathbf{k}\pm\mathbf{q}\) due to electronphonon coupling in a more accurate manner. However, the tutorial can be carried out also using the simpler LOE method, which will give qualitatively similar results. You can find more details about the theory underlying the two methods in the Notes sections of the InelasticTransmissionSpectrum module in the ATK manual.
Creating the silicon pn junction¶
To create the device, you can follow the instruction in the Create device section of the Silicon pn junction tutorial. However, in the present tutorial you will create a shorter device to speed up the calculations. In the Builder, contruct the device with the following changes:
When using the Surface (Cleave) plugin to create a <100>oriented structure, set the thickness to 24 layers instead of 52 layers;
When using the Device from Bulk plugin, choose 5.4306 Å as the electrode length;
Click the Doping plugin in \(p\) and \(n\) regions. In this case, select the half of device to set the \(p\) doping region and select the remainder part using the
to set the doping of thectrl+i
which will select the inversion part to set the \(n\) doping region. In both regions, set a doping of \(2 \cdot 10^{19} \mathrm{e/cm^{3}}\).
You can download the resulting device configuration from here: Si_pn_junction.py
.
Transmission calculation without electronphonon interactions¶
In this section, you will calculate and analyze the electronic structure of the \(p\)\(n\) junction at a reverse bias voltage of \(V_\mathrm{bias} = 0.4\ \mathrm{V}\), and the associated transmission spectrum without electronphonon interactions included, using the conventional Landauer formula as implemented in ATK in the TransmissionSpectrum module.
Setting up the calculation¶
From the Builder, send the structure to the Script generator using the button. Add the following blocks and set their parameters as follows:

Set the Calculator to ‘ATKSE: SlaterKoster (Device)‘
In Basic device:
 untick the ‘No SCF iteration‘ option
 set the Left electrode voltage to \(0.2\ \mathrm{V}\)
 set the Right electrode voltage to \(0.2\ \mathrm{V}\)
In Numerical accuracy parameters, set the kpoint Sampling to:
 \(\mathrm{k_A} = 7\)
 \(\mathrm{k_B} = 7\)
 \(\mathrm{k_C} = 101\)
In Iteration control parameters, set the Tolerance to ‘4e05‘

Set the kpoint Sampling to:
 \(\mathrm{k_A} = 21\)
 \(\mathrm{k_B} = 21\)

Set the kpoint Sampling to:
 \(\mathrm{k_A} = 21\)
 \(\mathrm{k_B} = 21\)
Finally, in the Global IO options, change the name of the Default output file to
transmission_V0.4.hdf5
.
Send the script to the Job manager using the button,
save it as transmission_V0.4.py
and press the button to run the
calculation. It will take less than 2 minutes on 4 CPUs. You can also download the full script
from here: transmission_V0.4.py
.
Analysis of the results¶
In the LabFloor, select the
ProjectedLocalDensityOfStates object from the file transmission_V0.4.hdf5
, and
click on the Projected Local Density of States plugin. In the plugin window, select
Spectral Current from the dropdown menu on the top left, and set the maximum value of the
Data range to 0.1. In the panel on the lefthand side of the screen, use the Zoom button
to make sure that the minimum value of the current is about \(10^{24} \mathrm{A/eV}\). The features
below this value can be associated to numerical noise.
From the left panel of the figure above, it can be seen that the maximum of the spectral current occurs inside the bias window (\(0.2\ \mathrm{eV} \leq \mathrm{Energy} \leq 0.2\ \mathrm{eV}\)). Furthermore, from the right panel, which shows the density of states across the device, it can be seen that within the bias window, there are no electronic states in central region (\(20\ \mathrm{Å} \leq \mathrm{z} \leq 45\ \mathrm{Å}\)) of the device, so that tunnelling between the left and the right electrodes is the only possible transport mechanism. This indicates that electronic transport in the device in reverse bias is dominated by tunnelling.
To analyze the probability of tunnelling in \(\mathbf{k}\)space, select the Transmission Spectrum object from the same file and click on the Transmission Analyzer plugin. In the plugin window, set the maximum value of the Data range to 0.1.
From the figure above, it can be seen that the tunnelling probability is strongly peaked at the \(\Gamma\)point of the Brilluoin zone at the Fermi level.
Transmission calculation with electronphonon interactions¶
In this section, you will calculate the transmission spectrum of the silicon \(p\)\(n\) junction at a reverse bias voltage of \(V_\mathrm{bias} = 0.4\ \mathrm{V}\) including the effect of the electronphonon interactions within the XLOE approximation [LCF+14].
To calculate the InelasticTransmissionSpectrum, you first need to calculate the DynamicalMatrix and the HamiltonianDerivatives of the system.
Setting up the dynamical matrix calculation¶
In the LabFloor, drag and drop the DeviceConfiguration object contained in
the file transmission_V0.4.hdf5
to the Script generator,
and modify the New Calculator block as follows:
 Select ATKForceField from the list of available Calculators
 Select ‘StillingerWeber_Si_1985‘ from the list of available classical potentials in Potential settings
Add the following blocks and set their parameters as follows:

Set Repetitions to ‘Custom‘
Set the number of repetitions to:
 \(\mathrm{n_A} = 3\)
 \(\mathrm{n_B} = 3\)
 \(\mathrm{n_C} = 1\)
Untick the ‘Acooustic sum rule‘ option
Tick the ‘Constrain electrodes‘ option
Untick the ‘Equivalent bulk‘ option
Finally, in the Global IO options, change the name of the Default output file to
dynmat.hdf5
.
Send the script to the Job manager using the button,
save it as dynmat.py
and press the button to run the calculation. You can
also download the full script from here: dynmat.py
.
Note
Alternatively, you can also use the special thermal displacement  Landauer method (STDLandauer) to include part of the electronphonon coupling effects in the transmission calculation at a much lower computational cost. The STDLandauer case study treats a system which is similar to that discussed in this tutorial.
Setting up the Hamiltonian derivatives calculation¶
In the LabFloor, drag and drop the DeviceConfiguration object contained in
the file transmission_V0.4.hdf5
to the Script generator,
and modify the New Calculator block as follows:
 In Basic device tick the ‘No SCF iteration‘ option.
Add an block and set the parameters as follows:
 Set Repetitions to ‘Custom‘
 Set the number of repetitions to:
 \(\mathrm{n_A} = 3\)
 \(\mathrm{n_B} = 3\)
 \(\mathrm{n_C} = 1\)
 In Constrains, click Add and constrain the electrode repetitions with Fixed constrains. The resulting list of constrained atoms should be: 0,1,2,3,44,45,46,47.
 Untick the ‘Equivalent Bulk‘ option
In the Global IO options, change the name of the Default output file to
dHdR_V0.4.hdf5
.
Send the script to the Job manager using the button,
save it as dHdR_V0.4.py
and press the button to run the
calculation. The calculation will take around 20 minutes on 8 CPUs. You can also download
the full script from here: dHdR_V0.4.py
.
Note
Here, the Hamiltonian derivatives are calculated nonselfconsistently to speed up the calculations.
Calculating the inelastic transmission spectrum¶
To calculate the inelastic part of the current, click the Script generator, and add the following blocks:
Open the block, select the
transmission_V0.4.hdf5
file and load the DeviceConfiguration object
contained in the file .
You will notice that two additional blocks have also been added:
In our case, the dynamical matrix and the Hamiltonian derivatives have already been calculated previously and can be reused. Delete the two blocks, add two more blocks and arrange them as shown in the figure below:
Click on the second block from the top.
Select the file dynmat.hdf5
and load the DynamicalMatrix object. Uncheck the
DeviceConfiguration object included in the same file. The panel should look like as in
the figure below:
Then, click on the third block from the top. Select
the file dHdR_V0.4.hdf5
and load the HamiltonianDerivatives object. The panel should
look like as in the figure below:
Finally, set the parameters for the block as shown below:
In the Energy part, set:
 \(\mathrm{E_0} = 0.5\ \mathrm{eV}\)
 \(\mathrm{E_1} = 0.5\ \mathrm{eV}\)
 Points = 25
In the kpoint Sampling part:
set Grid type to ‘Regular kpoint grid‘
set the grid range for \(\mathrm{k_A}\) and \(\mathrm{k_B}\) to \([0.2 : 0.2]\)
set the Sampling to:
 \(\mathrm{k_A} = 3\)
 \(\mathrm{k_B} = 3\)
In the qpoint Sampling part:
set Grid type to ‘Regular qpoint grid‘
set the grid range for \(\mathrm{q_A}\) and \(\mathrm{q_B}\) to \([0.5 : 0.5]\)
set the Sampling to:
 \(\mathrm{q_A} = 9\)
 \(\mathrm{q_B} = 9\)
In the Global IO options, change the name of the Default output file to
xloe_V0.4.hdf5
.
Send the script to the Job manager using the button,
save it as xloe_V0.4.py
and press the button to run the
calculation. The calculation will take around 1.5 hour on 24 CPUs.
You can also download the full script from here: xloe_V0.4.py
.
Analysis of the inelastic transmission spectrum¶
The results of the inelastic transmission spectrum can be analyzed in detail using the Inelastic Transmission Spectrum Analyzer.
In the LabFloor, select the
InelasticTransmissionSpectrum object contained in the file xloe_V0.4.hdf5
,
and click on the Inelastic Transmission Spectrum Analyzer plugin in the plugins panel
on the righthand side of the screen.
First of all, you will analyze the \(\mathbf{k}\) and \(\mathbf{q}\)dependency of the current. In the main window of the analyzer, set the following parameters:
 Plot type : ‘Current vs. kpoints‘
 Bias voltage : \(0.4\ \mathrm{V}\)
You will obtain a figure similar to the one below, showing the \(\mathbf{k}\)dependency of the current within the sampled range of kpoints. From the figure, it is evident that the main contribution to the current comes from the \(\Gamma\)point of the twodimensional Brillouin zone defined by \(k_\mathrm{A}\) and \(k_\mathrm{B}\), with a nonnegligible contribution from finite \(\mathbf{k}\)vectors.
Next, change the Plot type to ‘Current vs. qpoints‘. You will obtain a figure similar to the one below, showing the \(\mathbf{q}\)dependency of the current. From the figure, it is evident that also in this case the main contribution to the current comes from the \(\Gamma\)point of the twodimensional Brillouin zone defined by \(q_\mathrm{A}\) and \(q_\mathrm{B}\), , with a nonnegligible contribution from finite \(\mathbf{q}\)vectors.
The two figures above indicate that the current is mainly associated with transitions occurring at the \(\Gamma\)point of the twodimensional Brillouin zone, both in \(\mathbf{k}\)space and in \(\mathbf{q}\)space. To analyze which of the phonon modes at the \(\Gamma\)point is the one that contributes the most, set the following parameters in the main window of the analyzer:
Plot type : ‘Current vs. phonon mode‘
kpoints : ‘(0.00, 0.00)‘
qpoints : ‘(0.00, 0.00)‘
From the figure, it can be seen that the phonon that contributes the most to the current is that with index \(66\), and energy \(\hbar \omega = 63.01\ \mathrm{meV}\).
This phonon mode can be visualized by selecting the VibrationalMode object in the file
dynmat.hdf5
calculated in Section Setting up the dynamical matrix calculation and using the
Vibration Visualizer plugin:
Note
In the figures presented above, the values of the current are negative, because the simulations are performed in reverse bias conditions.
Calculation of the current with and without electronphonon interactions¶
In order to calculate the current without and with electronphonon interactions, download the script
current.py
and run it from the terminal as:
atkpython current.py transmission_V0.4.hdf5 xloe_V0.4.hdf5
.
Note
On windows, it is not convenient to run current.py
script from the terminal. In this script, you can replace filename1 and filename2 to the transmission_V0.4.hdf5
and xloeV0.4.hdf5
. In order to run it, drag it on the Job manager.
The script will calculate the elastic (\(I_\mathrm{el}\)) and inelastic (\(I_\mathrm{inel}\)) components of the current. The total current without and with electronphonon interactions will then be calculated as:
The output of the calculation will show the name of the files used as input, as well as the values of the calculated currents (highlighted in yellow):
++
 
 Atomistix ToolKit 2017.2 [Build 997d195] 
 
++

TransmissionSpectrum read from file :
transmission_V0.4.hdf5
InelasticTransmissionSpectrum read from file :
xloe_V0.4.hdf5

Current (w/o interactions) = 3.71e09 nA
Current (w interactions) = 9.03e06 nA

Timing: Total Per Step %

Loading Modules + MPI : 1.33 s 1.33 s 2.58% =============

Total : 51.58 s
Note
The above calculated currents are different with the currents before analysis of the InelasticTransmissionSpectrum because of different kpoints samplings. Also in this case, the value of the current is negative, because the simulations are performed in reverse bias conditions.
It can be seen that electronphonon scattering leads to an increase in the reverse bias current of about three orders of magnitude!
Speeding up the calculations¶
ATK implements several methods which can be used to speed up the calculations considerably and enable the calculation of the inelastic transmission for devices comprised of thousands of atoms. These methods are:
 The ‘Phonon energy intervals’ method
 The possibility of using the bulk dynamical matrix and Hamiltonian derivatives
For a more detailed description of these methods, see the section Large Device Calculations of the InelasticTransmissionSpectrum module in the ATK manual.
The ‘Phonon energy intervals’ method¶
This method allows us to group \(3N\) phonon modes of a device with \(N\) vibrating atoms into \(M\) energy intervals to form \(M\) new effective phonon modes, with \(M << 3N\). The inelastic transmission spectrum will therefore be calculated only for these \(M\) effective phonon modes, greatly reducing the computational cost of the calculation.
In the LabFloor, drag and drop the DeviceConfiguration object contained in
the file transmission_V0.4.hdf5
to the Script generator.
Setup the calculation of the inelastic transmission spectrum as in the Calculating the inelastic transmission spectrum
Section, but in this case set the Phonon Method in the Phonons part to ‘Phonon energy intervals‘.
Note
The default parameters of the Phonon energy intervals are fine for the present calculation, but in general one should ensure that the energy range \([E_\mathrm{0}:E_\mathrm{1}]\) spans that of the phonon modes of the calculated dynamical matrix.
In the Global IO options, change the name of the Default output file to
xloe_pheint_V0.4.hdf5
.
Send the script to the Job manager using the button,
save it as xloe_pheint_V0.4.py
and press the button to run the
calculation. The calculation will take only 20 minutes on 24 CPUs.
You can also download the full script from here: xloe_pheint_V0.4.py
.
Once the calculation is finished, run the script current.py
as:
atkpython current.py transmission_V0.4.hdf5 xloe_pheint_V0.4.hdf5
.
If you are a Windows user, see how to run the script.
The resulting output will be:
++
 
 Atomistix ToolKit 2017.2 [Build 997d195] 
 
++

TransmissionSpectrum read from file :
transmission_V0.4.hdf5
InelasticTransmissionSpectrum read from file :
xloe_pheint_V0.4.hdf5

Current (w/o interactions) = 3.71e09 nA
Current (w interactions) = 9.25e06 nA

Timing: Total Per Step %

Loading Modules + MPI : 1.34 s 1.34 s 19.27% =============

Total : 6.94 s
One can see that the calculated values for the current are essentially the same as those calculated in the Section Calculation of the current with and without electronphonon interactions, despite the fact that the present calculations are much faster.
Using the bulk dynamical matrix and Hamiltonian derivatives¶
Another option to speed up the calculations is to use the dynamical matrix and Hamiltonian derivatives of the bulk electrodes, instead of those of the device. However, this is possible only for devices with a structure which is translationally invariant along the \(\mathrm{C}\)direction, apart from doping, electrostatic regions and applied bias voltage.
The scripts needed to calculate the dynamical matrix and the Hamiltonian
derivatives of the bulk electrodes can be downloaded from here: dynmat_bulk.py
, dHdR_bulk.py
.
Running the two calculations on 4 CPUs will take less than 2 minutes.
In order to calculate the current inelastic transmission spectrum using the bulk dynamical matrix and Hamiltonian derivatives, repeat the steps followed in Section Calculating the inelastic transmission spectrum, with the changes described in the following.
First, click the Script generator, add the following blocks, and modify their parameters as follows:

 Load the DeviceConfiguration object from the
tranmission_V0.4.hdf5
file.
 Load the DeviceConfiguration object from the

 Set parameters as described in the Inelastic transmission spectrum section.
Remove the and blocks, add two more blocks right after the already present, and set their properties as follows:
 Click on the second block from the top, select the
dynamic_bulk.hdf5
file and load the DynamicalMatrix object contained in the file.  Click on the second block from the top, select the
dHdR_bulk.hdf5
file and load the HamiltonianDerivatives object contained in the file.
In the Global IO options, set the name of the Default output file to xloe_V0.4_bulk.hdf5
.
In this case, you need to modify the Inelastic Transmission Spectrum block in the Editor as follows:
1 2 3 4 5 6 7 8 9 10 11 12 13 14 inelastic_transmission_spectrum = InelasticTransmissionSpectrum( configuration=device_configuration, bulk_dynamical_matrix=dynamical_matrix, bulk_hamiltonian_derivatives=hamiltonian_derivatives, energies=numpy.linspace(0.5, 0.5, 25)*eV, kpoints=kpoints, qpoints=qpoints, self_energy_calculator=RecursionSelfEnergy(), energy_zero_parameter=AverageFermiLevel, infinitesimal=1e06*eV, phonon_modes=All, method=XLOE, spectral_representation=True, electrode_extensions=[0, 0],
Send the script to the Job manager using the button,
save it as xloe_V0.4_bulk.py
and press the button to run the
calculation. Also in this case, the calculation will take only 20 minutes on 24 CPUs.
You can also download the full script from here: xloe_V0.4_bulk.py
.
Once the calculation is finished, run the script current.py
as:
atkpython current.py transmission_V0.4.hdf5 xloe_V0.4_bulk.hdf5
.
If you use Windows, look at the how to run the script.
The resulting output will be:
++
 
 Atomistix ToolKit 2017.2 [Build 997d195] 
 
++

TransmissionSpectrum read from file :
transmission_V0.4.hdf5
InelasticTransmissionSpectrum read from file :
xloe_V0.4_bulk.hdf5

Current (w/o interactions) = 3.71e09 nA
Current (w interactions) = 6.42e06 nA

Timing: Total Per Step %

Loading Modules + MPI : 1.33 s 1.33 s 2.17% =============

Total : 61.15 s (1m01.15s)
It can be seen that the calculated values for the current is similar to that calculated using the dynamical matrix and Hamiltonian derivatives of the device configuration, although there are differences due to the fact that in the present case, the loss translational invariance of the Hamiltonian derivatives due to the applied bias is not taken into account.
References¶
[FPBJ07]  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 wideband approximation. Phys. Rev. B, 89:081405, Feb 2014. doi:10.1103/PhysRevB.89.081405. 