Photocurrent

class Photocurrent(configuration, energies=None, photon_energies=None, kpoints=None, photon_polarization=None, energy_resolution=None, calculate_all_transmissions=None)

Class for calculating the photocurrent and photon-mediated transmission in a device.

Parameters:
  • configuration (DeviceConfiguration) – The device configuration with an attached LCAOCalculator for which to calculate the photocurrent.

  • energies (list of PhysicalQuantity of type energy | Automatic) – The energies for which to calculate the photon-mediated transmission.
    Default: Automatic. Automatically setup energies based on the electrode Fermi levels and the photon energies.

  • photon_energies (list of PhysicalQuantity of type energy) – The photon energies to consider. Must be positive.
    Default: [1.0] * eV

  • kpoints (sequence (size 3) of int | MonkhorstPackGrid | KpointDensity | RegularKpointGrid) – The k-points for which to identify the band ladders. Note that the k-points must be in the xy-plane.
    Default: [0.0, 0.0, 0.0] (Gamma-point)

  • photon_polarization (list (size 3) of complex float) – The polarization of the photon illuminating the device configuration.
    Default: [0.0, 0.0, 1.0]

  • energy_resolution (PhysicalQuantity of type energy) – The distance between two energy points, when the energies are determined automatically. This parameter is only used, when energies=None or energies=Automatic
    Default: 0.05 * eV

  • calculate_all_transmissions (bool) – Boolean to specify if the photo-induced transmissions functions should be calculated at all possible energies. By default this is false, implying that the transmission functions only are calculated at the energies and photon energies, that would contribute to the total photo-current. If post-processing by shifting the Fermi level is wanted, the calculate_all_transmissions parameter can be set to True.
    Default: False

calculateAllTransmissions()
Returns:

Boolean controlling whether all transmissions are calculated.

Return type:

bool

electrodeFermiLevels()
Returns:

The electrode Fermi levels in absolute energy.

Return type:

list (size 2) of PhysicalQuantity of type energy

energies()
Returns:

The energies for which the photon-mediated transmission is calculated.

Return type:

list of PhysicalQuantity of type energy

energyResolution()
Returns:

The energy resolution used to setup the energies.

Return type:

PhysicalQuantity of type energy

energyZero()
Returns:

The energy zero used for the energy scale.

Return type:

PhysicalQuantity of type energy

evaluate(photon_modes=None, bias=None, temperature=None, fermi_shift=None, relative_permittivity=None, relative_permeability=None, photon_flux=None, spin=None, emission=None)

Evaluate the photocurrent.

Parameters:
  • photon_modes (All | int | list of int) – The indices of the photons for which to calculate the photocurrent. Must be non-negative.
    Default: All

  • bias (PhysicalQuantity of type voltage) – The voltage for which the current should be calculated.
    Default: The bias of the DeviceConfiguration.

  • temperature (PhysicalQuantity of type temperature) – The temperature fow which the current should be calculated.
    Default: 300*Kelvin

  • fermi_shift (PhysicalQuantity of type energy) – The shift of the Fermi level such that the current is evaluated at the energy E = fermi_energy + fermi_shift.
    Default: 0.0*eV

  • relative_permittivity (float) – The permittivity of the central region compared to vacuum. Must be positive.
    Default: 1.0

  • relative_permeability (float) – The permeability of the central region compared to vacuum. Must be positive.
    Default: 1.0

  • photon_flux (PhysicalQuantity of type 1/(area*time)) – The photon flux of the light incident on the central region.
    Default: 1.0*Second**-1*Angstrom**-2

  • spin (Spin.Up | Spin.Down | Spin.Sum) – The spin component for which to calculate the current.
    Default: Spin.Sum

  • emission (bool) – Boolean controlling whether photon emissions are included.
    Default: False

Returns:

The photocurrent as a list with the current for each of the photon energies.

Return type:

list of PhysicalQuantity of type current

kpoints()
Returns:

