# Scanning Tunneling Microscopy {#sec:stm}

## Introduction

In this section, we present the steps to obtain the current flowing from
a sample to a tip using Bardeen's formalism. We refer the readers to
\[[GW06](#GW06){.citation}, [L14](#L14){.citation}\] for a basic review
of Bardeen's theory. The main assumptions underlying the theory:

1.  Tunneling is weak enough to employ first-order time-dependent
    perturbation theory;
2.  Tip and sample states are nearly orthogonal.;
3.  The electron-electron interaction can be ignored;
4.  The occupation of the tip and sample are independent of each other;
5.  The tip and sample are each in electrochemical equilibrium.

Using the above hypotheses, the following current equation is derived

[$$I = \frac{2\pi e}{\hbar}\sum\limits_{\mu\nu} \left[ f_\text{FD}(E^S_\mu-E^S_F) - f_\text{FD}(E^T_\mu-E^T_F)\right] |M_{\mu\nu}|^2 \delta(E^T_\mu - E^S_\nu - eV)$$]{label="eq:bardeenI"}

where $e$ is the charge of the electron, $\hbar$ is the reduced Planck
constant, $\mu$, $\nu$ are band indices, $f_\text{FD}$ is the
Fermi-Dirac function, $E^{S,T}_\mu$ is an eigenvalue of the sample or
tip, $E^{S,T}_F$ is Fermi level of the sample or tip and $V$ is the bias
applied between the sample and tip. Finally, $M_{\mu\nu}$ are the
tunneling matrix elements defined as

[$$M_{\mu\nu} = \frac{\hbar^2}{2m}\int\limits_{z=z_0} dx dy 
\left[ 
\psi^S_\mu \frac{\partial \psi^T_\nu}{\partial z} 
- 
\frac{\partial \psi^S_\mu}{\partial z} \psi^T_\nu 
\right]$$]{label="eq:bardeenM"}

In `eq:bardeenM`{.interpreted-text role="eq"}, the domain of the surface
integral is limited to a single unit cell.

RESCU evaluates the *continuum* version of
`eq:bardeenI`{.interpreted-text role="eq"},
`eq:bardeenJ`{.interpreted-text role="eq"}.

[$$I = \frac{2\pi e}{\hbar} \int d\epsilon 
\left[ 
f_\text{FD}(\epsilon + eV - E^S_F) 
- f_\text{FD}(\epsilon-E^T_F)
\right] 
n^S(\epsilon + eV) n^T(\epsilon) M(\epsilon+eV,\epsilon)$$]{label="eq:bardeenJ"}

This is done in an `stm-current` calculation, during which RESCU
computes `eq:bardeenJ`{.interpreted-text role="eq"} as follows:

1.  Solve the KS equation of the sample plus tip system (as a whole
    system). This provides the eigenvalues $E^{S,T}_\mu$, which are the
    same for sample and tip since the whole system is calculated, and
    eigenstates $\psi_\mu$.
2.  Compute the PDOS of the sample $n^S(\epsilon)$ and tip
    $n^T(\epsilon)$.
3.  Split the whole system wavefunctions into sample and tip
    wavefunctions $\psi = a^S \psi^S + a^T \psi^T$ by computing
    $\psi^{S,T} = \sum a^{S,T}_\eta \phi^{S,T}_\eta$ where
    $\phi^{S,T}_\eta$ are atomic orbital centred at atoms in the sample
    and tip respectively, and
    $a^{S,T}_\eta = \langle \phi_\eta | \psi \rangle$.
4.  Compute the tunneling matrix elements $M_{\mu\nu}$ using
    `eq:bardeenM`{.interpreted-text role="eq"}.
5.  Evaluate `eq:bardeenJ`{.interpreted-text role="eq"}, interpolating
    $M_{\mu\nu}$ to obtain $M(\epsilon+eV,\epsilon)$.

This new feature of RESCU has limited functionality, here are the main
constraints:

1.  It is implemented for atomic orbital calculations
    (`LCAO.status = 1`).
2.  The tip is above the sample in the +z direction.
3.  Their is a plane $z = z_0$ which separates the sample and tip.
4.  The k-sampling is limited to $\Gamma$ for now, because the feature
    was develop to study large realistic samples and tips.

## Example: graphene sample/tungsten tip

We now turn to a simple sample tip system to introduce the basic
keywords relevant to Bardeen current calculations. As customary, we
begin with the calculation of the ground state density. Let's create
input file `c2_lcao_scf.input` as follows

``` matlab
LCAO.status = 1
info.calculationType = 'self-consistent'
info.savepath = 'results/c2_lcao_scf'
info.outfile = 'fc2_lcao_scf.out'
la = 4.648725932
atom.element = [1 1 2]
atom.fracxyz = 'c2.xyz'
domain.latvec = [[1/2  sqrt(3)/2     0]*5
[1/2 -sqrt(3)/2     0]*5
[  0          0    27]/la]*la
domain.lowres = 0.3
element.species = 'C'
element.path = './C_LDA_TM.mat'
element(2).species = 'W'
element(2).path = './W_LDA_TM.mat'
domain.boundary = [1,1,2]
mixing.type = 'density'
```

The atomic coordinates can be generated using the following code

``` matlab
da = 3.0425/27;
xyz = [2,1,(0.5-da)*3;1,2,(0.5-da)*3]/3;
repxyz = [];
for ii = 0:4
for jj = 0:4
repxyz = [repxyz; bsxfun(@plus, xyz, [ii,jj,0])];
end
end
repxyz = repxyz/diag([5,5,1]);
fid = fopen('c2.xyz','w');
fprintf(fid,'%d\n\n',size(repxyz,1)+1);
fprintf(fid,'C %1.7e %1.7e %1.7e\n', repxyz'); 
fprintf(fid,'W %1.7e %1.7e %1.7e\n',[0,0,0.5+da]); 
fclose(fid);
```

You can visualize the system to make sure everything looks reasonable
and then execute RESCU typing

``` shell
rescu -s c2_lcao_scf.input
rescu -i c2_lcao_scf.input
```

Next, the input file of the `stm-current` calculation is defined as
`c2_lcao_stm.input` as follows

``` matlab
info.calculationType = 'stm-current'
info.savepath = 'results/c2_lcao_stm'
info.outfile = 'fc2_lcao_stm.out'
rho.in = 'results/c2_lcao_scf'
stm.z0 = 14 ! units of Bohr
stm.voltage = 1.2
units.energy = 'eV'
```

The important keywords of `c2_lcao_stm.input` are:

-   `info.calculationType` When set to `stm-current`, RESCU computes the
    current passing from the sample to the tip using
    `eq:bardeenJ`{.interpreted-text role="eq"}.
-   `rho.in` Provides the path the to files containing the ground state
    density.
-   `stm.z0` If `stm.z0` is set to $z_0$, then the tunneling matrix
    elements are computed on the plane $z = z_0$ using Eq.
    `eq:bardeenM`{.interpreted-text role="eq"}. Moreover, the system is
    decomposed into orbitals below (sample) and above (tip) the
    $z_0$-plane. Thence, the PDOS of the sample and tip ($n^S$ and
    $n^T$), and the sample and tip wavefunctions ($\psi^S$ and $\psi^T$)
    depend on the value of `stm.z0`.
-   `stm.voltage` Bias between the sample and tip.

Execute RESCU typing

``` shell
rescu -i c2_lcao_stm.input
```

Upon completion of the calculation, the results may be read from
`c2_lcao_stm.mat`. RESCU may also save temporary results to
`c2_lcao_stm_stm_tmp.mat`, or more generally to whereever the `savepath`
keyword points plus`stm_tmp` appended. In this way, computationally
expensive results can be reused to compute the current for different
bias voltages; note that you should begin by calculating the highest
bias voltage for that to work however. In particular, one will find the
structure `stm` which contains

-   `I` The current, in ampere (amp).
-   `M` The tunneling matrix elements.
-   `voltage` The bias between the sample and tip, converted to atomic
    units.
-   `z0` The separation surface $z$ value, converted to atomic units.

In the above example, you should get an answer around
$-1.46\times 10^{-5}$ amp.

::: {#citations}

[GW06]{#GW06 .citation-label}

:   Alex Gottlieb and Lisa Wesoloski. [Bardeen's tunnelling theory as
    applied to scanning tunnelling microscopy: A technical guide to the
    traditional interpretation.
    \<10.1088/0957-4484/17/8/R01\>]{.title-ref} Nanotechnology -
    NANOTECHNOL 17 (Apr. 2006).

[L14]{#L14 .citation-label}

:   Samir Lounis. [Theory of Scanning Tunneling Microscopy. 2014.
    \<10.1016/S0076-695X(08)60006-X\>]{.title-ref} Methods in
    Experimental Physics, Academic Press, Volume 27, 1993, Pages 1-29.
:::
