Electron transport in ultra-thin films and ballistic electron emission microscopy

Electron transport at the nanoscopic scale is a field of high interest both at a fundamental and an applied level [1–5]. First, controlling electron motion at the nanometer lateral resolution is an important challenge. Second, spin-resolved electron currents are currently thoroughly investigated in epitaxial structures for their applications in giant magnetoresistance (GMR) devices. Significant advances in the study of magnetotransport phenomena at the nanometer scale through thin magnetic heterostructures have recently been obtained using ballistic electron emission microscopy (BEEM) [6, 7]. BEEM is an extension of scanning tunneling microscopy (STM) conceived to study metal-semiconductor interfaces [8, 9]. A remarkable demonstration of the STM potential to explore not only the surface region but interfaces buried beneath them too.

Initially, BEEM was devoted to measuring Schottky barriers at metal-semiconductor interfaces with nanometric lateral resolution. A basic theoretical description was presented by Kaiser and Bell based on a ballistic-electron model [8]. Ludeke and Prietsche subsequently improved the description of the metal-semiconductor interface by introducing the energy dependence of the metal-semiconductor transmission coefficient [10]. Both approaches take advantage of the free-electron energy dispersion for electrons propagating in the metal (${{\varepsilon}_{{\vec{k}}}}\propto {{k}^{2}}$). However, quickly became apparent that such a formalism could not explain the nanometric resolution found for some interfaces buried under thick layers. Such a key experimental feature was incorporated into the theory by imposing the propagation of carriers taking into account the periodic crystal lattice in the metallic layer and its corresponding electronic band structure [11].

In this formalism, the nanometric resolution is due to the existence of regions in the Fermi surface that force the carriers to propagate in focused beams. It is interesting to notice, however, that these focused beams coexist with others propagating accordingly to different patches of the same Fermi surface and resulting in extended s-like waves that spread over large regions. BEEM electrons are ultimately found in the detector in one state or the other depending on the experimental setup [12].