The fractional k-points used in the calculation.

Return type:

list of list (size 3) of float

maximumTemperature()
Returns:

The maximum trmperature used to setup the energies.

Return type:

PhysicalQuantity of type temperature

metatext()
Returns:

The metatext of the object or None if no metatext is present.

Return type:

str | None

nlprint(stream=None)

Print a string containing an ASCII table useful for plotting the AnalysisSpin object.

Parameters:

stream (python stream) – The stream the table should be written to.
Default: NLPrintLogger()

photonEnergies()
Returns:

The photon energies to consider.

Return type:

list of PhysicalQuantity of type energy

photonMediatedTransmission(photon_mode=None, carrier_type=None, initial_lead=None, final_lead=None, relative_permittivity=None, relative_permeability=None, photon_flux=None, spin=None)

Evaluate the photon-mediated transmission.

Parameters:
  • photon_mode (int) – The indices of the photons for which to calculate the photon-mediated transmission. Must be non-negative.
    Default: 0

  • carrier_type (Electron | Hole) – The carrier type for which to calculate the transmission.
    Default: Electron

  • initial_lead (Left | Right) – The lead that the unexcited carrier enters from.
    Default: Left

  • final_lead (Left | Right) – The lead that the excited carrier exits to.
    Default: Right

  • relative_permittivity (float) – The permitivity of the central region compared to vacuum. Must be positive.
    Default: 1.0

  • relative_permeability (float) – The permability of the central region compared to vacuum. Must be positive.
    Default: 1.0

  • photon_flux (PhysicalQuantity of type 1/(area*time)) – The photon flux of the light incident on the central region.
    Default: 1.0*Second**-1*Angstrom**-2

  • spin (Spin.Up | Spin.Down | Spin.Sum) – The spin component for which to calculate the transmission.
    Default: Spin.Sum

Returns:

The photon-mediated transmission.

Return type:

list (size len(kpoints), len(energies)) of float

photonPolarization()
Returns:

The photon polarization.

Return type:

list (size 3) of complex float

setMetatext(metatext)

Set a given metatext string on the object.

Parameters:

metatext (str | None) – The metatext string that should be set. A value of “None” can be given to remove the current metatext.

uniqueString()

Return a unique string representing the state of the object.

Usage Examples

The script below shows how to setup a Photocurrent calculation for a silicon p-n junction.

# -*- coding: utf-8 -*-
# -------------------------------------------------------------
# Two-probe Configuration
# -------------------------------------------------------------

# -------------------------------------------------------------
# Left Electrode
# -------------------------------------------------------------

# Set up lattice
vector_a = [3.84001408591, 0.0, 0.0]*Angstrom
vector_b = [0.0, 3.84001408591, 0.0]*Angstrom
vector_c = [0.0, 0.0, 10.8612]*Angstrom
left_electrode_lattice = UnitCell(vector_a, vector_b, vector_c)

# Define elements
left_electrode_elements = [Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon,
                           Silicon]

# Define coordinates
left_electrode_coordinates = [[  1.920007042956,  -0.            ,   0.678825      ],
                              [  1.920007042956,   1.920007042956,   2.036475      ],
                              [  0.            ,   1.920007042956,   3.394125      ],
                              [  0.            ,   0.            ,   4.751775      ],
                              [  1.920007042956,  -0.            ,   6.109425      ],
                              [  1.920007042956,   1.920007042956,   7.467075      ],
                              [  0.            ,   1.920007042956,   8.824725      ],
                              [  0.            ,   0.            ,  10.182375      ]]*Angstrom

# Set up configuration
left_electrode = BulkConfiguration(
    bravais_lattice=left_electrode_lattice,
    elements=left_electrode_elements,
    cartesian_coordinates=left_electrode_coordinates
    )

# Add external potential
external_potential = AtomicCompensationCharge([
    ('doping_0', 0.0200195107106),
    ('doping_1', -0.0200195107106)
    ])

