5.1. Photocurrent in WSe2 monolayer

In this example, we will learn how to treat the phonon-assisted photon absorption [ZGC14]. We assume that the user is familiar with the structure builder from Device Studio and is able to set and perform the basic NEGF+DFT calculations. Otherwise, these instructions can be found in the previous tutorial WSe2 monolayer.

In this tutorial, we’ll demonstrate how to calculate the photocurrent through the WSe2 device in the armchair direction, as shown in Fig. 5.1.1. The WSe2 structure exhibit a honeycomb lattice and, similar to graphene, presents two well-separated local energy extrema (valleys) in the Brillouin zone, which are commonly labeled as K and K’ [MHS12]. Theoretical and experimental works have already demonstrated that a finite valley polarization could be created and detected by circularly polarized light (\(\sigma_{\pm}\)) over a WSe2 monolayer [CWH12]. Additional information could be found in the listed references of this tutorial.


Fig. 5.1.1 Schematic plot of the two-probe monolayer WSe2 along the armchair direction.

To compute the photocurrent calculation proceeds in three steps:

  1. Perform the transport calculation at zero-bias;

  2. Perform transport calculations under a finite source-drain bias voltage (\(V_{ds}\)) and/or gate voltage (\(V_{g}\));

  3. Consider a monochromatic light impinging on the device channel, where the electron-photon interaction is treated perturbatively.

5.1.1. The transport calculation at zero-bias.

In this step, we’ll perform the calculations for the electrodes and central region of the two-probe device in the armchair direction. Since the source and drain electrodes from the WSe2 junction are equal, we just have to calculate the electrode Hamiltonian for only one of them. The files are as follows:


In WSe2 Monolayer tutorial, the electronic transport was calculated along the zigzag direction (z-axis). Thus, the calculation along the armchair direction (x-axis) should be performed considering front and back semi-periodic electrodes relative to the WSe2 central scattering region. Please check the NanoDCAL input reference.

  • Front: [-1 0 0] means the direction of - x

  • Back: [ 1 0 0] means the direction of + x

Based on the knowledge and understanding that you already have gained from the previous examples you should perform the front electrode (source) and the two-probe transport calculation. The user should also be able to guarantee that the electronic properties from the electrode are well converged before the transport calculations, as we have done before.


After each calculation, the user has to go through the output file to ensure that everything is well behaved. Do not start any bias voltage (\(V_{ds}\)) calculation until you are perfectly sure that all results are well converged. Please check the set values at calculation.SCF.convergenceCriteria.

5.1.2. Finite-bias voltage calculations

Now we can explore the WSe2-armchair device under a finite bias voltage between the source and drain electrodes. For the purpose of this example, we’ll consider the following \(V_{ds}\) values of 0.1 to 0.5 V. The highest bias value has to be started using the closest bias calculations. This can be ensured by including additional keywords in the WSe2-CentralRegion.scf.input indicating the directory with previous voltage value, as showed below:

system.voltageOfLead1 = 0.0
system.voltageOfLead2 = 0.5
calculation.SCF.donatorObject = ../0.4Vb/NanodcalObject.mat

After running the bias voltage, the output file (NanodcalObject.mat) for each voltage value will be necessary for the next step.

5.1.3. Photocurrent calculation

Now, the calculated hamiltonian’s (\(H_{e}\)) will be used as the basis for the photocurrent calculation. The electron-photon interaction will be considered as a perturbation over the calculated \(H_{e}\):

\[H = H_{e}+\frac{e}{m}\mathbf{A} \cdot \mathbf{P},\]

where \(\mathbf{A}\) is the electromagnetic vector potential and \(\mathbf{P}\) the momentum of the electron [H02]. NanoDCAL calculates a normalized and dimensionless photocurrent:

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

where \(\alpha\) labels the source and drain contributions, \(\tau\) is the valley index defined in the reciprocal space, \(s\) is the spin index and \(T_\alpha (E,k,s)\) is the effective transmission coefficient.

