# 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 (list of list (size 3) of float) – The k-points (in fractional coordinates) for which to identify the band ladders. 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. bool
electrodeFermiLevels()
Returns: The electrode Fermi levels in absolute energy. list (size 2) of PhysicalQuantity of type energy
energies()
Returns: The energies for which the photon-mediated transmission is calculated. list of PhysicalQuantity of type energy
energyResolution()
Returns: The energy resolution used to setup the energies. PhysicalQuantity of type energy
energyZero()
Returns: The energy zero used for the energy scale. 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 The photocurrent as a list with the current for each of the photon energies. list of PhysicalQuantity of type current
kpoints()
Returns: The fractional k-points used in the calculation. list of list (size 3) of float
maximumTemperature()
Returns: The maximum trmperature used to setup the energies. PhysicalQuantity of type temperature
metatext()
Returns: The metatext of the object or None if no metatext is present. str | unicode | 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. 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 The photon-mediated transmission. list (size len(kpoints), len(energies)) of float
photonPolarization()
Returns: The photon polarization. list (size 3) of complex float
setMetatext(metatext)

Set a given metatext string on the object.

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

## 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
)

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
)

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
)

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,
)

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. [Hen02] and the implementation described in Refs. [CHG12] and [ZGC+14].

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:

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

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

$\begin{split}\mathbf{G}^{>/<} = \mathbf{G}_0^r(\mathbf{\Sigma}^{>/<}_{L} + \mathbf{\Sigma}^{>/<}_{R} + \mathbf{\Sigma}_{ph}^{>/<})\mathbf{G}_0^a .\end{split}$

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 [ZGC+14]:

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

 [CHG12] Jingzhe Chen, Yibin Hu, and Hong Guo. First-principles analysis of photocurrent in graphene PN junctions. Phys. Rev. B, 85:155441, 2012. URL: https://link.aps.org/doi/10.1103/PhysRevB.85.155441, doi:10.1103/PhysRevB.85.155441.
 [Hen02] Lindor E. Henrickson. Nonequilibrium photocurrent modeling in resonant tunneling photodetectors. Journal of Applied Physics, 91(10):6273–6281, 2002. URL: http://aip.scitation.org/doi/abs/10.1063/1.1473677, doi:10.1063/1.1473677.
 [ZGC+14] (1, 2) Lei Zhang, Kui Gong, Jingzhe Chen, Lei Liu, Yu Zhu, Di Xiao, and Hong Guo. Generation and transport of valley-polarized current in transition-metal dichalcogenides. Phys. Rev. B, 90:195428, 2014.