left_electrode.setExternalPotential(external_potential)

# -------------------------------------------------------------
# Right Electrode
# -------------------------------------------------------------

# Set up lattice
vector_a = [3.84001408591, 0.0, 0.0]*Angstrom
vector_b = [0.0, 3.84001408591, 0.0]*Angstrom
vector_c = [0.0, 0.0, 10.8612]*Angstrom
right_electrode_lattice = UnitCell(vector_a, vector_b, vector_c)

# Define elements
right_electrode_elements = [Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon,
                            Silicon]

# Define coordinates
right_electrode_coordinates = [[  1.920007042956,  -0.            ,   0.678825      ],
                               [  1.920007042956,   1.920007042956,   2.036475      ],
                               [  0.            ,   1.920007042956,   3.394125      ],
                               [  0.            ,   0.            ,   4.751775      ],
                               [  1.920007042956,  -0.            ,   6.109425      ],
                               [  1.920007042956,   1.920007042956,   7.467075      ],
                               [  0.            ,   1.920007042956,   8.824725      ],
                               [  0.            ,   0.            ,  10.182375      ]]*Angstrom

# Set up configuration
right_electrode = BulkConfiguration(
    bravais_lattice=right_electrode_lattice,
    elements=right_electrode_elements,
    cartesian_coordinates=right_electrode_coordinates
    )

# Add external potential
external_potential = AtomicCompensationCharge([
    ('doping_0', 0.0200195107106),
    ('doping_1', -0.0200195107106)
    ])

right_electrode.setExternalPotential(external_potential)

# -------------------------------------------------------------
# Central Region
# -------------------------------------------------------------

# Set up lattice
vector_a = [3.84001408591, 0.0, 0.0]*Angstrom
vector_b = [0.0, 3.84001408591, 0.0]*Angstrom
vector_c = [0.0, 0.0, 76.0284]*Angstrom
central_region_lattice = UnitCell(vector_a, vector_b, vector_c)

# Define elements
central_region_elements = [Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon,
                           Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon,
                           Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon,
                           Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon,
                           Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon,
                           Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon,
                           Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon,
                           Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon]