The PhotoCurrent.input file can be created using Device Studio, following the steps below:

  1. Select the WSe2 device, go the the menu bar, click on SimulatorNanoDCALAnalysis to set physical quantity, as shown in Fig. 5.1.2:

  • In the Analysis window, on the top left choose the Device option;

  • Navigate to the Analysis panel (left side) ➟ select PhotoCurrent ➟ click on the Right arrow button to transfer PhotoCurrent calculator to the Calculation Selected panel.

  1. The PhotoCurrent options are shown in the right side of the Analysis window. Edit the parameters to fit the studied system:

  • Define the incident light propagation direction (-z direction), setting \(\mathbf{e_{1}}=[0,1,0]\) and \(\mathbf{e_{2}}=[1,0,0]\);

  • To model a circularly polarized light, navigate to the Analysis in Polarization vector ➟ select Oblique elliptical polarization ➟ set \(\theta = 90\) and \(\phi = 45\);

  • In Photon energy range set the values of \(E_{max} = E_{min} = 1.368\);

  • Go to the k-point sampling choose n1=50, n2=1, and n3=1;

  1. In the Analysis window, click over Edit irradiated atom button, select the atoms exclusively located at the central region (the atoms that comprises the electrode layers should not be illuminated). Check the selected atoms on Atoms indices ➟ click over select button ➟ ok and Generate files.

NanoDcal module simulator

Fig. 5.1.2 PhotoCurrent module simulator.

Now, we are ready to start the photocurrent calculation. Before doing that, we’ll briefly present the simulation parameters specified in the generated files.

PhotoCurrent input file

system.object = NanodcalObject.mat
calculation.name = photocurrent
%Parameters of the calculation
calculation.photocurrent.photonPolarizationBasis =[1 0 0; 0 1 0]
calculation.photocurrent.photonPolarizationAngleTheta = 90
calculation.photocurrent.photonPolarizationAnglePhi = 45
calculation.photocurrent.photonEnergyRange = [1.368821272021110 1.368821272021110]
calculation.photocurrent.photonEnergyInterval = 0.1
calculation.photocurrent.kSpaceGridNumber = [50 1 1]
calculation.photocurrent.irradiatedAtoms = [1:48]
%calculation.photocurrent.partitionScheme = {{}, {}}
calculation.photocurrent.plot = true

The keywords specify the following:

  • calculation.photocurrent.photonPolarizationBasis This input parameter specifies the two orthogonal real space direction vectors \(\mathbf{e_{1}}\) and \(\mathbf{e_{2}}\), which are used as basis to characterize the polarization of the photon. The lengths of \(\mathbf{e_{1}}\) and \(\mathbf{e_{2}}\) are meaningless,and they will be normalized by NanoDCAL.

  • calculation.photocurrent.photonPolarizationAngleTheta A angle used to characterize the photon polarization.

  • calculation.photocurrent.photonPolarizationAnglePhi A angle used to characterize the photon polarization. It is in the unit of degree.

  • calculation.photocurrent.photonEnergyRange Photocurrent is dependent on photon energy. This parameter defines a photon energy range [\(E_{1}\), \(E_{2}\)], and the photocurrent curve over the range will be calculated. Normally, \(E_{1}\) <\(E_{2}\), but if \(E_{1}= E_{2}\), the photocurrent is only calculated at a photon energy point of \(E_{2}\). The unit is \(eV\).

  • calculation.photocurrent.photonEnergyInterval Energy interval used to determine the photon energy points, and at each photon energy point the corresponding photocurrent will be calculated. The unit is \(eV\).

  • calculation.photocurrent.kSpaceGridNumber number of small \(k\)-space grids in each direction which, together with kSpaceGridShift, are used to produce the parameter kSpacePoints.

  • calculation.photocurrent.irradiatedAtoms Index of those atoms which are irradiated by the light.

  • calculation.photocurrent.partitionScheme Atoms in the system are divided into small groups so that atoms in each group interact only with atoms in the next group. This parameter gives the atom index of each of the atom groups.


  • The photon propagation direction is given by \(\mathbf{e_{1}} \times \mathbf{e_{2}}\), while the photon polarization is defined by,


where \(\theta\) and \(\phi\) are given as input parameters.

  • In this example we have calculated the photocurrent with a circularly polarized light, using the general elliptical polarization option. However, the user should note that different polarization vectors could be set for the incident light with NanoDCAL. Please, look at the Input Reference.

Run the simulation with WSe2-hpc.input file. After the calculation finish, the following output files are generated: CalculatedResults.mat, log.txt, Photocurrent.fig, and Photocurrent.mat.

