19. Scanning Tunneling Microscopy

19.1. 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, L14] 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

(19.1.1)\[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)\]

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

(19.1.2)\[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]\]

In (19.1.2), the domain of the surface integral is limited to a single unit cell.

RESCU evaluates the continuum version of (19.1.1), (19.1.3).

(19.1.3)\[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)\]

This is done in an stm-current calculation, during which RESCU computes (19.1.3) 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 (19.1.2).

  5. Evaluate (19.1.3), 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.

19.2. 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

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

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

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

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 (19.1.3).

  • 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. (19.1.2). 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

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 plusstm_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.

[GW06]

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> Nanotechnology - NANOTECHNOL 17 (Apr. 2006).

[L14]

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