# Define coordinates
central_region_coordinates = [[  1.920007042956,  -0.            ,   0.678825      ],
                              [  1.920007042956,   1.920007042956,   2.036475      ],
                              [  0.            ,   1.920007042956,   3.394125      ],
                              [  0.            ,   0.            ,   4.751775      ],
                              [  1.920007042956,  -0.            ,   6.109425      ],
                              [  1.920007042956,   1.920007042956,   7.467075      ],
                              [  0.            ,   1.920007042956,   8.824725      ],
                              [  0.            ,   0.            ,  10.182375      ],
                              [  1.920007042956,  -0.            ,  11.540025      ],
                              [  1.920007042956,   1.920007042956,  12.897675      ],
                              [  0.            ,   1.920007042956,  14.255325      ],
                              [  0.            ,   0.            ,  15.612975      ],
                              [  1.920007042956,  -0.            ,  16.970625      ],
                              [  1.920007042956,   1.920007042956,  18.328275      ],
                              [  0.            ,   1.920007042956,  19.685925      ],
                              [  0.            ,   0.            ,  21.043575      ],
                              [  1.920007042956,  -0.            ,  22.401225      ],
                              [  1.920007042956,   1.920007042956,  23.758875      ],
                              [  0.            ,   1.920007042956,  25.116525      ],
                              [  0.            ,   0.            ,  26.474175      ],
                              [  1.920007042956,  -0.            ,  27.831825      ],
                              [  1.920007042956,   1.920007042956,  29.189475      ],
                              [  0.            ,   1.920007042956,  30.547125      ],
                              [  0.            ,   0.            ,  31.904775      ],
                              [  1.920007042956,  -0.            ,  33.262425      ],
                              [  1.920007042956,   1.920007042956,  34.620075      ],
                              [  0.            ,   1.920007042956,  35.977725      ],
                              [  0.            ,   0.            ,  37.335375      ],
                              [  1.920007042956,  -0.            ,  38.693025      ],
                              [  1.920007042956,   1.920007042956,  40.050675      ],
                              [  0.            ,   1.920007042956,  41.408325      ],
                              [  0.            ,   0.            ,  42.765975      ],
                              [  1.920007042956,  -0.            ,  44.123625      ],
                              [  1.920007042956,   1.920007042956,  45.481275      ],
                              [  0.            ,   1.920007042956,  46.838925      ],
                              [  0.            ,   0.            ,  48.196575      ],
                              [  1.920007042956,  -0.            ,  49.554225      ],
                              [  1.920007042956,   1.920007042956,  50.911875      ],
                              [  0.            ,   1.920007042956,  52.269525      ],
                              [  0.            ,   0.            ,  53.627175      ],
                              [  1.920007042956,  -0.            ,  54.984825      ],
                              [  1.920007042956,   1.920007042956,  56.342475      ],
                              [  0.            ,   1.920007042956,  57.700125      ],
                              [  0.            ,   0.            ,  59.057775      ],
                              [  1.920007042956,  -0.            ,  60.415425      ],
                              [  1.920007042956,   1.920007042956,  61.773075      ],
                              [  0.            ,   1.920007042956,  63.130725      ],
                              [  0.            ,   0.            ,  64.488375      ],
                              [  1.920007042956,  -0.            ,  65.846025      ],
                              [  1.920007042956,   1.920007042956,  67.203675      ],
                              [  0.            ,   1.920007042956,  68.561325      ],
                              [  0.            ,   0.            ,  69.918975      ],
                              [  1.920007042956,  -0.            ,  71.276625      ],
                              [  1.920007042956,   1.920007042956,  72.634275      ],
                              [  0.            ,   1.920007042956,  73.991925      ],
                              [  0.            ,   0.            ,  75.349575      ]]*Angstrom

# Set up configuration
central_region = BulkConfiguration(
    bravais_lattice=central_region_lattice,
    elements=central_region_elements,
    cartesian_coordinates=central_region_coordinates
    )

# Add external potential
external_potential = AtomicCompensationCharge([
    ('doping_0', 0.0200195107106),
    ('doping_1', -0.0200195107106)
    ])

central_region.setExternalPotential(external_potential)

device_configuration = DeviceConfiguration(
    central_region,
    [left_electrode, right_electrode],
    equivalent_electrode_lengths=[10.8612, 10.861199999999997]*Angstrom,
    )

# Add tags
device_configuration.addTags('doping_0', [ 0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12,
                                          13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25,
                                          26, 27])
device_configuration.addTags('doping_1', [28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40,
                                          41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53,
                                          54, 55])

# -------------------------------------------------------------
# Calculator
# -------------------------------------------------------------
#----------------------------------------
# Basis Set
#----------------------------------------
basis_set = [
    BasisGGASG15.Silicon_Low,
    ]

#----------------------------------------
# Exchange-Correlation
#----------------------------------------
exchange_correlation = GGAHalf.PBE

#----------------------------------------
# Numerical Accuracy Settings
#----------------------------------------
left_electrode_k_point_sampling = MonkhorstPackGrid(
    na=4,
    nb=4,
    nc=87,
    )
left_electrode_numerical_accuracy_parameters = NumericalAccuracyParameters(
    k_point_sampling=left_electrode_k_point_sampling,
    density_mesh_cutoff=45.0*Hartree,
    )

right_electrode_k_point_sampling = MonkhorstPackGrid(
    na=4,
    nb=4,
    nc=87,
    )
right_electrode_numerical_accuracy_parameters = NumericalAccuracyParameters(
    k_point_sampling=right_electrode_k_point_sampling,
    density_mesh_cutoff=45.0*Hartree,
    )

device_k_point_sampling = MonkhorstPackGrid(
    na=4,
    nb=4,
    nc=87,
    )
