General atomistic approach for modeling metal-semiconductor interfaces using density functional theory and non-equilibrium Green's function

3 MB·

by  Daniele Stradi

APS/123-QED
General atomistic approach for modeling metal-semiconductor interfaces using density
functional theory and non-equilibrium Green’s function
Daniele Stradi,1,2,∗ Umberto Martinez,2 Anders Blom,2 Mads Brandbyge,1 and Kurt Stokbro2
1Center for Nanostructured Graphene (CNG), Department of Micro- and Nanotechnology (DTU Nanotech),
Ørsteds Plads, Building 345B, DK-2800 Kongens Lyngby, Denmark
2Quantumwise A/S, Fruebjergvej 3, Postbox 4, DK-2100 Copenhagen, Denmark
Metal-semiconductor contacts are a pillar of modern semiconductor technology. Historically,
6 their microscopic understanding has been hampered by the inability of traditional analytical and
1 numerical methods to fully capture the complex physics governing their operating principles. Here
0 weintroduceanatomisticapproachbasedondensityfunctionaltheoryandnon-equilibriumGreen’s
2 function,whichincludesalltherelevantingredientsrequiredtomodelrealisticmetal-semiconductor
r interfaces and allows for a direct comparison between theory and experiments via I-Vbias curves
a simulations. We apply this method to characterize an Ag/Si interface relevant for photovoltaic
M
applications and study therectifying-to-Ohmictransition as function of thesemiconductor doping.
We also demonstrate that the standard “Activation Energy” method for the analysis of I-V
bias
1 data might be inaccurate for non-ideal interfaces as it neglects electron tunneling, and that finite-
size atomistic models have problems in describing these interfaces in the presence of doping, due

] to a poor representation of space-charge effects. Conversely, the present method deals effectively
i
c with both issues, thus representing a valid alternative to conventional procedures for the accurate
s characterization of metal-semiconductor interfaces.
-
l
r PACSnumbers: 72.10.-d,73.30.+y,73.40.-c,71.15.Mb
t
m
t. I. INTRODUCTION atomistic aspect of the interface, although it is nowa-
a
days accepted that chemistry plays a dominant role in
m
Metal-semiconductor(M-S)contactsplayapivotalrole determiningtheelectroniccharacteristicsoftheinterface
d- inalmostanysemiconductor-basedtechnology. Theyare [7, 9, 11–13]. These ambiguities complicate the assign-
n anintegralpartofabroadrangeofdeviceswithapplica- mentofthe featuresobservedinthe measuredspectrato
o tion as diverse as photovoltaics (PV) [1], transistors and specific characteristics of the M-S interface.
c diodes [2, 3], and fuel-cells [4, 5]. Conversely,atomisticelectronicstructuremethods[14]
[
The requirement of M-S interfaces with tailored char- are an ideal tool for the characterization of M-S inter-

3 acteristics, such as a specific resistance at the contact, faces,andhavebeensuccessfullyemployedovertheyears
v has fueled research on the topic for decades [6, 7]. Nev- for their analysis [15–22]. However, due to their compu-
1 ertheless,despitethehighdegreeofsophisticationofcur- tational cost, these studies have focussed on model in-
5
rentsemiconductortechnology,theunderstandingofM-S terfaces described using finite-size models formed by few
6
interfaces at the microscopic levelstill constitutes a con- atomiclayers(e.g.,slabs),thevalidityofwhichisjustified
4
0 siderablechallenge[8,9]. Eventhestructureoftheinter- intermsofthelocalnatureoftheelectronicperturbation
. faceitself, whichisburiedinthe macroscopicbulkmetal due to the interface. For similar reasons, most studies
1
and semiconductor materials, represent a serious imped- have considered non-doped interfaces, as the models re-
0
6 iment, as it makes the direct explorationof the interface quired to describe a statistically meaningful distribution
1 properties cumbersome. ofdopantsinthesemiconductorwouldbeexcessivelyde-
: A measure of the device current I as a function of the manding [23, 24]. Last but not least, these model calcu-
v
i applied bias Vbias is a standardprocedure to probe a M- lations only describe the system at equilibrium (i.e., at
X S interface [2], despite the drawback that the measured V = 0 V), thereby missing a direct connection with
bias
r I-Vbias curves do not provide any direct information on the I-Vbias measurements.
a
the interface itself, but rather on the full device charac-
Here, we developa general frameworkwhich attempts
teristics. As such, it is common practice to interpret the
to overcome the limitations inherent in conventional
I-V curves by fitting the data with analytical mod-
bias electronic structure methods for simulating M-S inter-
els, which are then used to extract the relevant interface
faces. We employ density functional theory (DFT)
parameters such as the Schottky barrier height Φ [10].
[25] together with the non-equilibrium Green’s func-
As no general analytical model exists, the accuracy of
tion (NEGF) method [26] to describe the infinite, non-
this procedure critically depends on whether the model
periodic interface exactly. The DFT+NEGF scheme al-
describes well the physical regime of the interface un-
lows us to predict the behavior of the M-S interface un-
der scrutiny. Furthermore, most models disregard the
der working conditions by simulating the I-V charac-
bias
teristics of the interface at zero and at finite V . To
bias
describe correctly the electronic structure of the doped
∗ [e-mail] semiconductor, we employ an exchange-correlation (xc)
2
functional designed ad-hoc to reproduce the experimen- and xc contributions of the KS Hamiltonian has been
tal semiconductor band gap [27], and a novel spatially set to 150 Ry. Monkhorst-Pack [43] grids of k-points
dependent effective scheme to account for the doping on have been used to sample the 3D (2D) Brillouin zone in
the semiconductor side. the DFT (DFT+NEGF) simulations. We have used an
We apply this novel DFT+NEGF approach to study 11×11×11 grid of k-points for the bulk calculations,
the characteristics of a Ag/Si interface relevant for PV and a k-points grid of 18×9×1 (18×9) for the DFT
applications [12, 23, 28–38]. Specifically, we focus on (DFT+NEGF) simulations of the interface. Geometry
the (100)/(100)interface [28, 37] and on the dependence optimizations have been performed by setting the con-
of I-V characteristics on the semiconductor doping – vergence threshold for the forces of the moving atoms to
bias
notice that the method is completely general and can 2 × 10−2 eV/˚A. In all the simulations, periodic bound-
be used to describe other M-S interfaces with different aryconditions(PBCs)wereusedtodescribetheperiodic
crystalline orientations. We consider a range of doping structure extending along the directions parallel to the
densities for which the interface changes from rectifying interface plane. In the slab DFT simulation, Dirichelet
toOhmic. Wedemonstratethatthe“ActivationEnergy” and Neumann boundary conditions were applied in the
(AE)methodroutinelyemployedtoanalyseM-Scontacts direction normal to the interface on the silver and sili-
systematically overestimate the value of Φ, with an er- con sides of the simulation cell, respectively, whereas in
ror that is both bias and doping dependent, due to the the DFT+NEGF simulations the same direction wasde-
assumption of a purely thermionic transport mechanism scribed using Dirichelet boundary conditions at the two
across the barrier. Conversely, we show how an analysis boundaries between the interface and the bulk-like elec-
of the I-V characteristics based on the DFT+NEGF trodes.
bias
electronic structure data provides a coherent picture of
the rectifying-to-Ohmic transition as the doping is var-
ied. Finally, we also show that a slab model does not A. “Spill-in” terms
provide a good representation of the interface electronic
structurewhendopinginthesemiconductoristakeninto AsdescribedinRef. 26,theDFT+NEGFmethodused
account. Thisisduetotheinabilityofthesemiconductor tosimulatetheinfinite,non-periodicAg(100)/Si(100)in-
side of the slab to screen the electric field resulting from terface relies on a two-probe setup, in which a left (L)
the formation of the interface. and a right (R) semi-infinite electron reservoirs are con-
The paper is organized as follows. Section II and Sec- nectedthroughacentral(C)regioncontainingthe inter-
tionIIIdescribethecomputationalmethodsandthesys- face. Once the chemicalpotentials µ of the reservoirs
L,R
tem models employed in this work, respectively. Section have been defined, a self-consistent (SCF) KS procedure
IVA presents the calculated I-Vbias characteristics and is used to obtain the electronic density in the C region.
the validationofthe AE methodbasedonthe calculated The main quantity being evaluated in the SCF cycle is
data. Section IVB deals with the analysis of the I-Vbias thedensitymatrixrequiredtoexpresstheelectronicden-
curves in terms of the electronic structure of the inter- sity of the C region in the basis of PAOs centered in the
face as obtained from the DFT+NEGF calculations. In same region, D¯CC. Assuming µ > µ , D¯CC takes the
L R
Section IVC, the simulations are compared to finite-size form
slab models. The main conclusions are drawn in Section
V.
1 µR
D¯CC =− Im[G¯CC]dǫ
π
Z−∞ (1)
II. COMPUTATIONAL METHODS 1 µL
− G¯CCIm[Σ¯LL]G¯†CCdǫ,
π
ZµR
The Ag(100)/Si(100) interface has been simulated us-
ingKohn-Sham(KS)DFTasimplementedinatomistix where Σ¯LL is the self-energy matrix describing the cou-
toolkit [39] (ATK). DFT [25] and DFT+NEGF [26] plingofthecentralregiontothesemi-infiniteLreservoir,
simulationshavebeenperformedusingaformalismbased andtheGreen’sfunctionofthecentralregionG¯CC isob-
on a non-orthogonal pseudo-atomic orbitals [40] (PAOs) tained by
basis set.
The one-electronKSvalenceorbitalsareexpandedus-
ing a linear combinationof double-ζ PAOs including po- G¯CC(ǫ)=[(ǫ+iδ)S¯CC −H¯CC −Σ¯LL−Σ¯RR]−1, (2)
larizationfunctions(DZP).Theconfinementradiir em-
c
ployed are 4.39 Bohr, 7.16 Bohr, 7.16 Bohr for the Ag with S¯CC and H¯CC being the overlap and Hamiltonian
4d, 5s and 5p orbitals, and 5.40 Bohr, 6.83 Bohr, 6.83 matrices associated with the PAOs centered at the C
Bohr for the Si 3s, 3p, 3d orbitals, respectively. The region, and Σ¯RR being the self-energy matrix of the R
ionic cores have been described using Troullier-Martins reservoir.
[41] norm-conserving pseudo-potentials [42]. The energy However, even if the DFT+NEGF method provides
cutoffforthereal-spacegridusedtoevaluatetheHartree anelegantscheme to evaluate D¯CC, it shouldbe noticed
3
“spill-in”. For the Hamiltonian, corrective terms are ap-
plied to both the two-center and three-center integrals.
Specifically, if two PAOs φ and φ centered in the C
i j
regionlie close to a boundary, e.g. the L/C one, the cor-
rected Hamiltonian in Eq. 2 will include also the matrix
element H′ = hφ |VLL(r)|φ i associated with the tail
i,j i eff j
of the PAOs protruding into the L region, VLL(r) being
eff
the periodic KS potential of the semi-infinite L reservoir
(Fig. 1a). Similar arguments hold also for the Hamil-
tonian three-center non-local terms. For the electronic
density of the C region, additional contributions are in-
cluded for each pair of PAOs φ and φ located close to
i j
a boundary when at least one of them is centered at the
neighboringreservoirregion. In total, twonew contribu-
tions must be added to the electronic density evaluated
using Eqs. 1-2 for each pair of PAOs at each boundary.
For the L/C boundary, these are (Fig. 1b):
nLL = DLLφLφL,
i,j i j
i,j
X (3)
nLC = DLCφLφC,
i,j i j
i,j
X
which can be further distinguished based on whether
both (D¯LL) or just one (D¯LC) of the twoPAOs involved
is centeredatthe L region. Inthe calculationspresented
in this work, the “spill-in” terms are independent of the
applied bias voltage. This is justified as we checked
that the non-periodic KS potential at the boundary of
the C region for each value of the applied bias matches
smoothly with the periodic KS potential of the neigh-
boringreservoir,i.e. thatthe“screeningapproximation”
is verified – see Ref. 26 for additional details. We stress
thatincludingthese“spill-in”termsisessentialtoensure
a stable and well-behaved convergence behavior of the
SCF cycle, which turns out to be especially important
for heterogeneous systems such as the Ag(100)/Si(100)
interface investigated in this work.
FIG.1. Schemeshowingthetwo-center“spill-in”termsused
for the evaluation of the Hamiltonian (a) and the electronic
density (b) of the C region for a pair of s type PAOs φ
j
and φ located close to the L/C boundary. In (a,b) the blue B. Exchange-correlation potential
i
shaded regions indicate the integrals performed in the C re-
gion,whereastheredandgreenregionsindicatethe“spill-in” Further complications in describing the
terms for the Hamiltonian and for the electronic density, re-
Ag(100)/Si(100) interface arise from the fact that
spectively.
one of its sides is semiconducting. In fact, a majorprob-
lem affecting the description of metal-semiconductor
interfaces is the severe underestimation of the semicon-
that solving Eqs. 1-2 is not sufficient to obtain the cor- ductor band gap in DFT calculations using (semi-)local
rect Hamiltonian and the electronic density of the C re- xc-functionals based on the local density approximation
gion. The reason is that the relevant integrals involved (LDA) or on the generalized gradient approximation
in Eqs. 1-2 are evaluated only in the region of space en- (GGA) [45]. For model calculations based on few-layer
compassingtheCregion,andonlyfortheatomslocalized thick fully periodic systems, such an underestimation
in that region. As a consequence, the tails of the PAOs has been shown to result in unrealistically low Schottky
located close to both sides of the L/C and R/C bound- barriersat the interface [15, 46]. In orderto remedy this
aries, which penetrate into the neighboring regions, are drawback, we have evaluated the electronic structure
not accounted for (see Fig. 1). To correct this behavior, of the LDA-optimized interface geometries using the
we introduce additional corrective terms, that we name Tran-Blaha meta-GGA xc-functional (TB09) [27]. The
4
a
where Ω is the volume ofthe simulation celland the two
empirical parameters α = −0.012 (dimensionless) and
β = 1.023 Bohr21 have been fitted to reproduce the ex-
perimental band gaps of a large set of semiconductors
[27]. To obtain a description as accurate as possible of
the semiconductor band gap at the Si(100) side of the
interface, we have tuned the value of the c parameter
in orderto reproduce the experimentally measuredband
gap of bulk silicon, Eexp = 1.17 eV [44]. This has been
gap
accomplishedby calculating the bandgapof bulk silicon
at fixed values of the c parameter in a range around the
b self-consistently computed value in which the variation
of ETB09 with c is linear. Then, the optimal value of c
gap
has been determined as the intersect between the value
of Eexp and a linear fit to the computed values of ETB09
gap gap
EgTaBp0 9 - o = 1.169 eV (Fig. 2a). Using the TB09 xc-functional with the c pa-
rameterfixed atthe optimal value determined using this
procedure (hereafter, TB09-o), we calculate the indirect
bandgapofbulk siliconto be ETB09−o = 1.169eV (Fig.
gap
2b), in excellent agreement with the value 1.17 eV. The
TB09-o functional has been used for all the electronic
structure and transport analyses of the Ag(100)/Si(100)
interfacereportedinthiswork. Wehavecheckedthatthe
band structure of bulk silver calculated using the TB09-
FIG. 2. (a) Fitting procedure for the TB09 xc-functional o functional is very similar to that calculated using the
c parameter. Squares (blue): calculated indirect band gap LDA, which is known to perform well for noble metals.
of bulk silicon ETB09 obtained for different values of the c
gap
parameter. Dashed line (blue): fit to the computed data of
ETgaBp09vs. cobtainedbylinearregression. Dottedline(black): C. Semiconductor doping
experimentallymeasuredbulksilicon bandgap[44]. Dashed-
dotted line (orange): optimal value of the c parameter (opt-
The last requirement to describe realistically the elec-
c), obtained as the intersect between ETB09 and the fit to
gap tronic structure of the Si(100)/Ag(100) interface is to
the ETB09 data. (b) Region around the indirect band gap in
gap account for the doping on the silicon side of the inter-
the bulk silicon band structure calculated using the optimal
face. Here, doping is achieved in an effective scheme by
c parameter determined from (a).
introducinglocalizedchargesboundtotheindividualsil-
icon atoms. More specifically, in ATK [39] the total self-
consistent electronic density ρ (r) is defined as [40]:
TB09 xc-functional has been shown to provide band tot
gaps in excellent agreement with the experiments for
a wide range of semiconductors including silicon, at a Natoms
computational cost comparable to that of conventional ρ (r)=δρ(r)+ ρ (r), (6)
tot I
(semi-)local functionals. In the TB09 xc-functional,
I
the exchange potential υTB09(r) depends explicitly on X
x
electron kinetic energy τ(r), where Natomsρ (r) is the sum of the atomic densities
I I
of the individual neutral atoms of the system. As each
atomicPdensity ρ (r) is a constant term, it can be aug-
I
3c−2 4τ(r) mented with a localized “compensation” charge having
υxTB09(r)=cυxBR(r)+ π s6ρ(r), (4) the opposite sign of the desired doping density, which
acts as a carrier attractor by modifying the electrostatic
potential on the atom. Using these “compensation”
with τ(r)=1/2 Ni=1|∇ψi(r)|2, N being the total num- charges, an effective doping can be achieved both in the
ber of KS orbitals, ψi(r) the i-th orbital, ρ(r) the elec- DFTandintheDFT+NEGFsimulations. Intheformer,
tronic density aPnd υxBR(r) the Becke-Roussel exchange the “compensation”chargeaddedto eachsiliconatomis
potential [47]. The parameter c in equation (4) is evalu- neutralized by explicitly adding a valence charge of the
ated self-consistently and takes the form oppositesign,sothatthesystemremainschargeneutral.
In the latter, the “compensation” charge is neutralized
implicitly by the carriersprovidedby the reservoirs,and
1 |∇ρ(r)| 21 thesystemismaintainedchargeneutralunderthecondi-
c=α+β dr , (5)
Ω ρ(r) tion that the intrinsic electric field in the system is zero.
(cid:20) Z∞ (cid:21)
5
a
Left (μL) Right (μR)
W
Si(100)
b
te nn
lhciriD.c.b amueN.c.b
Wslab
Si(100)
FIG.3. GeometriesemployedtosimulatetheAg(100)/Si(100) interfaceusingtwo-probemodels(a) orslab models(b). Silver,
silicon and hydrogen atoms are shown in grey,beige and white, respectively.
Thiseffectivedopingschemehastheadvantageof(i)not the Si(100) dangling bonds sitting above the “hollow”
depending on the precise atomistic details of the doping fcc sites of the Ag(100) surface, has then been used as a
impurities, and (ii) being completely independent of the representative model of the interface.
size and exact geometry of the system. Startingfromthelowestenergyconfigurationobtained
using the 15-layers slab, we have then constructed more
realistic models of the interface. Specifically, we have
expanded the bulk regions of the 15-layer slab to cre-
III. SYSTEM
ate two-probe setups effectively describing the infinite,
non-periodic interface (Fig. 3a). A final geometry opti-
In order to obtain a reliable description of the mization has been carried out using a two-probe setup
Ag(100)/Si(100) interface, we have followed a stepwise in which the C region has been described by 8 Ag(100)
procedure. Initially, we have carried out a preliminary layers and an undoped silicon layer having a total width
screeningoftheinterfacegeometriesandbondingconfig- WCC = 47.84 ˚A. The optimized geometry has been
urations by using a 2×1 slab model formed by a 6-layers Si(100)
used to construct two-probe setups in which the doping
Ag(100) slab interfaced with a 9-layers unreconstructed
of the silicon side has been taken into account using the
Si(100)slab. The calculatedbulk lattice constantsofsil-
effective doping method described in Section IIC. We
icon(a = 5.41˚A)andsilver (a = 4.15˚A)are ingood
Si Ag have considered doping densities of n-type carriers (n )
d
agreement with those reported in the literature [23]. To in the experimentally relevant range [1018 cm−3 – 1020
match the Ag(100) and the Si(100) surface, we have ap- cm−3]. As discussed in more detail in Section IV, the
pliedanisotropiccompressivestrainǫ =ǫ =–0.0793
xx yy width of the Si(100) layer needed to describe accurately
along the surface lattice vectors v of the Ag(100) sur-
1,2 theinterfaceinthetwo-probesimulationsdependsonthe
face. We have checked that in the compressed Ag(100)
size of the depletion region (W ) on the silicon side of
D
structure, the dispersion of the s-band and its position
the interface. The relation between W and n is 1/W
D d D
with respect to the d-band are very similar to those cal-
∝ n1/2, so that progressively narrower C regions can be
culated using the equilibrium value of a . The Si(100) d
Ag
used as the doping level is increased without any loss in
surfaceoppositetotheinterfacehasbeenpassivatedwith
accuracy. Therefore, in the following, the results pre-
hydrogenatoms. Thegeometryoftheresulting15-layers
sented for n = 1020 cm−3, n = 1019 cm−3 and n =
slab has then been optimized using the LDA by keeping d d d
1018 cm−3 refertocalculationsperformedwithCregions
the farthest (with respect to the interface plane) 4 lay-
of widths WCC = 47.84 ˚A, WCC = 197.436 ˚A and
ers of the Ag(100) surface frozen, and by allowing the Si(100) Si(100)
farthest (with respect to the interface plane) 4 layers of WCC = 447.92˚A, respectively. We have checked that
Si(100)
the Si(100) slab to move as a rigid body, thereby freez- reducing the width of the C regiondoes not haveany ef-
ing only the interatomic distances and angles. All the fect on the results, as long as all the space-chargeeffects
remainingatomshavebeenallowedtofullyrelax. Differ- dueto the presenceofthe interfacetakeplacewithinthe
entstarting guessesfor the interface structure havebeen screening region. Furthermore, we notice how using a
tested, corresponding to different configurations of the two-probe setup also allows to simulate the characteris-
Si(100) dangling bonds with respect to the high symme- tics of the interface when the L and R reservoirs are set
try fcc sites of the Ag(100) surface. The lowest energy at two different chemical potentials µ 6= µ due to an
L R
configuration among those considered, corresponding to applied bias voltage qV = µ −µ . As will become
bias R L
6
a
IV. RESULTS
A. Device characteristics and validation of the
activation energy model
×1
×10
×100
Fig. 4a shows the current–voltage (I-V ) character-
bias
istics calculated for the two-probe setup at low (n =
d
1018 cm−3), intermediate (n = 1019 cm−3) and high
d
(n = 1020 cm−3) doping densities of the Si(100) side of
d
the interface. A strong dependence on the doping con-
centrationis evident. At low doping,the interface shows
a well-defined Schottky diode-like behavior: the forward
b
bias (V > 0 V) current increases about six orders of
bias
magnitudeintherangeofV [0.02V–0.5V],whereas
bias
n = 2.40 the reverse bias one (V < 0 V) varies only within
bias
one order of magnitude in the corresponding range. The
diode-likeasymmetryintheI-V curvespersistsatin-
bias
termediate doping, although it is less pronounced than
n = 1.82 at low doping; the current at forward bias and reverse
bias varying within three and two orders of magnitude,
respectively. The scenario changes qualitatively at high
doping as the I-V curve becomes highly symmetric,
bias
n = 1.09
suggesting an Ohmic behavior of the interface.
According to thermionic emission theory, the I-V
bias
characteristics of a Schottky diode can be described by
[2]
qVbias
I =I0 enkBT −1 , (7)
h i
where q is the elementary charge, k is the Boltzmann
FIG. 4. Calculated I-V (a) and forward bias I/(1 − B
bias constant, T is the temperature, I is the saturation cur-
eq|Vbias|/kBT)-Vbias (b) characteristics at nd = 1018 cm−3 rent and n is the so-called ideal0ity factor. The latter
(blue triangles), n = 1019 cm−3 (green squares), n = 1020
cm−3 (red dots). dIn (a), the values of I at n = 10d18 cm−3 accounts for the deviation of the I-Vbias characteristics
d
andn =1019 cm−3 havebeenmultiplied byafactor 10and from those of an ideal diode, for which n = 1. Fitting
d
100, respectively. The solid lines in (b) are fit to the data in the simulated data at forward bias to Eq. 7 allows to
therange0.02≤V ≤0.08VusingEq. 7. Theidealityfac- extract n from the slope of the fitted curves. In Fig. 4b
bias
tornextractedfromtheslopeofeachfittedcurveisreported the fitted curves are compared to the forward bias data.
using thesame color as thecorresponding curve. Thelatter arepresentedusinganalternativeformofEq.
7,
clear later, this allows for a direct comparison to exper-
qVbias −qVbias
iments and for analyzing the electronic structure of the I =I0enkBT 1−e kBT , (8)
interface under working conditions.
(cid:16) (cid:17)
Finally, to understand to which extent the slab model which allows for a better comparison with the fitted
is able to describe accurately the electronic structure of curves as I/(1−e−qVbias/kBT) varies exponentially with
the infinite, non-periodic interface, we have also consid- V in the fitting interval considered, viz. 0.02 V ≤
bias
ered slab models having a similar interface geometry as V ≤ 0.08 V.
bias
that used in the two probe setup (Fig. 3b). Both short At low doping, n = 1.09, indicating that the system
and long slab models have been constructed, in which behaves essentially as an ideal Schottky diode. At in-
thewidthoftheSi(100)layerusedtodescribethesilicon termediate doping, n = 1.82, and the system deviates
side of the interface has been set to either Wslab(short) = significantly from the ideal behavior. At high doping, n
Si(100)
=2.40,consistentlywiththeobservationthatthesystem
38.33 ˚A or Wslab(long) = 98.62 ˚A. Notice how these val-
Si(100) does not behave anymore as a Schottky diode.
uesofWslab aremanytimes largerthanthoseusedfor The I–V simulation allows to test the reliability of
Si(100) bias
similar studies of the Ag(100)/Si(100)interface reported the experimental procedures used to extract the Schot-
in the literature [37]. tky barrier Φ. In particular, we focus on the so-called
7
a
c n = 1019 cm-3 d
d
n = 1018 cm-3
d
n = 1018 cm-3
d
e f
b
nd = 1018 cm-3 nd = 1019 cm-3
n = 1019 cm-3
d
FIG. 5. (a,b) Empty dots: calculated I-T data at different bias voltages for n = 1018 cm−3 (a) and n = 1019 cm−3 (b).
d d
Solid lines: fit to the simulated data using Eq. 10. (c,d) Left-hand side (filled dots) and right-hand side (dashed line) of Eq.
10 as a function of V . Thevalues of the left-hand side havebeen extracted from the slope of thefitted I-T curvesin (a,b).
bias
Thesolid linesare linear fitstothedata. Theright-handside of Eq. 10 hasbeen plotted using thevalueof ΦAE calculated at
V =0.02 V,which approachesthevalueofΦat V = 0V.(e,f) SchottkybarrierheightΦAE evaluated usingEq. 10 asa
bias bias
function of V .
bias
“Activation-Energy” (AE) method, which does not re- self-consistentsimulationsperformedforselectedtemper-
quire any a priori assumption on the electrically active atures show that this approach is valid within the range
interface area A [2]. In the AE method the I–T depen- of T considered, 250 K ≤T ≤ 400 K.
dence is measuredat a small constantV . Overa lim-
bias Ideally, for a given doping the Schottky barrier de-
itedrangeofT aroundroomtemperature,assumingthat
pends exclusively on the M-S energy level alignment at
the RichardsonconstantA∗ andΦ are constant,the I-T
theinterfaceandtherefore,disregardingimage-forcelow-
characteristics can be described by the expression
ering effects, shouldremainconstantwith V [2]. This
bias
implies that in Eq. 10, the left-hand side should equal
IT−2 =AA∗e−qkΦBATE eq(VkbiBaTs/n). (9) tphreesreingthtc-ahsaendthsiisdecoantdaitniyonvailsuenootfVvbeiraisfi.edH,oawsevthere,ivnartihae-
tion of the right-handside term with V is largerthan
Following Eq. 9, the Schottky barrier height ΦAE can bias
that of the left-hand side term (see Fig. 5c,d). Indeed,
be extracted from the ln(I/T2) vs. 1/T data using for n = 1018 cm−3 (n = 1019 cm−3), a linear fit to the
d d
calculated values of the left-hand side of Eq. 10 gives a
k d[ln(I/T2)] V slope of –664 meV/V (–177 meV/V), whereas the slope
− B =ΦAE− bias, (10) associatedto the variationof the right-handside term is
q d(1/T)) n
–917 meV/V (–549 meV/V).
in which n is the ideality factor extracted above. Following the procedure in the AE method we use the
Fig. 5a,b shows the simulated AE plots (as Arrhenius value of n obtained from Fig. 4 to subtract the bias de-
plots) at different values of Vbias for low and intermedi- pendence. The result is shown in Fig. 5e,f and it can
ate doping densities, at which the interface still displays be seen that this leads to an unphysical increase of ΦAE
clear Schottky diode–like characteristics. The I-T de- with V . The error becomes more severe as n is in-
bias d
pendence has been evaluated in a linear response fash- creased. At low (intermediate) doping, ΦAE varies from
ion, using the Landauer-Bu¨ttiker expression for the cur- by 30% (325%) in the range of V considered, lead-
bias
rent,I = 2hq T(E,µL,µR)[f(Ek−BµTL)−f(Ek−BµTR)]dE with ing to a change∆ΦAE = 3.73 kBT (∆ΦAE = 7.31kBT).
the transmission coefficient T(E,µ ,µ ) evaluated self- Thus, the intrinsic accuracy of the AE method depends
R L R
consistently at an electron temperature of 300 K. Fully strongly on multiple factors. On the one side, the non-
8
a
n = 1018cm-3
d
b
n = 1019cm-3
d
c
n = 1020cm-3
d
FIG.6. Localdensityofstates(LDOS)ofthetwo-probesetupatequilibriumforn =1018 cm−3 (a),n =1019 cm−3 (b)and
d d
nd =1020 cm−3 (c). Theenergyontheverticalisrelativetothesystem chemicalpotentialµL,R. Regionsoflow(high)LDOS
are shown in dark (bright) color. The blue line in each panel indicates the macroscopic average of the Hartree potential hV i
H
subtracted the electron affinity of bulkSi and µL,R. The yellow vertical line in each panel indicates theassociated Φpot.
linear increase in ΦAE with V suggests that using a alongthedirectionnormaltotheinterfaceatthedifferent
bias
singlevalueofV isnotsufficienttoobtainanaccurate doping densities considered. Increasing the doping has a
bias
estimate of Φ. On the other side, the change in ∆ΦAE two-foldeffectonthe electronicpropertiesofthesystem:
with n atagivenV indicates thatthe AE method is onthe one side, W decreasesfrom∼200nmto ∼20nm
d bias D
unsuited for comparative analyses of the variation of Φ when the doping is increased from n = 1018 cm−3 to
d
with doping. These facts call for a more direct and gen- n =1020 cm−3,asadirectconsequenceoftheincreased
d
eral strategy for the characterization of M-S interfaces screeningofthen-dopedsilicon. Increasingn alsoshifts
d
under working conditions. the Fermi level towards the silicon conduction bands. In
particular,atn =1018 cm−3 theconductionbandmini-
d
mum(CBM)ofsiliconatZ>W liesatE−µ =+20
D L,R
meV,whereasatn =1019 cm−3 andn =1020 cm−3 it
B. Electronic properties of the interface d d
liesatE−µ =–40meVandE−µ =–100meV,re-
L,R L,R
spectively. Itis also worthnoticing how the macroscopic
A strongadvantageofthe DFT+NEGF simulations is
average [49] of the Hartree potential along the direction
that they allow the visualization of the electronic struc-
normal to the interface, hV i (blue lines in Fig. 6), fol-
H
tureoftheinterfaceandthedirecttrackingofitschanges
lowsthe profileof the siliconCBM close to as well as far
when n and V are varied. This makes it possible to
d bias away from the interface. Similarly to what happens for
analyze the calculated I-V characteristics in terms of
bias theelectronicbands,hV ibecomesconstantatZ>W ,
H D
the electronic structure of the interface.
indicatingthattheelectronicstructurestartstoresemble
Fig. 6showsthe localdensityofstates [48](LDOS)of that of the infinite periodic bulk. A closer analysis also
thetwo-probemodelatequilibrium(i.e., atV =0V)
bias
9
aa
a
b c
n = 1018 cm-3
d
n = 1019 cm-3 b
d
n = 1018 cm-3 n = 1018 cm-3
d d
d n = 1019 cm-3 e n = 1019 cm-3
d d
FIG.7. (a)SchemeoftheelectronicstructureoftheAg/Siin-
terfaceatforwardbiasvoltage. (b)ProfileofhV ifordifferent
H
V atlowdoping. Theenergyontheverticalaxisisrelative FIG. 8. (a) Filled circles: energy of maximum spectral cur-
bias
to the electron affinity χ of bulk Si and the metal chemical rent E(Imax) in Fig. 7c as a function of Vbias at low doping.
potentialµL. TheverticallinesindicateφF atVbias =0.02V Thesolidlineisaguidetotheeyes. Filled squares: variation
(blue,solid)andV =0.2V(red,dashed). (c)Solidcurves: of the slope-dependent term of Eq. 10 (same as in Fig. 5c).
bias
spectral current density I(E) for different Vbias at low dop- Thesolid lineisaguidetotheeyes. Filled triangles: φF asa
ing. The dashed line indicates the value of Φpot. hVHi and functionofVbias. Thedashedlineshowsthebiasdependence
I(E)curvescalculated atincreasingly higherVbias areshown Vbias/n from Eq. 10. The energy on the vertical axis is rel-
in blue→green→yellow→red color scale. (d,e) Same as (b,c), ative to the semiconductor chemical potential µR. (b) Same
but for intermediate doping. as (a), but for intermediate doping.
reveals that a finite density of states extends consider- For nd = 1020 cm−3 the barrier is considerably lower,
ably on the semiconductor side of the interface, due to Φpot = 133 meV, reflecting the more pronouncedOhmic
penetrationofthe metallicstatesintothesemiconductor behaviorobservedin the I-Vbias curves. Focusing on the
side [50–52]. lowandintermediatedopingcases,itcanbenoticedhow
the values of Φpot are considerably larger than those of
Due to the lack of a well-defined electronic separation
ΦAE at V → 0 V. In particular,atlow doping Φpot−
between the metal and the semiconductor, it is difficult bias
ΦAE = 112 meV, whereas at intermediate doping the
toprovideanunambiguousvalueforΦbasedontheelec-
difference is even larger, Φpot−ΦAE = 286 meV.
tronicstructuredataonly. However,dueto the factthat
hV i closely traces the CBM, it is still possible to esti- A consistent physical picture that rationalizes the
H
mate the Schottky barrier by defining Φpot as the dif- I-V curves can be obtained by studying the dop-
bias
ference between µ and the maximum of hV i on the ing dependence of the spectral current I(E) =
L H
semiconductor side of the interface, hVmaxi (see Fig. 6). 2qT(E,µ ,µ )[f(E−µL)−f(E−µR)]. Fig. 7b,d shows
H h L R kBT kBT
We calculate Φpot = 412 meV and Φpot = 342 meV the profiles of hV i obtained at forward bias in the bias
H
for n = 1018 cm−3 and n = 1019 cm−3, respectively. range 0.02 V < V < 0.2 V for low and intermediate
d d bias
10
doping densities. As V is increased, hVmaxi shifts to- a
bias H
wards higher energies due to image-force effects [2], and
hV ibecomes progressivelyflatter onthe semiconductor
H
side. The overall result of these changes is a decrease of
the barrier φ associated with the thermionic emission
F
processfromtheSi(100)conductionbandtoAg(100)(see
Fig. 7a):
φ =Φ−V /n. (11) n = 1018 cm-3
F bias d
The associated spectral currents I(E) are shown in Fig. b
7c,e. For an interface in which the only contribution to
transport comes from thermionic emission, I(E) should
be non-zero only at E − µ > Φpot. However, in the
L
presentcaseI(E)isfinitealsoatE−µ <Φpot,indicat-
L
ingthatelectrontunnelinghasanon-negligiblecontribu-
tion to I. This contribution is much larger for interme-
diate than for low doping densities. Indeed, at V →
bias n = 1019 cm-3
0 V, the position of E(Imax) lies very close to Φpot in d
the low doping case, as expected in the case of a nearly c
ideal Schottky diode. Conversely, at intermediate dop-
ingE(Imax)lieswellbelowΦpot,indicatingthatelectron
tunneling has become the dominant transport process.
The trend of I(E) with V is consistent with these
bias
considerations. At low doping, E(Imax) is pinned to
hVmaxi,whereastheonsetoffiniteI(E)atE−µ <Φpot
H L
movestowardshigherenergies,followingthe variationof
n = 1020 cm-3
hVHi. On the other hand, at intermediate doping the d
overall shape of I(E) remains the same as V is in-
bias
creased,andthevariationofE(Imax)followscloselythat
of hV i. We also notice the presence of a narrow reso-
H
nanceatE−µ =+0.395eV,whosepositionisindepen-
L
dentonn andV . Thisisduetoalocalizedelectronic
d bias FIG. 9. Profile of hV i along the direction Z normal to the
state at the interface which is pinned to µ . H
L interfaceplane,calculatedforthetwo-probesetup(bluesolid
The variation of φ with V can be related to the
F bias line)andfortheshort(greendottedline)andlong(reddashed
slope-dependent term of Eq. 10 by assuming Φ = ΦAE line)slabmodels,atdopingdensitiesn =1018 cm−3 (a),n
d d
inEqs. 10-11,thusallowingforadirectcomparisonwith =1019 cm−3 (b)andn =1020 cm−3 (c). Theverticalblack
d
the AE data (see Fig. 8). Independently of the value solid lineindicated theposition of theinterface. Thevertical
of V , the slope-dependent term lies always below φ , green (red) line indicates the position Si(100) layer farthest
bias F
due to the missing contribution of electron tunneling in from the interface in theshort (long) slab model.
the AE method: the latter assumes that the current has
a purely thermionic origin, and consequently predicts a
value of φ lower than the actual one. In agreement often hV i [49]. The perturbation of the bulk electronic
F H
withthepreviousanalyses,thisdeviationisconsiderably structure in each material due to the presence of the in-
larger in the intermediate doping case, due to the much terface is accounted for by either a slab [37] or a fully
larger contribution of electron tunneling. periodic [57] model. hV i is then used as a common
H
reference to align the electronic structure obtained from
independent calculations of the two bulk materials. De-
C. Comparison of the two-probe with the slab spite its widespreaduse, this strategyrelies on two dras-
model ticassumptions. Firstly,itisimplicitlyassumedthatthe
electronic properties of the interface are independent of
The results obtained using the two-probe model can thedopinglevelofthesemiconductor. Moreover,itisas-
be used as a reference to validate the use of finite-size sumedthattheelectronicpropertiesinthecentralpartof
models to describe the Ag(100)/Si(100) interface. Such eachsideoftheinterfacemodelareagoodapproximation
models are integral parts of the band alignment method to those of the two bulk materials.
oftenusedtoevaluateΦusingconventionalDFT[53–56]. Fig. 9 shows a comparison between hV i obtained at
H
The method relies on aligning the electronic band struc- the different doping densities consideredfor the two slab
tures of the two bulk materials forming the interface on models (short andlong, see Section III) and for the two-
an absolute energy scale by using a reference quantity, probesetup. Wenoticethatintroducinganeffectivedop-