Inelastic current in a silicon p-n junction

Version: 2017.0

Inelastic scattering effects due to electron-phonon 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 band-to-band 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 electron-phonon 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 electron-phonon 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 p-n junction

To create the device, you can follow the instruction in the Create device section of the Silicon_p_n tutorial. However, in the present tutorial you will create a shorter device to speed up the calculations. In the builder_icon 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;

  • Use the Doping plugin in Miscellaneous ‣ Doping to set the doping of the \(p\) and \(n\) regions to a value of \(2 \cdot 10^{19} \mathrm{cm^{-3}}\).

    ../../_images/doping.png

You can download the resulting device configuration from here: Si_pn_junction.py.

Transmission calculation without electron-phonon 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 electron-phonon interactions included, using the conventional Landauer formula as implemented in ATK in the TransmissionSpectrum module.

Setting up the calculation

From the builder_icon Builder, send the structure to the script_generator_icon Script generator using the sendto_icon button. Add the following blocks and set their parameters as follows:

  • calculator_icon New Calculator

    • Set the Calculator to ‘ATK-SE: Slater-Koster (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 k-point Sampling to:

      • \(\mathrm{k_A} = 7\)
      • \(\mathrm{k_B} = 7\)
      • \(\mathrm{k_C} = 100\)
    • In Iteration control parameters, set the Tolerance to ‘4e-05

  • analysis_icon Analysis ‣ ProjectedLocalDensityOfStates

    • Set the k-point Sampling to:

      • \(\mathrm{k_A} = 21\)
      • \(\mathrm{k_B} = 21\)
  • analysis_icon Analysis ‣ TransmissionSpectrum

    • Set the k-point 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_V-0.4.hdf5‘.

Send the script to the job_manager_icon Job manager using the sendto_icon button, save it as ‘transmission_V-0.4.py‘ and press the start_icon 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_V-0.4.py.

Analysis of the results

In the LabFloor, select the labfloor_projectedlocaldensityofstates_icon ProjectedLocalDensityOfStates object from the file ‘transmission_V-0.4.hdf5‘, and click on the Projected Local Density of States plugin. In the plugin window, select Spectral Current from the drop-down menu on the top left, and set the maximum value of the Data range to 0.1. The resulting plot should look like the figure below.

../../_images/pldos_spectralcurrent.png

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

../../_images/transmission_analyzer.png

From the figure above, it can be seen that the tunnelling probability is strongly peaked at the \(\Gamma\)-point of the Brilluoin zone.

Transmission calculation with electron-phonon 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 electron-phonon 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_V-0.4.hdf5‘ to the script_generator_icon Script generator, and modify the calculator_icon New Calculator block as follows:

  • Select ATK-ForceField 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:

  • analysis_icon Analysis ‣ DynamicalMatrix

    • 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

  • analysis_icon Analysis ‣ VibrationalMode

Finally, in the Global IO options, change the name of the Default output file to ‘dynmat.hdf5‘.

Send the script to the job_manager_icon Job manager using the sendto_icon button, save it as ‘dynmat.py‘ and press the start_icon button to run the calculation. You can also download the full script from here: dynmat.py.

Setting up the Hamiltonian derivatives calculation

In the LabFloor, drag and drop the DeviceConfiguration object contained in the file ‘transmission_V-0.4.hdf5‘ to the script_generator_icon Script generator, and modify the calculator_icon New Calculator block as follows:

  • In Basic device tick the ‘No SCF iteration‘ option.

Add an analysis_icon Analysis ‣ HamiltonianDerivatives 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_V-0.4.hdf5‘.

Send the script to the job_manager_icon Job manager using the sendto_icon button, save it as ‘dHdR_V-0.4.py‘ and press the start_icon 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_V-0.4.py.

Note

Here, the Hamiltonian derivatives are calculated non-selfconsistently to speed up the calculations.

Calculating the inelastic transmission spectrum

In the LabFloor, drag and drop the DeviceConfiguration object contained in the file ‘transmission_V-0.4.hdf5‘ to the script_generator_icon Script generator, add the following blocks and modify their parameters as follows:

  • analysis_icon Analysis ‣ InelasticTransmissionSpectrum

    • In the Energy part, set:

      • \(\mathrm{E_0} = -0.5\ \mathrm{eV}\)
      • \(\mathrm{E_1} = 0.5\ \mathrm{eV}\)
      • Points = 25
    • In the k-point Sampling part:

      • set Grid type to ‘Regular k-point 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 q-point Sampling part:

      • set Grid type to ‘Regular q-point 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\)
    ../../_images/inelastictransmissionsetup_phall.png