device_numerical_accuracy_parameters = NumericalAccuracyParameters(
    k_point_sampling=device_k_point_sampling,
    density_mesh_cutoff=45.0*Hartree,
    )

#----------------------------------------
# Poisson Solver Settings
#----------------------------------------
left_electrode_poisson_solver = FastFourier2DSolver(
    boundary_conditions=[[PeriodicBoundaryCondition(),PeriodicBoundaryCondition()],
                         [PeriodicBoundaryCondition(),PeriodicBoundaryCondition()],
                         [PeriodicBoundaryCondition(),PeriodicBoundaryCondition()]]
    )

right_electrode_poisson_solver = FastFourier2DSolver(
    boundary_conditions=[[PeriodicBoundaryCondition(),PeriodicBoundaryCondition()],
                         [PeriodicBoundaryCondition(),PeriodicBoundaryCondition()],
                         [PeriodicBoundaryCondition(),PeriodicBoundaryCondition()]]
    )

#----------------------------------------
# Contour Integral Settings
#----------------------------------------
equilibrium_contour = SemiCircleContour(
    integral_lower_bound=1.43322362959*Hartree,
    )
contour_parameters = ContourParameters(
    equilibrium_contour=equilibrium_contour,
    )

#----------------------------------------
# Electrode Calculators
#----------------------------------------
left_electrode_calculator = LCAOCalculator(
    basis_set=basis_set,
    exchange_correlation=exchange_correlation,
    numerical_accuracy_parameters=left_electrode_numerical_accuracy_parameters,
    poisson_solver=left_electrode_poisson_solver,
    )

right_electrode_calculator = LCAOCalculator(
    basis_set=basis_set,
    exchange_correlation=exchange_correlation,
    numerical_accuracy_parameters=right_electrode_numerical_accuracy_parameters,
    poisson_solver=right_electrode_poisson_solver,
    )

#----------------------------------------
# Device Calculator
#----------------------------------------
calculator = DeviceLCAOCalculator(
    basis_set=basis_set,
    exchange_correlation=exchange_correlation,
    numerical_accuracy_parameters=device_numerical_accuracy_parameters,
    contour_parameters=contour_parameters,
    electrode_calculators=
        [left_electrode_calculator, right_electrode_calculator],
    )

device_configuration.setCalculator(calculator)
nlprint(device_configuration)
device_configuration.update()
nlsave('silicon_pn_photocurrent.hdf5', device_configuration)

# -------------------------------------------------------------
# Photocurrent
# -------------------------------------------------------------
kpoints = MonkhorstPackGrid(
    na=1,
    nb=1,
    )

photocurrent = Photocurrent(
    configuration=device_configuration,
    energies=numpy.linspace(-0.1, 4.1, 51)*eV,
    photon_energies=numpy.linspace(1, 4, 5)*eV,
    kpoints=kpoints,
    photon_polarization=[0+0j, 0+0j, 1+0j],
    )
nlsave('silicon_pn_photocurrent.hdf5', photocurrent)

Notes

Notice that in the script above the Photocurrent is only calculated for a single k-point and a few photon energies in order to have a simple example. For converged calculations one needs to increase the k-point sampling to something like:

kpoints = MonkhorstPackGrid(
  na=21,
  nb=21,
  )

It is important that the energies cover all possible transitions. Internally, the program will perform calculations at energy points \(E\) and \(E - \hbar\omega\), which means that the lower energy limit should be slightly below (due to thermal energy smearing) the lowest electrode chemical potential (in the example above this is 0.0*eV). The highest energy should be slightly larger than the highest photon energy, i.e., something like:

energies = numpy.linspace(- 0.1, max(photon_energies) + 0.1, number_of_points) * eV

The calculation time for a Photocurrent object depends on the number of k-points, energy points, and photon energy points. For a single photon energy, the calculation will take around two to three times the time for a normal TransmissionSpectrum calculation with the same energies and k-point sampling. However, for many photon energies the calculations are often rather time-consuming.