The .mat results can be loaded using the MATLAB platform.

5.1.4. Photocurrent analysis

Note that the finite-bias calculations were restricted to \(V_{ds} < E_{gap}\) (\(E_{gap}= 1.37eV\)), which avoids the existence of current without the photon incidence. The valley current is obtained by summing the current contributions from the two spin index (\(s=\uparrow,\downarrow\))


The total current is obtained by a further summation over K and K’,


The user should note that the total current must be conserved (\(I_{S}=-I_{D}\)). However, a net valley-polarized current \(I_{S/D}^{\tau}\) is delivered to the source and drain. An additional key quantity could also be obtained from those current values, the valley polarization \(\eta\). At nonequilibrium \(V_{ds} \neq 0\), valley polarization becomes a function of the bias defined as

\[\eta(V_{ds})= \left| \frac{I_{\alpha,K} - I_{\alpha,K'}}{I_{\alpha,K} + I_{\alpha,K'}} \right|.\]

The script below calculates the \(I_{S/D,K}\), \(I_{S/D,K'}\), \(I_{S/D}^{\tau}\) and the valley polarization \(\eta\).

PhotoCurrent components script

% Script for obtain different component photocurrent
% NanoAcademic Technologies

clear all;clc

% Load photocurrent data
load ./Photocurrent.mat

% Get current information

photocurrent  = reshape(sum(data.kResolvedPhotocurrent, 4), 2, 50);

dk = data.coordinatesOfKPoints(1,2) - data.coordinatesOfKPoints(1,1);

photocurrentDKP = sum(photocurrent(2, 1:25), 2)*dk;
photocurrentDK = sum(photocurrent(2, 26:end), 2)*dk;

fprintf('I_{D, K} component is equal to %f \n', photocurrentDK)
fprintf(['I_{D, K''' '}component is equal to %f \n'], photocurrentDKP)

photocurrentSKP = sum(photocurrent(1, 1:25), 2)*dk;
photocurrentSK = sum(photocurrent(1, 26:end), 2)*dk;

fprintf('I_{S, K} component is equal to %f \n', photocurrentSK)
fprintf(['I_{S, K''' '}component is equal to %f \n'], photocurrentSKP)

photocurrentSTau = (photocurrentSK - photocurrentSKP)/2;

photocurrentDTau = (photocurrentDK - photocurrentDKP)/2;

fprintf('Valley current I_{S, tau} component is equal to %f \n', photocurrentSTau)
fprintf('Valley current I_{D, tau} component is equal to %f \n', photocurrentDTau)

Eta = abs((photocurrentDK - photocurrentDKP)/(photocurrentDK + photocurrentDKP));

fprintf(['Valley polarization is equal to %f \n'], Eta)

Fig. 5.1.3 Different valley components and total photocurrent in the source and drain.

From the Fig. 5.1.3, the net valley currents are nonzero and satisfy the continuity relation \(I_{S}^{\tau}=-I_{D}^{\tau}\). Since for \(\sigma_{-}\) light \(I_{S/D,K} \thickapprox 0\), the net valley current \(I_{S/D}^{\tau}\) delivered by the transistor is composed of electrons having quantum index K’. Similarly, the valley current generated by \(\sigma_{+}\) light is composed of electrons having index K. Hence, the transistor under a moderate bias (0.3 V) is efficient in operating as a valleytronic device [ZGC14].

  1. Zhang, K. Gong, J. Chen, L. Liu, Y. Zhu, D. Xiao, and H. Guo. Generation and transport of valley-polarized current in transition-metal dichalcogenides Phys. Rev. B 90 (2014), p. 195428.

    1. Henrickson. Nonequilibrium photocurrent modeling in resonant tunneling photodetectors J. Appl. Phys. 91 (2002), p. 6273.

    1. Mak, K. He, J. Shan, and T. F. Heinz. Control of valley polarization in monolayer MoS2 by optical helicity Nat. Nanotech. 7 (2012), p. 494.

  1. Cao, G. Wang, W. P. Han, H. Q. Ye, C. R. Zhu, J. R. Shi, Q. Niu, P. H. Tan, E. G. Wang, B. L. Liu, and J. Feng. Valley-selective circular dichroism of monolayer molybdenum disulphide Nat. Commun. 3 (2012), p. 887.