In the Global IO options, change the name of the Default output file to ‘xloe_V-0.4.hdf5‘.

Send the script to the editor_icon Editor using the sendto_icon button, and modify the Dynamical Matrix and Hamiltonian Derivatives blocks as follows:

1
2
3
4
5
6
7
8
9
# -------------------------------------------------------------
# Dynamical Matrix
# -------------------------------------------------------------
dynamical_matrix = nlread('dynmat.hdf5',DynamicalMatrix)[-1]

# -------------------------------------------------------------
# Hamiltonian Derivatives
# -------------------------------------------------------------
hamiltonian_derivatives = nlread('dHdR_V-0.4.hdf5',HamiltonianDerivatives)[-1]

Send the script to the job_manager_icon Job manager using the sendto_icon button, save it as ‘xloe_V-0.4.py‘ and press the start_icon 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_V-0.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 labfloor_inelastictransmissionspectrum_icon InelasticTransmissionSpectrum object contained in the file ‘xloe_V-0.4.hdf5‘, and click on the Inelastic Transmission Spectrum Analyzer plugin in the plugins panel on the right-hand 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. k-points
  • 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 k-points. From the figure, it is evident that the main contribution to the current comes from the \(\Gamma\)-point of the two-dimensional Brillouin zone defined by \(k_\mathrm{A}\) and \(k_\mathrm{B}\), with a non-negligible contribution from finite \(\mathbf{k}\)-vectors.

../../_images/current_vs_kpoints.png

Next, change the Plot type to ‘Current vs. q-points‘. 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 two-dimensional Brillouin zone defined by \(q_\mathrm{A}\) and \(q_\mathrm{B}\), , with a non-negligible contribution from finite \(\mathbf{q}\)-vectors.

../../_images/current_vs_qpoints.png

The two figures above indicate that the current is mainly associated with transitions occurring at the \(\Gamma\)-point of the two-dimensional 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

  • k-points : ‘(0.00, 0.00)

  • q-points : ‘(0.00, 0.00)

    ../../_images/current_vs_phononmode.png

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:

../../_images/vibration_visualizer.png
../../_images/mode66.gif

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 electron-phonon interactions

In order to calculate the current without and with electron-phonon interactions, download the script current.py and run it from the terminal as:

atkpython current.py transmission_V-0.4.hdf5 xloe_V-0.4.hdf5.

The script will calculate the elastic (\(I_\mathrm{el}\)) and inelastic (\(I_\mathrm{inel}\)) components of the current. The total current without and with electron-phonon interactions will then be calculated as:

\[I (\mathrm{w\ interactions}) = I_\mathrm{el}\]\[I (\mathrm{w/o\ interactions}) = I_\mathrm{el} + I_\mathrm{inel}\]

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.0 [Build 15d690e]                                     |
 |                                                                              |
 +------------------------------------------------------------------------------+
 --------------------------------------------------------------------------------
  TransmissionSpectrum read from file :
  transmission_V-0.4.hdf5
  InelasticTransmissionSpectrum read from file :
  xloe_V-0.4.hdf5
 --------------------------------------------------------------------------------
  Current (w/o interactions)    = -3.71e-09 nA
  Current (w interactions)      = -9.03e-06 nA
 --------------------------------------------------------------------------------

 Timing:                          Total     Per Step        %

 --------------------------------------------------------------------------------

 Loading Modules + MPI   :       1.33 s       1.33 s       2.58% |=============|
 --------------------------------------------------------------------------------
 Total                   :      51.58 s