Theory

The photocurrent is calculated using first-order perturbation theory within the 1st Born approximation. The implementation follows the theory in Ref. [1] and the implementation described in Refs. [2] and [3].

In brief, the electron-light interaction is added to the Hamiltonian as:

\[\hat{H} = \hat{H}_0 + \frac{e}{m_0}\mathbf{A}\cdot\hat{\mathbf{p}},\]

where \(\hat{H}_0\) is the Hamiltonian without the electron-light interaction, \(e\) is the electron charge, \(m_0\) is the free electron mass, \(\hat{p}\) is the momentum operator, and \(\mathbf{A}\) is the electromagnetic vector potential. Assuming a single-mode monochromatic light source, we have:

\[\hat{A} = \left(\frac{\hbar\sqrt{\mu_r\epsilon_r}}{2N\omega\epsilon c} I_\omega \right)^{1/2}\left(\mathbf{a} \hat{b}e^{-i\omega t} + \mathbf{a}^*\hat{b}^\dagger e^{i\omega t} \right),\]

where \(\mathbf{a}\) is the (possibly complex) polarization vector, \(\mu_r\) is the relative magnetic susceptibility and \(\epsilon_r\) the relative dielectric constant, \(\epsilon\) is the dielectric constant, and \(\omega\) and \(c\) are the frequency and speed of light, respectively. \(\hat{b}^\dagger\) and \(\hat{b}\) are the photon creation and annihilation operators. \(I_\omega\) is the photon flux, defined as the number of photons (\(N\)) per unit time per unit area:

\[I_\omega = \frac{N c}{V \sqrt{\mu_u\epsilon_r}}.\]

The first-order coupling matrix now becomes:

\[M_{ln} = \frac{e}{m_0}\left(\frac{\hbar\sqrt{\mu_r\epsilon_r}}{2N\omega\epsilon c} I_\omega \right)^{1/2}\langle l | \hat{\mathbf{p}}\cdot \mathbf{a}| n \rangle ,\]

where \(|n\rangle\) is an LCAO basis function.

The 1st Born electron-photon self-energies are:

\[\mathbf{\Sigma}_{ph}^> = [N \mathbf{M}^\dagger \mathbf{G}_0^>(E^+)\mathbf{M} + (N + 1) \mathbf{M} \mathbf{G}_0^>(E^-)\mathbf{M}^\dagger]\]
\[\mathbf{\Sigma}_{ph}^< = [N \mathbf{M} \mathbf{G}_0^<(E^-)\mathbf{M}^\dagger + (N + 1) \mathbf{M}^\dagger \mathbf{G}_0^>(E^-)\mathbf{M}] ,\]

where \(E^{\pm} = E \pm \hbar\omega\). The Green’s function including electron-photon interactions to first order is then

\[\mathbf{G}^{>/<} = \mathbf{G}_0^r(\mathbf{\Sigma}^{>/<}_{L} + \mathbf{\Sigma}^{>/<}_{R} + \mathbf{\Sigma}_{ph}^{>/<})\mathbf{G}_0^a .\]

In the above two equations \(\mathbf{G}_0^{r,>,<}\) denote the non-interacting Green’s functions, and \(\mathbf{\Sigma}^{>/<}_{L,R}\) are the lesser and greater self-energies due to coupling to the electrodes. The current in electrode \(\alpha\) (left or right) with spin \(\sigma\) is calculated as:

\[I_{\alpha, \sigma} = \frac{e}{\hbar}\int \frac{dE}{2\pi}\sum_{k} T_\alpha(E, k, \sigma),\]

where the effective transmission coefficients are given by [3]:

\[T_\alpha(E, k, \sigma) = {\rm Tr}\left\{i\Gamma_\alpha(E, k) [1 - f_\alpha] G^< + f_\alpha G^> \right\}_{\sigma\sigma}.\]