A different reason to introduce band structure effects arises from the conservation of the component of the wave vector parallel to the interface, ${{\overset{{}}{{k}}}_{//}}$, after transmission through an epitaxial interface with a good match between both materials. Pure ballistic electrons cannot pick up the high ${{\overset{{}}{{k}}}_{//}}$ value necessary to be injected through interfaces like, e.g. Au(1 1 1)/Si(1 1 1). These interfaces act as a filter for ${{\overset{{}}{{k}}}_{//}}$ values [6]. In the hypothesis of ${{\overset{{}}{{k}}}_{//}}$ conservation, band structure effects are the only possible mechanism to acquire the needed ${{\overset{{}}{{k}}}_{//}}$ for injection. Certainly, the ${{\overset{{}}{{k}}}_{//}}$ conservation depends on the quality of the interface; our aim is to analyze samples grown at the Institute of Physics of Rennes that have been assessed from an experimental point of view as of good enough quality to ensure a reasonable amount of ${{\overset{{}}{{k}}}_{//}}$ conservation [13].

Recent advances in the BEEM magnetotransport experiments call for a revision of the calculation scheme based on a decimation technique devised to target thick or semi-infinite slabs [14, 15]. Instead, a typical GMR device consists of intercalated magnetic and non-magnetic layers of few nanometers and is not well described by a semi-infinite slab [7]. To overcome this difficulty, we have developed a model for epitaxial heterostructures based on a layer-by-layer construction of the Hamiltonian, which deals with layered heterostructures of any composition and width. The present work is therefore devoted to the description of such a model adapted to ultra-thin slabs. As a future development, we shall extend this approach to take into account spin-polarization effects, which are necessary to describe magnetic materials.

The paper is organized as follows: in section 2 we introduce a layer-by-layer theoretical approach based on the Keldysh's non-equilibrium Green's function formalism. In particular, we write the current in terms of non-equilibrium Green functions, and we implement a layer-by-layer Green's functions approach based on Dyson's equation. Therefore, we can deal with metallic slabs of any thickness, even monolayers, as well as with slabs of alternating composition, e.g. Au/Fe/Au/Ag, etc. Currently, only elastic interactions have been taken into account—electron–phonon scattering and other losses mechanisms (impurities, etc) will be introduced in future work.

As an example of the method given in section 2, we present in section 3 numerical calculations for ultra-thin Au(1 1 1) layers. We highlight the new effects brought by the layer-by-layer formalism by comparing the BEEM current in ultra-thin Au(1 1 1) metal layers to propagation in thick Au(1 1 1) layers. We discuss the appearance of a BEEM current at the center of the reciprocal space ($\overline{\Gamma}$) along the Au(1 1 1) direction for a ten Au layers slab. This current is forbidden in thicker films due to the bulk band gap characteristic of coinage metals. Furthermore, we discuss the time-reversal symmetry of the Keldysh nonequilibrium approach and its breakdown as a function of the degree of decoherence. Finally, we analyze the analogies and differences of the full non-equilibrium Keldysh formalism with the semiclassical formalism.

Section 4 is devoted to a comparison with the semi-classical approach. Finally, we comment on the time-reversal symmetry of current distributions and its relation to quantum coherence.

A typical BEEM setup consists of a STM tip injecting carriers in a metal slab adsorbed on a semiconducting substrate that we assume thick enough to be considered semi-infinite, as depicted in figure 1.

Zoom In Zoom Out Reset image size

The Hamiltonian describing the tip (T), the sample (S, both the metal slab and the semiconducting substrate) and their interaction (I), considered at the level of electron hopping between the two, is given by:

$$\begin{eqnarray}\begin{array}{*{35}{l}} {} & H={{H}_{T}}+{{H}_{S}}+{{H}_{I}} \\ {} & {{H}_{T}}=\underset{\alpha}{\sum}{{\varepsilon}_{\alpha}}{{n}_{\alpha}} \\ {} & {{H}_{S}}=\underset{i}{\sum}{{\varepsilon}_{i}}{{n}_{i}}+\underset{i,j}{\sum}{{t}_{i,j}}c_{i}^{\dagger}{{c}_{j}}+\text{h}.\text{c}. \\ {} & {{H}_{I}}=\underset{\alpha,i}{\sum}{{t}_{\alpha,i}}c_{\alpha}^{\dagger}{{c}_{i}}+\text{h}.\text{c}. \end{array}\end{eqnarray} \tag{ 1 }$$

We use Greek letters for the tip and Latin letters for the sample (either the metal and/or the semiconductor). Creation, annihilation and number operators, ${{\hat{c}}^{\dagger}}$, $\hat{c}$ and $\hat{n}$, are defined in the usual way on a localized basis so the different Hamiltonians can be written in a straightforward way via a parametrized tight-binding formalism or via a more accurate ab initio density functional formalism using a suitable local combination of atomic orbitals as a basis set [16, 17]. Notice that the time dependence of all operators has been dropped in order to lighten the notation, and it is understood that each label i stands for the three cartesian coordinates necessary for a proper identification of the atomic site. Finally, we assume that the STM tip and the sample are separately at equilibrium and described by Hamiltonians H_T and H_S. We assume that the chemical potential of the tip, ${{\mu}_{T}}$, is different from the chemical potential of the sample, ${{\mu}_{S}}$, and we arbitrarily choose ${{\mu}_{T}}>{{\mu}_{S}}$, i.e. the tip is at a higher electric potential.

When the tip and the sample are coupled by the interaction, H_I, electrons flow from the tip to the sample and a BEEM current is established. The interaction term is treated as a perturbation through the hopping matrices ${{t}_{\alpha,i}}$ that link tunneling active atoms in the tip, α, with the corresponding ones in the sample, i. It should be noticed that, formally, this perturbative approach amounts to first establishing a DC bias (${{\mu}_{T}}>{{\mu}_{S}}$) and later the coupling between the tip and the sample, which is the correct time sequence for a BEEM experiment [18].

To orderly introduce the symbols we use in figure 1 we first re-express the formalism developed in [19] up to equation (8). Subsequently, we develop a layer-to-layer formalism in equations (9) to (15), which is the aim and the original contribution of this work.

The expression of the BEEM current at the plane n  +  1 in the heterostructure (see figure 1) can be written in terms of electrons and holes flowing in opposite directions as [18]:

$$\begin{eqnarray}&&\begin{array}{*{20}{l}}\langle {{\hat{J}}_{n+1}}\rangle =\frac{\text{i}e}{\hbar}\underset{m\leqslant n}{\sum}{{t}_{n+1,m}}\langle \hat{c}_{n+1}^{\dagger}\left(t+{{0}^{+}}\right){{\hat{c}}_{m}}(t)\rangle \\ \qquad \quad \ \;\,-\frac{\text{i}e}{\hbar}\underset{m\leqslant n}{\sum}{{t}_{m,n+1}}\langle \hat{c}_{m}^{\dagger}\left(t+{{0}^{+}}\right){{\hat{c}}_{n+1}}(t)\rangle\end{array} \end{eqnarray} \tag{ 2 }$$

Taking into account the definition of the non-equilibrium Green functions of the complete setup,

$$\begin{eqnarray}&&\mathcal{G}_{mn}^{+-}\left(t-{{t}^{\prime}}\right)=\text{i}\langle \hat{c}_{n}^{\dagger}\left({{t}^{\prime}}\right){{\hat{c}}_{m}}(t)\rangle \end{eqnarray}$$

we rewrite equation (2) as

$$\begin{eqnarray}\begin{array}{*{35}{l}} \langle {{{\hat{J}}}_{n+1}}\rangle & =\frac{e}{\hbar}\underset{m\leqslant n}{\sum}{{t}_{n+1,m}}\mathcal{G}_{m,n+1}^{+-}-{{t}_{m,n+1}}\mathcal{G}_{n+1,m}^{+-} \\ {} & =\frac{2e}{\hbar}\underset{m\leqslant n}{\sum}\Re \left({{t}_{n+1,m}}\mathcal{G}_{m,n+1}^{+-}\right) \end{array}\end{eqnarray} \tag{ 3 }$$

Green functions are evaluated in the limit ${{t}^{\prime}}\to {{t}^{+}}$ and in the last step we have used the property, valid at equal times and derived from the definition above (${{t}_{n+1,m}}\in R$)

$$\begin{eqnarray}&&\mathcal{G}_{mn}^{+-}=-{{\left(\mathcal{G}_{nm}^{+-}\right)}^{\ast}}\end{eqnarray}$$

In all these equations we have assumed a stationary condition. Therefore, the non-equilibrium Green functions depend only on the difference of the time variables. This hypothesis would not be true in general, under non-stationary and non-equilibrium conditions.

The Green function $\mathcal{G}_{m,n+1}^{+-}$ refers to the full Hamiltonian H and it can be evaluated applying perturbation theory with H_I as a perturbation by the usual Fourier transformed Dyson's expansion in the frequency domain (so that the time convolution becomes diagonal) to obtain:

$$\begin{eqnarray}&&\mathcal{G}_{m,n+1}^{+-}=G_{m,n+1}^{+-}+{{\left({{G}_{m,i}}{{t}_{i,\alpha}}{{\mathcal{G}}_{\alpha,n+1}}\right)}^{+-}}\end{eqnarray} \tag{ 4 }$$

Sum over repeated indices of consecutive matrices has been assumed, and for the meaning of the different symbols used in the Green functions we refer to the caption of figure 1.

Multiple products of non-equilibrium Green functions, like the term ${{\left({{G}_{m,i}}{{t}_{i,\alpha}}{{\mathcal{G}}_{\alpha,n+1}}\right)}^{+-}}$, represent complex-time contour integrations and can be re-expressed in terms of equilibrium and non-equilibrium Green functions by means of Langreth's theorem [20]. Here, the perturbation is instantaneous, implying that the two branches of the complex-time contour are not linked, i.e. the hopping t⁺⁻ is zero. Therefore:

$$\begin{eqnarray}&&{{\left({{G}_{m,i}}{{t}_{i,\alpha}}{{\mathcal{G}}_{\alpha,n+1}}\right)}^{+-}}=G_{m,i}^{R}{{t}_{i,\alpha}}\mathcal{G}_{\alpha,n+1}^{+-}+G_{m,i}^{+-}{{t}_{i,\alpha}}\mathcal{G}_{\alpha,n+1}^{A}\end{eqnarray} \tag{ 5 }$$

To close Dyson's equation (4) we need one more iteration to express the Green's function ${{\mathcal{G}}_{\alpha,n+1}}$, which relates the tip and the sample, in terms of the Green's function ${{\mathcal{G}}_{m,n+1}}$, which is fully internal to the sample (Latin indices). We have:

$$\begin{eqnarray}&&\mathcal{G}_{\alpha,n+1}^{+-}={{\left({{G}_{\alpha,\beta}}{{t}_{\beta,i}}{{\mathcal{G}}_{i,n+1}}\right)}^{+-}}=G_{\alpha,\beta}^{R}{{t}_{\beta,i}}\mathcal{G}_{i,n+1}^{+-}+G_{\alpha,\beta}^{+-}{{t}_{\beta,i}}G_{i,n+1}^{A}\end{eqnarray} \tag{ 6 }$$

We notice that the zeroth-order term (${{G}_{\alpha,n+1}}$) is missing because the tip and the sample are decoupled in the absence of H_I; equation (6) is obtained from Langreth's theorem, equation (5). Finally, the current at layer n  +  1 can be written as an integral over a window of energies from Schottky's barrier ${{\phi}_{\text{SB}}}$ to the applied voltage, $\text{eV}$ [21]:

$$\begin{eqnarray}&&\begin{array}{*{20}{l}}{{J}_{n+1}}(\text{eV},{{\vec{k}}_{//}})=\frac{4e\pi}{\hbar}(\,{{f}_{T}}-{{f}_{S}})\underset{m\ne n+1}{\sum}\Im \\ \quad \times {\int}_{{{\phi}_{\text{SB}}}}^{\text{eV}}\frac{\text{d}E}{2\pi}\text{Tr}\left[{{t}_{n+1,m}}\mathcal{G}_{m,1}^{R}{{t}_{1,0}}{{\rho}_{0,0}}{{t}_{0,1}}\mathcal{G}_{1,n+1}^{A}\right]\end{array}\end{eqnarray} \tag{ 7 }$$

Here f_T and f_S are the Fermi factors for the tip and the sample, respectively, and ${{\rho}_{0,0}}$ is the density of states of the uncoupled tip. This equation re-expresses the one given in [21] in terms of our notation of figure 1.

We can further simplify equation (7) with the assumption that the coupling matrices t_0,1 are much smaller than hopping matrices in the metal. This usually applies to the tunneling conditions of BEEM, where a typical distance between the sample and the tip is  ∼5 Å. In this case, we can work at the lowest-order perturbation theory and use for the sample the un-renormalized Green functions g instead of G (i.e. we neglect the influence of the tip on the sample). Finally, assuming zero temperature:

$$\begin{eqnarray}&&\begin{array}{*{20}{l}}{{J}_{n+1}}(\text{eV},{{\vec{k}}_{//}})=\frac{4e\pi}{\hbar}\underset{m\ne n+1,m<n+1}{\sum}\\ \qquad \qquad \quad \times \Im {\int}_{{{\phi}_{\text{SB}}}}^{\text{eV}}\frac{\text{d}E}{2\pi}\text{Tr}\left[{{t}_{n+1,m}}G_{m,1}^{R}{{t}_{1,0}}{{\rho}_{0,0}}{{t}_{0,1}}G_{1,n+1}^{A}\right]\end{array}\end{eqnarray} \tag{ 8 }$$

The non-equilibrium current has been finally expressed in terms of equilibrium Green's functions of the uncoupled systems, quantities which can be calculated with the usual equilibrium formalism. In order to apply this formalism to ultra-thin slabs, we have chosen to evaluate equilibrium Green's functions $G_{m,1}^{R}$ and $G_{1,n+1}^{A}$ of equation (8) via a layer-by-layer procedure, therefore avoiding the decimation technique used previously [14].

To this purpose, Green's functions are calculated iteratively by coupling together the $(n+1)\text{th}$-layer of unperturbed hamiltonian H_n+1 with the first n layers, described by the hamiltonian H_n:

$$\begin{eqnarray}&&{{H}_{0}}=\overset{{{H}_{n}}}{{\overbrace{{\underset{{{{\vec{k}}}_{//}}_{{}}}{\sum}\underset{i=1}{\overset{n}{\sum}}\varepsilon _{{{{\vec{k}}}_{//}}}^{(i)}\hat{c}_{i{{{\vec{k}}}_{//}}}^{\dagger}{{{\hat{c}}}_{i{{{\vec{k}}}_{//}}}}+\underset{{{{\vec{k}}}_{//}}}{{{\sum}^{}}}\underset{i\ne j}{\overset{n}{\sum}}{{t}_{ij{{{\vec{k}}}_{//}}}}\hat{c}_{i{{{\vec{k}}}_{//}}}^{\dagger}{{{\hat{c}}}_{j{{{\vec{k}}}_{//}}}}}}}}\quad +\overset{{{H}_{n+1}}}{{\overbrace{{\underset{{{{\vec{k}}}_{//}}}{{{\sum}^{}}}\varepsilon \hat{c}_{n+1{{{\vec{k}}}_{//}}}^{\dagger}{{{\hat{c}}}_{n+1{{{\vec{k}}}_{//}}}}}}}}\end{eqnarray} \tag{ 9 }$$

Here n can be as small as 1 for a bi-layer. Their interaction up to next-nearest-layer hopping, H_I, is described by:

$$\begin{eqnarray}&&{{H}_{I}}=\underset{H_{I}^{\left(n,n+1\right)}}{{\underbrace{{\underset{{{{\vec{k}}}_{//}}}{{{\sum}^{}}}{{t}_{n,n+1,{{{\vec{k}}}_{//}}}}\hat{c}_{n,{{{\vec{k}}}_{//}}}^{\dagger}{{{\hat{c}}}_{n+1,{{{\vec{k}}}_{//}}}}}}}}+\underset{H_{I}^{\left(n-1,n+1\right)}}{{\underbrace{{\underset{{{{\vec{k}}}_{//}}}{{{\sum}^{}}}{{t}_{n-1,n+1,{{{\vec{k}}}_{//}}}}\hat{c}_{n-1,{{{\vec{k}}}_{//}}}^{\dagger}{{{\hat{c}}}_{n+1,{{{\vec{k}}}_{//}}}}}}}}\end{eqnarray} \tag{ 10 }$$

A Fourier transform in ${{\vec{k}}_{//}}$ has been performed in the above equations to take into account the two-dimensional intra-layer periodicity. We notice that the calculation of Green's function of order n  +  1 proceeds iteratively, starting from the Green's function of two uncoupled layers, coupling them through H_I and using the result for the two coupled layers as an input for the calculation of n  +  1  =  3. And so on. The complete procedure, implemented in a FORTRAN code (BEEM v3.0), shall be published elsewhere.

In this scheme, we want to evaluate the current that flows from the last layer of a slab (n, the injector) to the first layer of another slab (n  +  1, the collector). To this aim we write the advanced and retarded Green's functions appearing in equation (8), $G_{1,n+1}^{A}$ and $G_{m,1}^{R}$, in terms of the separate Green's functions for the injector, g⁽ⁱ⁾, and for the collector, g^(c). The corresponding Dyson's equations can be written as follows:

$$\begin{eqnarray}\left\{\begin{array}{r l} G_{1,n+1}^{A} & =g_{1,n}^{A,(i)}\,{{t}_{n,n+1}}\,G_{n+1,n+1}^{A} \\ G_{n+1,n+1}^{A} & =g_{n+1,n+1}^{A,(c)}+g_{n+1,n+1}^{A,(c)}\,{{t}_{n+1,n}}\,G_{n,n+1}^{A} \\ G_{n,n+1}^{A} & =g_{n,n}^{A,(i)}\,{{t}_{n,n+1}}\,G_{n+1,n+1}^{A} \end{array}\right. \end{eqnarray} \tag{ 11 }$$

and

$$\begin{eqnarray}\left\{\begin{array}{r l} G_{n,1}^{R} & =g_{n,1}^{R,(i)}+g_{n,n}^{R,(i)}\,{{t}_{n,n+1}}\,G_{n+1,1}^{R} \\ G_{n+1,1}^{R} & =g_{n+1,n+1}^{R,(c)}\,{{t}_{n+1,n}}\,G_{n,1}^{R} \end{array}\right. \end{eqnarray} \tag{ 12 }$$

Solving the above equation system, equations (11) and (12), we obtain (I is the identity matrix):

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} G_{1,n+1}^{A}=g_{1,n}^{A,(i)}\,{{t}_{n,n+1}}\,\times \overset{D_{n+1,n+1}^{A}}{{\overbrace{{{{\left[I-g_{n+1,n+1}^{A,(c)}\,{{t}_{n+1,n}}\,g_{n,n}^{A,(i)}\,{{t}_{n,n+1}}\right]}^{-1}}}}}}g_{n+1,n+1}^{A,(c)} \end{array}\end{eqnarray} \tag{ 13 }$$

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} G_{n,1}^{R}=\underset{D_{n,n}^{R}}{{\underbrace{{{{\left[I-g_{n,n}^{R,(i)}\,{{t}_{n,n+1}}\,g_{n+1,n+1}^{R,(c)}\,{{t}_{n+1,n}}\right]}^{-1}}}}}}g_{n,1}^{R,(i)} \end{array}\end{eqnarray} \tag{ 14 }$$

To lighten the notation, in equations (11)–(14) we have considered only nearest-neighbor hopping. Generalization of these expressions to longer range hopping, however, is straightforward and for completeness it has been included in the following equations. We can now substitute these propagators in equation (8), to obtain the final formula for the current:

$$\begin{eqnarray}&&\begin{array}{*{20}{l}}{{J}_{n+1}}(\text{eV},{{\vec{k}}_{//}})\\ \quad \,=\frac{4e\pi}{\hbar}\underset{(l,m)<n+1}{\sum}\Im {\int}_{{{\phi}_{\text{SB}}}}^{\text{eV}}\frac{\text{d}E}{2\pi}\text{Tr}[{{t}_{n+1,l}}\,D_{l,l}^{R}g_{l,1}^{R,(i)}\,{{t}_{1,0}}{{\rho}_{00}}{{t}_{0,1}}\\ \quad \,g_{1,m}^{A,(i)}{{t}_{m,n+1}}D_{n+1,n+1}^{A}g_{n+1,n+1}^{A,(c)}]\end{array}\end{eqnarray} \tag{ 15 }$$

This equation can be written in a more symmetric way by making the following transformations: (1) ${{\Gamma}_{1,1}}={{t}_{1,0}}{{\rho}_{00}}{{t}_{0,1}}$. (2) $\Im [\text{i}\pi g_{n+1,n+1}^{A(c)}]={{\rho}_{n+1,n+1}}$. (3) Expanding the denominator in equation (13) and rearranging terms it can be shown that ${{t}_{m,n+1}}D_{n+1,n+1}^{A}=D_{m,m}^{A}{{t}_{m,n+1}}$. (4) ${{\Gamma}_{m,l}}={{t}_{m,n+1}}{{\rho}_{n+1,n+1}}{{t}_{n+1,l}}$. Then,

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} {{J}_{n+1}}(\text{eV},{{{\vec{k}}}_{//}})=\frac{4e\pi}{\hbar}\underset{(l,m)<n+1}{\sum}{\int}_{{{\phi}_{\text{SB}}}}^{\text{eV}}\frac{\text{d}E}{2\pi}\text{Tr} \\ \qquad \qquad \,\,\,\,\times \,\left[D_{l,l}^{R}g_{l,1}^{R(i)}\,{{\Gamma}_{1,1}}\,g_{1,m}^{A(i)}D_{m,m}^{A}{{t}_{m,n+1}}{{\rho}_{n+1,n+1}}{{t}_{n+1,l}}\right] \\ \,\,\,\qquad \quad \ =\frac{4e}{\hbar}\underset{(l,m)<n+1}{\sum}{\int}_{{{\phi}_{\text{SB}}}}^{\text{eV}}\frac{\text{d}E}{2\pi}\text{Tr}\left[D_{l,l}^{R}g_{l,1}^{R(i)}\,{{\Gamma}_{1,1}}\,g_{1,m}^{A(i)}D_{m,m}^{A}\,{{\Gamma}_{m,l}}\right] \end{array}\end{eqnarray} \tag{ 16 }$$

This equation is similar to Fisher–Lee's equation if we interpret the $\Gamma$ as injection rate matrices. In fact, Fisher–Lee's equation is simply recovered by substituting the multiple-scattering renormalized propagators by their zeroth order approximation, i.e. by replacing D by the identity matrix, which is valid if the perturbation between the collector and the injector is weak [22, 23]

$$\begin{eqnarray}&&{{J}_{n+1}}(\text{eV},{{\vec{k}}_{//}})=\frac{4e}{\hbar}{\int}_{{{\phi}_{\text{SB}}}}^{\text{eV}}\frac{\text{d}E}{2\pi}\text{Tr}\left[g_{n,1}^{R(i)}\,{{\Gamma}_{1,{{1}^{\prime}}}}\,g_{{{1}^{\prime}},{{n}^{\prime}}}^{A(i)}\,{{\Gamma}_{{{n}^{\prime}},n}}\right]\end{eqnarray} \tag{ 17 }$$

where we have used primed indices to remind the reader about the contraction of indices in consecutive matrices.

In this expression, if we assume to have one orbital per layer only, since ${{\Gamma}_{1,1}}$ is a term that can be taken as a constant, it is clear that the current arriving at the collector (layer n  +  1) is basically governed by the probability density derived from the propagator in the injector $|{{g}_{n,1}}{{|}^{2}}$. Finally, by taking

$$\begin{eqnarray}&&{{T}_{n+1,0}}=2\pi \rho _{n+1,n+1}^{\frac{1}{2}}{{t}_{n+1,n}}g_{n,1}^{R}{{t}_{1,0}}\rho _{0,0}^{\frac{1}{2}}\end{eqnarray}$$

it can be seen that our expression reduces to Landauer's for the conductance [24]:

$$\begin{eqnarray}&&\frac{\partial J\left(\text{eV}\right)}{\partial \left(\text{eV}\right)}=\frac{2{{e}^{2}}}{h}~\text{Tr}\left[T{{T}^{\dagger}}\right]\left(\text{eV}\right)\end{eqnarray} \tag{ 18 }$$

Our next section is devoted to testing and illustrating this formalism. A previous version of the BEEM code has been based on a transfer-matrix formalism and the decimation technique; it propagated carriers inside the metal base using Green's functions computed for a semi-infinite slab [25]. On the other hand, the current version of the BEEM code (v3.0) can be used to describe propagation in ultra-thin layers. To illustrate the differences between both approaches we discuss three results: (1) the appearance of BEEM current around ${{\overset{{}}{{k}}}_{//}}=0$ in the (1 1 1) direction, which is forbidden for thick layers due to the bulk band gap characteristic of coinage metals. (2) The effect of the loss of time-reversal symmetry as the number of planes are increased. (3) Finally, we compare and test the simple semi-classical approach with the full quantum-mechanical we have developed here.

The formalism we are presenting here improves upon the former in two important points: (i) it describes thin films with finite widths, and (ii) it describes the propagation through the interface between the metal and the semiconductor taking into account the atomic structure of the interface. The latter is important for comparing with experiments but introduces the obvious complication of describing the metal-semiconductor interface at an atomic level. As an example, the Au/Si interface between two (1 1 1) surfaces, known to be of fair quality, implies maching a $3\times 3$ (Si) with a $4\times 4$ (Au) 2D unit cell, which makes a tedious calculation only necessary if a direct comparison with experiments is desired. Our goal merely is to present and test this formalism. Therefore we make approximations to simplify the description of the interface, while we keep focusing on the former one. To work with a $1\times 1$ 2D unit cell, we describe the semiconductor by its density of states ρ, and we take ${{\Gamma}_{{{n}^{\prime}},n}}$ as ${{T}^{2}}\rho {{\delta}_{{{n}^{\prime}},n}}$, T being an effective hopping between the injector and the collector. This approximation is reminiscent of a similar one that has been tested before in the context of the STM with good results [26]. Incorporating this simplification, the BEEM current—equation (17)– is proportional to $g_{n,1}^{R(i)}{{\Gamma}_{1,{{1}^{\prime}}}}g_{{{1}^{\prime}},n}^{A(i)}$; this indicates that the propagator in the metal is the relevant quantity to describe the transport.

As commented above, this equation is valid when the interaction between the metal and the semiconductor is weak, and the multiple scattering in the interface can be neglected. Our results in figures 2 and 3 are related to such scenario. The following level of sophistication can be achieved by reference to figure 4, where the projected density of states of the semiconductor has been plotted. The current distribution on the last plane of the metal can then be propagated to the semiconductor by taking into account the conduction bands of the semiconductor and an appropriate energy and ${{\vec{k}}_{//}}$ dependent transmission coefficient [19]. A word of caution is in order here. Figure 4 shows that to work with the Au(1 1 1) 1x1 periodicity there is a price to pay. The mismatch between Au and Si projected lattices forces us to consider patches of density in the semiconductor inside the 1st Brillouin zone of Au, ${{\rho}_{n+1,n+1}}({{\vec{k}}_{//}}+\vec{G})$, where $\vec{G}$ is a vector of the projected 2D reciprocal lattice of Si. These regions are integrated with a transmission coefficient that carries the ${{\vec{k}}_{//}}$ dependence computed for Au, which is a valid approximation only as long as not a too high intensity is injected this way. In fact, this is a practical point that does not concern us for the cases we are presenting since we are not yet attempting to include the semiconductor in an experimentally useful way.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Naturally, the way to improve this description is to move to equation (16) to account for multiple-scattering effects via the denominators in equations (13) and (14). In a real system, the effect of these factors is not only to renormalize the propagators, but also to modify the ${{\vec{k}}_{//}}$ condition at the interface by mixing the $\vec{G}$ vectors of the different periodicities at both sides of the interface.

In all the examples given in this paper, Hamiltonians in equation (1) have been obtained from a suitable tight-binding parametrization for Au [27], an approximation accurate enough for our purposes.

The study of the Au(1 1 1)/Si(1 1 1) and Au(1 1 1)/Si(0 0 1) interfaces was instrumental to demonstrating the critical role of the electronic band structure in BEEM propagation, as opposed to the initial use of a free-electron dispersion [11, 21]. In this section, we revisit the Au(1 1 1)/Si heterostructures to describe the evolution of the BEEM current for ultra-thin slabs as a function of the number of layers by using the new layer-by-layer approach developed in section 2.

Figure 2 shows the evolution of the current profiles along the $\bar{M}\bar{\Gamma}\bar{M}$ direction of the Au(1 1 1) 2D Brillouin zone when the number of layers is progressively doubled at each step from 6 to 768 layers. As the number of layers is increased, more poles in the Green's functions contribute to the current and the profile is further smoothed. Notice that whereas the calculation of the surface DOS converged after 20 layers only (not shown here), the convergence of the current is only reached between 700 and 800 layers. We remind that we consider the convergence as reached when our layer-by-layer construction reproduces the results of the semi-infinite decimation technique. The reason for the different convergence properties can be understood as follows: the surface DOS is extracted from the surface Green function $\hat{G}_{1,1}^{(n)}$ of a n-layer slab that after n  =  20 layers is effectively uncoupled from the other side of the slab. On the contrary, the calculation of the current is based on Green's functions of the type $\hat{G}_{1,n}^{(n)}$ because the current propagates throughout the material (from layer 1 to layer n) and therefore is more sensitive to deep regions in the slab.

Interestingly enough, in figure 2 it is observed that a peak appears at ${{\vec{k}}_{//}}=0$ for twelve layers, but it is not present in any of the other cases. To clarify the conditions that origin this peak we have performed layer-by-layer calculations from seven to eighteen layers, see figure 3. In this figure, we use a smaller broadening than in figure 2, $\eta =0.005$ eV, to reduce the smoothing of the peaks which obscures their identification. The ${{\vec{k}}_{//}}=0$-peak under study appears for ten layers, reaches its maximum for eleven layers and disappears after twelve layers under these conditions. This result is significant since it can alter the filtering effect of the metal base. If the available density of states in the semi-conductor is located mainly around ${{\vec{k}}_{//}}=0$ then a slab with the right thickness would permit the transport of a BEEM current to the Si substrate. We associate this current to a quantum finite-size effect; therefore, it must be sensitive not only to the actual size of the slab but also to the precise boundary conditions at the interface, in particular, the matching with the semiconductor surface.

This effect can be better understood by comparing the BEEM current for Au(1 1 1)/Si(1 1 1) and Au(1 1 1)/Si(0 0 1) for a single energy (+1 eV above the Fermi energy) for two thicknesses of the metal base: 10 and 800 layers of Au, respectively. This concept can be illustrated in the simplest approximation; we follow the procedure in [21], and we compare these currents with the available density of states in the projected 2D Brillouin zones for Si(0 0 1) and Si(1 1 1) respectively. The result is shown in figure 4. The 800 layer case is very similar to what had been previously obtained through the decimation technique and BEEM currents for both orientations of silicon are qualitatively similar to each other as they overlap high current areas in an analogous way. However, for 10 layers the result is different (figures 4(c) and (d)). The current distribution for these thin layers is quasi-discrete, and there is practically no overlap with the available density of states, except near ${{\vec{k}}_{//}}=0$ for the (0 0 1) orientation of silicon. In other words, for ultra-thin layers electrons would only be injected in through the epitaxial Si(1 0 0) interface, not through the Si(1 1 1) one. This is an interesting example showing how thin films may behave very differently from thick films regarding ${{\vec{k}}_{//}}$-filtering properties. A word of caution is in order here: as we have commented above the critical thickness for which this effect might be observed may depend on the detailed matching of electronic states at the interface, and on the hypothesis that inelastic losses can be neglected. Such a condition that can be favored by low temperatures (to avoid scattering with phonons), high material purity (to avoid scattering with defects), and energies close to the Schottky barrier (where the main contribution to the BEEM current is the elastic component of the electron beam).

4.1. Semi-classical

As mentioned in the introduction, the BEEM current was initially modelled by Kaiser and Bell as a ballistic process based on a quadratic free-electron dispersion law for the carriers [8]. This picture was refined by the work of Garcia-Vidal it et al that argued that it is not possible to neglect the periodic potential of the crystal, responsible for the main contribution to the elastic scattering of BEEM electrons [11]. In that paper, the Keldysh formalism has been used to derive an expression for the current (similar to equation 8), but Green's functions have been computed in a semi-classical approximation assuming the specular reflection of carriers at the surface [28].

In this section, we are interested in a direct comparison between the semi-classical and the full quantum approach. Strictly speaking, at equilibrium there is no net flow of current. Therefore, in the semi-classical theory of the metal we define the meaning of current density at point $\vec{r}$ and time t as:

$$\begin{eqnarray}&&\vec{j}(\vec{r},t)=q\rho (\vec{r},t){{\vec{v}}_{g}}(\vec{r},t)\end{eqnarray} \tag{ 19 }$$

with a density of electrons per volume unit, ρ, carrying an effective charge q, and moving at group velocity ${{\vec{v}}_{g}}(\vec{r},t)$. As the group velocity of electrons is proportional to the reciprocal-space gradient of their energies, after space and time Fourier transforms, we obtain:

$$\begin{eqnarray}&&{{\vec{j}}_{{{{\vec{k}}}_{//}}}}\propto \underset{n}{\sum}{\int}^{}{{\vec{\nabla}}_{{\vec{k}}}}\varepsilon _{{\vec{k}}}^{n}\cdot \delta ({{\varepsilon}^{n}}-\varepsilon _{{\vec{k}}}^{n})\,\text{d}{{\vec{k}}_{z}}\end{eqnarray} \tag{ 20 }$$

where $\delta ({{\varepsilon}^{n}}-\varepsilon _{{\vec{k}}}^{n})$ is the density of states at energy ε, $\varepsilon _{{\vec{k}}}^{n}$ are the eigenvalues of the $n\text{th}$ band, ${{\vec{k}}_{z}}$ is the component of the wave-vector parallel to the epitaxy direction and ${{\overset{{}}{{k}}}_{//}}$ is the component parallel to the interfaces (orthogonal to the epitaxy direction).

Therefore, in this approximation the current can be obtained from a simple bulk band-structure calculation. In particular, we have evaluated this expression in the tight-binding approach, using the same Hamiltonian as the one that was employed for non-equilibrium calculations. It is clear that this integral is zero for electrons which propagate in all directions. For such a current to be meaningful, we have to consider only electrons that propagate in the k_z  >  0 direction. This point of view has the advantage of a simpler physical interpretation with respect to the full non-equilibrium Keldysh approach of section 2 [6]. Of course, the choice k_z  >  0 states is artificial, but it can be justified it a posteriori by the final discussion of the next section, where the time-reversal breaking due to the loss of coherence is discussed. Notice the fact that choosing k_z  >  0 automatically implies 3-fold propagation symmetry around Au(1 1 1) direction.

Once the current-density vector is calculated for each $\vec{k}$-point, all ${{\vec{j}}_{{\vec{k}}}}$ have to be summed on the epitaxial components of the current-density. In other words, the full 3D Brillouin zone is projected into a 2D Brillouin Zone. Then the current is projected and it can be compared to the accessible density of states (DOS) in the semiconductor. We assume that the parallel component of the wave-vector is conserved at the interface; therefore, a current arriving at the last layer of the metal base (collector) crosses the interface only if it matches a corresponding available density of states in the semiconductor.

The comparison with the full non-equilibrium Keldysh calculation of [11], represented by the red curve in figure 5, is quite instructive. On one side, it is clear that the two results are different, which justifies the need for the full non-equilibrium Keldysh formalism. On the other side, the two results are similar; in fact, indistinguishable if the present experimental sensitivity is taken into account since only relative-size effects can be detected. Therefore, working with the much simpler semi-classical equilibrium procedure can be justified for many purposes (see, e.g. [13] or [6]); the conclusion being that most of the electrons propagate in narrow cones near $-{{24}^{\circ}}$ and ${{27}^{\circ}}$ with respect to the (1 1 1) direction, in the $\Gamma WWX \Gamma$ and $\Gamma KLUX \Gamma$ planes, respectively.

Zoom In Zoom Out Reset image size

This comparative analysis, however, needs a further specification. Both approaches assume that the BEEM current is represented by the overlap in the projected 2D Brillouin zone between regions of large currents in the last layer of the injector, n, and regions of available density of states in the first layer of the collector, n  +  1:

$$\begin{eqnarray}&&{{J}_{\text{B}}}(\text{eV})\propto {\int}^{}\text{d}{{\vec{k}}_{//}}{\int}_{{{\phi}_{\text{SB}}}}^{eV}\,{{J}_{n}}(E,{{\vec{k}}_{//}})\,{{T}_{n,n+1}}(E,{{\vec{k}}_{//}})\,{{\rho}_{n+1,n+1}}(E,{{\vec{k}}_{//}})\,\text{d}E\end{eqnarray} \tag{ 21 }$$

Here ${{T}_{n,n+1}}(E,{{\vec{k}}_{//}})$ is a transmission coefficient that describes the propagation from the last layer of the injector, usually a metal, to the first layer of the collector, either a different metal or a semiconductor. This transmission coefficient can be computed in a variety of ways, e.g. using a nearly-free electron model to match states in the metal to the semiconductor using the Jones-zone approximation [19].

On the other hand, assuming that the BEEM current, computed by simply summing up the elastic current enclosed by the surface DOS, and the transmission coefficient do not significantly depend on energy we may write:

$$\begin{eqnarray}&&{{J}_{\text{B}}}(\text{eV})\propto {\int}^{}\text{d}{{\vec{k}}_{//}}{{J}_{n}}(\text{eV},{{\vec{k}}_{//}}){{T}_{n,n+1}}(\text{eV},{{\vec{k}}_{//}}){\int}_{{{\phi}_{\text{SB}}}}^{\text{eV}}{{\rho}_{n+1,n+1}}(E,{{\vec{k}}_{//}})\,\text{d}E\end{eqnarray} \tag{ 22 }$$

Such a hypothesis cannot a priori be taken for granted and it must be checked before simply projecting the DOS. Figure 6 proves that for thick layers the distribution of the elastic electrons does not substantially vary with the energy. On the contrary, we have found that in thin films the dependence on the energy cannot be neglected. That is why, for structures of few atomic layers, a quantitative analysis of the BEEM current should pass through the evaluation of the integral equation (21). Working with the simpler projection scheme of the equilibrium calculation can at best lead to qualitative conclusions.

Zoom In Zoom Out Reset image size

4.2. Time-reversal symmetry: coherence

To summarize, in this paper:

• (1)  
  We have introduced a new calculation scheme for the BEEM current based on a layer-by-layer construction (Dyson's equation) and a Keldysh's non-equilibrium Green's function formalism. This scheme has been motivated by the availability of experimental data on ultra-thin samples, like those reported in [6, 7].
• (2)  
  We have implemented the calculation scheme in a new BEEM computer code that will be published elsewhere. The new code has been benchmarked against the old one for cases where similar results are expected.
• (3)  
  We have highlighted new features expected for very few atomic-layer systems, in particular quantum-size effects.
• (4)  
  Finally, we have studied the conditions that would allow using the simpler semi-classical approach successfully instead of the more elaborated fully quantum mechanical one.

This work has been supported by the MINECO (MAT-2014/54231-C4-1-P and MAT-2014/59966/R), the EU (ERC-2013-SYG-610256), and the program Maria de Maeztu (MDM/2014/0377).