It can be seen that electron-phonon scattering leads to an increase in the reverse bias current of about three orders of magnitude!

Note

Also in this case, the value of the current is negative, because the simulations are performed in reverse bias conditions.

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:

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_V-0.4.hdf5‘ to the script_generator_icon 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‘.

../../_images/inelastictransmissionsetup_phenergyintervals.png

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_V-0.4.hdf5‘.

Send the script to the job_manager_icon Job manager using the sendto_icon button, save it as ‘xloe_pheint_V-0.4.py‘ and press the start_icon 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_V-0.4.py.

Once the calculation is finished, run the script current.py as:

atkpython current.py transmission_V-0.4.hdf5 xloe_pheint_V-0.4.hdf5.

The resulting output will be:

 +------------------------------------------------------------------------------+
 |                                                                              |
 | Atomistix ToolKit 2017.0 [Build 15d690e]                                     |
 |                                                                              |
 +------------------------------------------------------------------------------+
 --------------------------------------------------------------------------------
  TransmissionSpectrum read from file :
  transmission_V-0.4.hdf5
  InelasticTransmissionSpectrum read from file :
  xloe_pheint_V-0.4.hdf5
 --------------------------------------------------------------------------------
  Current (w/o interactions)    = -3.71e-09 nA
  Current (w interactions)      = -9.25e-06 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 electron-phonon 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 using the bulk dynamical matrix and Hamiltonian derivatives, repeat the steps followed in Section Calculating the inelastic transmission spectrum, with the following changes:

  • In the Global IO options, set the name of the Default output file to ‘xloe_V-0.4_bulk.hdf5

  • In the editor_icon Editor, modify the Dynamical Matrix and Hamiltonian Derivatives blocks as follows:

    1
    2
    3
    4
    5
    6
    7
    8
    9
    # -------------------------------------------------------------
    # Dynamical Matrix
    # -------------------------------------------------------------
    dynamical_matrix = nlread('dynmat_bulk.hdf5',DynamicalMatrix)[-1]
    
    # -------------------------------------------------------------
    # Hamiltonian Derivatives
    # -------------------------------------------------------------
    hamiltonian_derivatives = nlread('dHdR_bulk.hdf5',HamiltonianDerivatives)[-1]
    
  • In the editor_icon Editor, modify the Inelastic Transmission Spectrum block 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=1e-06*eV,
        phonon_modes=All,
        method=XLOE,
        spectral_representation=True,
        electrode_extensions=[0, 0],
    

Send the script to the job_manager_icon Job manager using the sendto_icon button, save it as ‘xloe_V-0.4_bulk.py‘ and press the start_icon 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_V-0.4_bulk.py.

Once the calculation is finished, run the script current.py as:

atkpython current.py transmission_V-0.4.hdf5 xloe_V-0.4_bulk.hdf5.

The resulting output will be:

 +------------------------------------------------------------------------------+
 |                                                                              |
 | Atomistix ToolKit 2017.0 [Build 15d690e]                                     |
 |                                                                              |
 +------------------------------------------------------------------------------+
 --------------------------------------------------------------------------------
  TransmissionSpectrum read from file :
  transmission_V-0.4.hdf5
  InelasticTransmissionSpectrum read from file :
  xloe_V-0.4_bulk.hdf5
 --------------------------------------------------------------------------------
  Current (w/o interactions)    = -3.71e-09 nA
  Current (w interactions)      = -6.42e-06 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 wide-band approximation. Phys. Rev. B, 89:081405, Feb 2014. doi:10.1103/PhysRevB.89.081405.