# Photo-current tutorial

## Requirements

### Software components

- nanotools
- NanoDCAL+

### Pseudopotentials

- C_LDA_TM_DZP.mat

Copy the mat file from the pseudo database to the current working directory and `export NANODCALPLUS_PSEUDO=$PWD`.

## Theory

According to the first order perturbation theory, the charge current moving into a lead $L$ induced by optical excitation is calculated as

$$
J_L = \frac{ieN}{h} \int \textup{Tr}   \Gamma_L(E) (1-f_L(E)) G^r(E) M G^<(E-\hbar \omega) M^\dagger G^a(E) +  \textup{Tr}  \Gamma_L(E) f_L(E) G^r(E) M^\dagger G^>(E+\hbar \omega) M G^a(E)  dE
$$

where $N$ is number of photons, $\Gamma$ is broadening function of the lead, $f$ is Fermi-Dirac function,
$G$ is Green function, and $M$ is the matrix form of the polarization operator

$$
\hat{M} = \frac{e}{m_0} \sqrt{\frac{\hbar}{ 2 \omega \epsilon_r \epsilon_0 V}} \mathbf{\hat{e} \cdot \hat{p}}
$$

Here $\omega$ is the light frequency, $\epsilon_r\epsilon_0$ is permitivity of the material,
$V$ is volume of the optical cavity, $\hat{p}$ is quantum momentum operator, and $\mathbf{\hat{e}}$ is a complex polarization vector defined with an orthonormal pair of real vectors $\mathbf{(e_1,e_2)}$ and a pair of elliptical angles $(\theta,\phi)$:

$$
\mathbf{\hat{e}} = [\cos (\theta) \cos (\phi) - i \sin (\theta) \sin(\phi)]\mathbf{{e_1}}
+[\sin(\theta) \cos (\phi)+ i \cos(\theta) \sin(\phi)]\mathbf{{e_2}}
$$

To simplify the workflow of computing $J_L$, we define the dimensionless photo-transmission as follows

$$
\mathfrak{T}_L(E) =  \left ( \frac{\alpha c h a_0^2}{i \hbar \omega \epsilon_r V} \right )^{-1}  \textup{Tr}   [ \Gamma_L(E) (1-f_L(E)) G^r(E) M G^<(E-\hbar \omega) M^\dagger G^a(E) + \Gamma_L(E) f_L(E) G^r(E) M^\dagger G^>(E+\hbar \omega) M G^a(E) ]
$$

where $c$ is vacuum speed of light,
$a_0$ is Bohr radius,
$\alpha$ is fine-structure constant.
We further define the dimensionless photocurrent response coefficient

$$
\mathfrak{R}_L = \frac{\alpha}{\hbar \omega} \int \mathfrak{T}_L(E) dE
$$

Hence the photocurrent can be expressed as

$$ J*L = e a_0^2 \sqrt{\mu_r/\epsilon_r } I*\omega \mathfrak{R}\_L $$

where $I_\omega$ is the photon-flux in terms of particle density:

$$
I_\omega=\frac{N}{V}\frac{c}{\sqrt{\epsilon_r \mu_r}}
$$

For experimentalists, it may be more convenient to characterize photon-flux with the energy-flux parameter $P_\omega$,
i.e. optical power per area,
which has a simple relation with $I_\omega$: $P_\omega = \hbar \omega I_\omega$.

## Example

In this example, we calculate the photo-current through an electrically biased carbon chain.
Note that there should exist a charge current through the chain even without optical field.
Here the photo-current refers to the part of charge-current specifically induced by optical excitation.

python script

```python
from nanotools.photocurrent import PhotoCurrent
from nanotools import Atoms, Cell, System, TotalEnergy, TwoProbe
from nanotools.utils import ureg
import numpy as np

cell = Cell(avec=[12, 0.0, 0.0,
                  0.0, 12, 0.0,
                  0.0, 0.0, 16],
                  resolution = 0.12)
atoms= Atoms(positions="c8.xyz")
sys = System(cell=cell, atoms=atoms)
sys.kpoint.set_grid([1,1,40])
sys.pop.set_type("fd")
sys.pop.set_sigma(0.01)
center = TotalEnergy(sys)

# define the two-probe ground-state problem and solve
dev = TwoProbe(center=center, bias=0.5, transport_axis=2)
dev.set_mpi_command(mpi=f"mpiexec -n 4")
dev.left.solver.mpidist.kptprc = 4 # parallelize over k points
dev.left.solver.restart.densityPath = "left_out.h5"
dev.center.solver.restart.densityPath = "center_out.h5"
# dev.center.solver.mpidist.zptprc = 4 # parallelize over contour energy points
dev.solve()

#define the orthonormal basis vectors
base_vec = [  [-0.118915679855830, 0.316904693970085, -0.940973153721270],
              [-0.857210947821695, 0.445467497858298, 0.258356535211516] ]

#define the photo-current problem and solve
sct = PhotoCurrent.from_twoprobe(dev,
      polar_basis_vec = base_vec,
      theta = np.pi/2, phi= np.pi/4,
      elec_energy_resol = 0.01, #energy resolution for photo-transmission (the integrand)
      photon_energy = 0.7)
sct.solve()
sct.plot_photo_transmission(lead="left", save_to="phot_left.json")
sct.plot_photo_transmission(lead="right", save_to="phot_right.json")
r = sct.get_response()
c = sct.get_photo_current(
    flux=1.0 * ureg.watt / ureg.meter**2, # optical energy flux
    mu=1.0, epsilon=1.0) # relative material permitivity and permeability
print("photo-current response coefficient:")
print(r)
print("photo-current (density):")
print(c)
```

structure file c8.xyz, in Cartesian coordinates

```python
8
s x y z
C     6.00000000      6.00000000      1
C     6.00000000      6.00000000      3
C     6.00000000      6.00000000      5
C     6.00000000      6.00000000      7
C     6.00000000      6.00000000      9
C     6.00000000      6.00000000      11
C     6.00000000      6.00000000      13
C     6.00000000      6.00000000      15
```

On your screen, you should see the following output when the program is finished

```python
photo-current response coefficient:
[[[-14.87261766  14.87272017]]]
photo-current (density):
[[[-0.041317034921601804 0.04131731971157402]]] ampere / meter ** 2
```

Both response and current are 3d-arrays, with each dimension corresponding to photon-energy, spin, and lead.
In this example, the optical mode is assumed monocolor, and spin is degenerate.
Therefore, the output for photo-current here should be interpreted as
$J_L = -0.041317 A / m^2 \: J_R = 0.041317 A / m^2$.
The current (density) unit is $A / m^2$ because we apply periodic boundary condition across the transversal directions of the device structure.
The charge current through one single carbon chain should be calculated as $J_{L/R}$ multiplied by the cross section area of the simulation cell.

In addition, the following figure should appear in your work directory.
This figure plots $\mathfrak{T}_R$ as a function of electronic energies.

![Photo-transmission](../figs/photo_transmission.png)