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

Metal-semiconductor contacts are a pillar of modern semiconductor technology. Historically, their microscopic understanding has been hampered by the inability of traditional analytical and numerical methods to fully capture the complex physics governing their operating principles. Here we introduce an atomistic approach based on density functional theory and non-equilibrium Green's function, which includes all the relevant ingredients required to model realistic metal-semiconductor interfaces and allows for a direct comparison between theory and experiments via I-V bias curves simulations. We apply this method to characterize an Ag/Si interface relevant for photovoltaic applications and study the rectifying-to-Ohmic transition as function of the semiconductor doping.We also demonstrate that the standard"Activation Energy"method for the analysis of I-V bias 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 with both issues, thus representing a valid alternative to conventional procedures for the accurate characterization of metal-semiconductor interfaces.


I. INTRODUCTION
Metal-semiconductor (M-S) contacts play a pivotal role in almost any semiconductor-based technology.They are an integral part of a broad range of devices with application as diverse as photovoltaics (PV) [1], transistors and diodes [2,3], and fuel-cells [4,5].
The requirement of M-S interfaces with tailored characteristics, such as a specific resistance at the contact, has fueled research on the topic for decades [6,7].Nevertheless, despite the high degree of sophistication of current semiconductor technology, the understanding of M-S interfaces at the microscopic level still constitutes a considerable challenge [8,9].Even the structure of the interface itself, which is buried in the macroscopic bulk metal and semiconductor materials, represent a serious impediment, as it makes the direct exploration of the interface properties cumbersome.
A measure of the device current I as a function of the applied bias V bias is a standard procedure to probe a M-S interface [2], despite the drawback that the measured I-V bias curves do not provide any direct information on the interface itself, but rather on the full device characteristics.As such, it is common practice to interpret the I-V bias curves by fitting the data with analytical models, which are then used to extract the relevant interface parameters such as the Schottky barrier height Φ [10].As no general analytical model exists, the accuracy of this procedure critically depends on whether the model describes well the physical regime of the interface under scrutiny.Furthermore, most models disregard the * daniele.stradi@quantumwise.com atomistic aspect of the interface, although it is nowadays accepted that chemistry plays a dominant role in determining the electronic characteristics of the interface [7,9,[11][12][13].These ambiguities complicate the assignment of the features observed in the measured spectra to specific characteristics of the M-S interface.
Conversely, atomistic electronic structure methods [14] are an ideal tool for the characterization of M-S interfaces, and have been successfully employed over the years for their analysis [15][16][17][18][19][20][21][22].However, due to their computational cost, these studies have focussed on model interfaces described using finite-size models formed by few atomic layers (e.g., slabs), the validity of which is justified in terms of the local nature of the electronic perturbation due to the interface.For similar reasons, most studies have considered non-doped interfaces, as the models required to describe a statistically meaningful distribution of dopants in the semiconductor would be excessively demanding [23,24].Last but not least, these model calculations only describe the system at equilibrium (i.e., at V bias = 0 V), thereby missing a direct connection with the I-V bias measurements.
Here, we develop a general framework which attempts to overcome the limitations inherent in conventional electronic structure methods for simulating M-S interfaces.We employ density functional theory (DFT) [25] together with the non-equilibrium Green's function (NEGF) method [26] to describe the infinite, nonperiodic interface exactly.The DFT+NEGF scheme allows us to predict the behavior of the M-S interface under working conditions by simulating the I-V bias characteristics of the interface at zero and at finite V bias .To describe correctly the electronic structure of the doped semiconductor, we employ an exchange-correlation (xc) functional designed ad-hoc to reproduce the experimental semiconductor band gap [27], and a novel spatially dependent effective scheme to account for the doping on the semiconductor side.
We apply this novel DFT+NEGF approach to study the characteristics of a Ag/Si interface relevant for PV applications [12,23,[28][29][30][31][32][33][34][35][36][37][38].Specifically, we focus on the (100)/(100) interface [28,37] and on the dependence of I-V bias characteristics on the semiconductor dopingnotice that the method is completely general and can be used to describe other M-S interfaces with different crystalline orientations.We consider a range of doping densities for which the interface changes from rectifying to Ohmic.We demonstrate that the "Activation Energy" (AE) method routinely employed to analyse M-S contacts systematically overestimate the value of Φ, with an error that is both bias and doping dependent, due to the assumption of a purely thermionic transport mechanism across the barrier.Conversely, we show how an analysis of the I-V bias characteristics based on the DFT+NEGF electronic structure data provides a coherent picture of the rectifying-to-Ohmic transition as the doping is varied.Finally, we also show that a slab model does not provide a good representation of the interface electronic structure when doping in the semiconductor is taken into account.This is due to the inability of the semiconductor side of the slab to screen the electric field resulting from the formation of the interface.
The paper is organized as follows.Section II and Section III describe the computational methods and the system models employed in this work, respectively.Section IV A presents the calculated I-V bias characteristics and the validation of the AE method based on the calculated data.Section IV B deals with the analysis of the I-V bias curves in terms of the electronic structure of the interface as obtained from the DFT+NEGF calculations.In Section IV C, the simulations are compared to finite-size slab models.The main conclusions are drawn in Section V.
The one-electron KS valence orbitals are expanded using a linear combination of double-ζ PAOs including polarization functions (DZP).The confinement radii r c employed are 4.39 Bohr, 7.16 Bohr, 7.16 Bohr for the Ag 4d, 5s and 5p orbitals, and 5.40 Bohr, 6.83 Bohr, 6.83 Bohr for the Si 3s, 3p, 3d orbitals, respectively.The ionic cores have been described using Troullier-Martins [41] norm-conserving pseudo-potentials [42].The energy cutoff for the real-space grid used to evaluate the Hartree and xc contributions of the KS Hamiltonian has been set to 150 Ry.Monkhorst-Pack [43] grids of k-points have been used to sample the 3D (2D) Brillouin zone in the DFT (DFT+NEGF) simulations.We have used an 11 × 11 × 11 grid of k-points for the bulk calculations, and a k-points grid of 18 × 9 × 1 (18 × 9) for the DFT (DFT+NEGF) simulations of the interface.Geometry optimizations have been performed by setting the convergence threshold for the forces of the moving atoms to 2 × 10 −2 eV/ Å.In all the simulations, periodic boundary conditions (PBCs) were used to describe the periodic structure extending along the directions parallel to the interface plane.In the slab DFT simulation, Dirichelet and Neumann boundary conditions were applied in the direction normal to the interface on the silver and silicon sides of the simulation cell, respectively, whereas in the DFT+NEGF simulations the same direction was described using Dirichelet boundary conditions at the two boundaries between the interface and the bulk-like electrodes.

A. "Spill-in" terms
As described in Ref. 26, the DFT+NEGF method used to simulate the infinite, non-periodic Ag(100)/Si(100) interface relies on a two-probe setup, in which a left (L) and a right (R) semi-infinite electron reservoirs are connected through a central (C) region containing the interface.Once the chemical potentials µ L,R of the reservoirs have been defined, a self-consistent (SCF) KS procedure is used to obtain the electronic density in the C region.The main quantity being evaluated in the SCF cycle is the density matrix required to express the electronic density of the C region in the basis of PAOs centered in the same region, DCC .Assuming µ L > µ R , DCC takes the form where ΣLL is the self-energy matrix describing the coupling of the central region to the semi-infinite L reservoir, and the Green's function of the central region ḠCC is obtained by with SCC and HCC being the overlap and Hamiltonian matrices associated with the PAOs centered at the C region, and ΣRR being the self-energy matrix of the R reservoir.However, even if the DFT+NEGF method provides an elegant scheme to evaluate DCC , it should be noticed that solving Eqs.1-2 is not sufficient to obtain the correct Hamiltonian and the electronic density of the C region.The reason is that the relevant integrals involved in Eqs.1-2 are evaluated only in the region of space encompassing the C region, and only for the atoms localized in that region.As a consequence, the tails of the PAOs located close to both sides of the L/C and R/C boundaries, which penetrate into the neighboring regions, are not accounted for (see Fig. 1).To correct this behavior, we introduce additional corrective terms, that we name "spill-in".For the Hamiltonian, corrective terms are applied to both the two-center and three-center integrals.Specifically, if two PAOs φ i and φ j centered in the C region lie close to a boundary, e.g. the L/C one, the corrected Hamiltonian in Eq. 2 will include also the matrix element H ′ i,j = φ i |V LL ef f (r)|φ j associated with the tail of the PAOs protruding into the L region, V LL ef f (r) being the periodic KS potential of the semi-infinite L reservoir (Fig. 1a).Similar arguments hold also for the Hamiltonian three-center non-local terms.For the electronic density of the C region, additional contributions are included for each pair of PAOs φ i and φ j located close to a boundary when at least one of them is centered at the neighboring reservoir region.In total, two new contributions 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): which can be further distinguished based on whether both ( DLL ) or just one ( DLC ) of the two PAOs involved is centered at the L region.In the calculations presented 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 neighboring reservoir, i.e. that the "screening approximation" is verified -see Ref. 26 for additional details.We stress that including these "spill-in" terms is essential to ensure 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.

B. Exchange-correlation potential
Further complications in describing the Ag(100)/Si(100) interface arise from the fact that one of its sides is semiconducting.In fact, a major problem affecting the description of metal-semiconductor interfaces is the severe underestimation of the semiconductor band gap in DFT calculations using (semi-)local xc-functionals based on the local density approximation (LDA) or on the generalized gradient approximation (GGA) [45].For model calculations based on few-layer thick fully periodic systems, such an underestimation has been shown to result in unrealistically low Schottky barriers at the interface [15,46].In order to remedy this drawback, we have evaluated the electronic structure of the LDA-optimized interface geometries using the Tran-Blaha meta-GGA xc-functional (TB09) [27].The TB09 xc-functional has been shown to provide band gaps in excellent agreement with the experiments for a wide range of semiconductors including silicon, at a computational cost comparable to that of conventional (semi-)local functionals.In the TB09 xc-functional, the exchange potential υ TB09 x (r) depends explicitly on electron kinetic energy τ (r), with , N being the total number of KS orbitals, ψ i (r) the i-th orbital, ρ(r) the electronic density and υ BR x (r) the Becke-Roussel exchange potential [47].The parameter c in equation ( 4) is evaluated self-consistently and takes the form where Ω is the volume of the simulation cell and the two empirical parameters α = −0.012(dimensionless) and β = 1.023Bohr 1 2 have been fitted to reproduce the experimental 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 order to reproduce the experimentally measured band gap of bulk silicon, E exp gap = 1.17 eV [44].This has been accomplished by calculating the band gap of bulk silicon at fixed values of the c parameter in a range around the self-consistently computed value in which the variation of E TB09 gap with c is linear.Then, the optimal value of c has been determined as the intersect between the value of E exp gap and a linear fit to the computed values of E TB09 gap (Fig. 2a).Using the TB09 xc-functional with the c parameter fixed at the optimal value determined using this procedure (hereafter, TB09-o), we calculate the indirect band gap of bulk silicon to be E TB09−o gap = 1.169 eV (Fig. 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) interface reported in this work.We have checked that the band structure of bulk silver calculated using the TB09o functional is very similar to that calculated using the LDA, which is known to perform well for noble metals.

C. Semiconductor doping
The last requirement to describe realistically the electronic structure of the Si(100)/Ag(100) interface is to account for the doping on the silicon side of the interface.Here, doping is achieved in an effective scheme by introducing localized charges bound to the individual silicon atoms.More specifically, in ATK [39] the total selfconsistent electronic density ρ tot (r) is defined as [40]: where is the sum of the atomic densities of the individual neutral atoms of the system.As each atomic density ρ I (r) is a constant term, it can be augmented with a localized "compensation" charge having 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" charges, an effective doping can be achieved both in the DFT and in the DFT+NEGF simulations.In the former, the "compensation" charge added to each silicon atom is neutralized by explicitly adding a valence charge of the opposite sign, so that the system remains charge neutral.In the latter, the "compensation" charge is neutralized implicitly by the carriers provided by the reservoirs, and the system is maintained charge neutral under the condition that the intrinsic electric field in the system is zero.This effective doping scheme has the advantage of (i) not depending on the precise atomistic details of the doping impurities, and (ii) being completely independent of the size and exact geometry of the system.

III. SYSTEM
In order to obtain a reliable description of the Ag(100)/Si(100) interface, we have followed a stepwise procedure.Initially, we have carried out a preliminary screening of the interface geometries and bonding configurations by using a 2×1 slab model formed by a 6-layers Ag(100) slab interfaced with a 9-layers unreconstructed Si(100) slab.The calculated bulk lattice constants of silicon (a Si = 5.41 Å) and silver (a Ag = 4.15 Å) are in good agreement with those reported in the literature [23].To match the Ag(100) and the Si(100) surface, we have applied an isotropic compressive strain ǫ xx = ǫ yy = -0.0793along the surface lattice vectors v 1,2 of the Ag(100) surface.We have checked that in the compressed Ag(100) structure, the dispersion of the s-band and its position with respect to the d-band are very similar to those calculated using the equilibrium value of a Ag .The Si(100) surface opposite to the interface has been passivated with hydrogen atoms.The geometry of the resulting 15-layers slab has then been optimized using the LDA by keeping the farthest (with respect to the interface plane) 4 layers of the Ag(100) surface frozen, and by allowing the farthest (with respect to the interface plane) 4 layers of the Si(100) slab to move as a rigid body, thereby freezing only the interatomic distances and angles.All the remaining atoms have been allowed to fully relax.Different starting guesses for the interface structure have been tested, corresponding to different configurations of the Si(100) dangling bonds with respect to the high symmetry fcc sites of the Ag(100) surface.The lowest energy configuration among those considered, corresponding to the Si(100) dangling bonds sitting above the "hollow" fcc sites of the Ag(100) surface, has then been used as a representative model of the interface.
Starting from the lowest energy configuration obtained 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 create two-probe setups effectively describing the infinite, non-periodic interface (Fig. 3a).A final geometry optimization has been carried out using a two-probe setup in which the C region has been described by 8 Ag(100) layers and an undoped silicon layer having a total width W CC Si(100) = 47.84Å.The optimized geometry has been used to construct two-probe setups in which the doping of the silicon side has been taken into account using the effective doping method described in Section II C. We have considered doping densities of n-type carriers (n d ) in the experimentally relevant range Si(100) = 447.92Å, respectively.We have checked that reducing the width of the C region does not have any effect on the results, as long as all the space-charge effects due to the presence of the interface take place within the screening region.Furthermore, we notice how using a two-probe setup also allows to simulate the characteristics of the interface when the L and R reservoirs are set at two different chemical potentials µ L = µ R due to an applied bias voltage qV bias = µ R − µ L .As will become clear later, this allows for a direct comparison to experiments and for analyzing the electronic structure of the interface under working conditions.Finally, to understand to which extent the slab model is able to describe accurately the electronic structure of the infinite, non-periodic interface, we have also considered slab models having a similar interface geometry as that used in the two probe setup (Fig. 3b).Both short and long slab models have been constructed, in which the width of the Si(100) layer used to describe the silicon side of the interface has been set to either W = 98.62 Å.Notice how these values of W slab Si(100) are many times larger than those used for similar studies of the Ag(100)/Si(100) interface reported in the literature [37].

A. Device characteristics and validation of the activation energy model
Fig. 4a shows the current-voltage (I-V bias ) characteristics calculated for the two-probe setup at low (n d = 10 18 cm −3 ), intermediate (n d = 10 19 cm −3 ) and high (n d = 10 20 cm −3 ) doping densities of the Si(100) side of the interface.A strong dependence on the doping concentration is evident.At low doping, the interface shows a well-defined Schottky diode-like behavior: the forward bias (V bias > 0 V) current increases about six orders of magnitude in the range of V bias [0.02 V -0.5 V], whereas the reverse bias one (V bias < 0 V) varies only within one order of magnitude in the corresponding range.The diode-like asymmetry in the I-V bias curves persists at intermediate doping, although it is less pronounced than 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 bias curve becomes highly symmetric, 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] where q is the elementary charge, k B is the Boltzmann constant, T is the temperature, I 0 is the saturation current and n is the so-called ideality factor.The latter accounts for the deviation of the I-V bias characteristics from those of an ideal diode, for which n = 1.Fitting the simulated data at forward bias to Eq. 7 allows to extract n from the slope of the fitted curves.In Fig. 4b the fitted curves are compared to the forward bias data.The latter are presented using an alternative form of Eq. 7, which allows for a better comparison with the fitted curves as I/(1 − e −qV bias /kB T ) varies exponentially with V bias in the fitting interval considered, viz.0.02 V ≤ V bias ≤ 0.08 V.At low doping, n = 1.09, indicating that the system behaves essentially as an ideal Schottky diode.At intermediate doping, n = 1.82, and the system deviates significantly from the ideal behavior.At high doping, n = 2.40, consistently with the observation that the system does not behave anymore as a Schottky diode.
The I-V bias simulation allows to test the reliability of the experimental procedures used to extract the Schottky barrier Φ.In particular, we focus on the so-called The solid lines are linear fits to the data.The right-hand side of Eq. 10 has been plotted using the value of Φ AE calculated at V bias = 0.02 V, which approaches the value of Φ at V bias = 0 V. (e,f) Schottky barrier height Φ AE evaluated using Eq. 10 as a function of V bias .
"Activation-Energy" (AE) method, which does not require any a priori assumption on the electrically active interface area A [2].In the AE method the I-T dependence is measured at a small constant V bias .Over a limited range of T around room temperature, assuming that the Richardson constant A * and Φ are constant, the I-T characteristics can be described by the expression Following Eq. 9, the Schottky barrier height Φ AE can be extracted from the ln(I/T 2 ) vs. 1/T data using in which n is the ideality factor extracted above.Fig. 5a,b shows the simulated AE plots (as Arrhenius plots) at different values of V bias for low and intermediate doping densities, at which the interface still displays clear Schottky diode-like characteristics.The I-T dependence has been evaluated in a linear response fashion, using the Landauer-Büttiker expression for the current, ]dE with the transmission coefficient T (E, µ L , µ R ) evaluated selfconsistently at an electron temperature of 300 K. Fully self-consistent simulations performed for selected temperatures show that this approach is valid within the range of T considered, 250 K ≤ T ≤ 400 K.
Ideally, for a given doping the Schottky barrier depends exclusively on the M-S energy level alignment at the interface and therefore, disregarding image-force lowering effects, should remain constant with V bias [2].This implies that in Eq. 10, the left-hand side should equal the right-hand side at any value of V bias .However, in the present case this condition is not verified, as the variation of the right-hand side term with V bias is larger than that of the left-hand side term (see Fig. 5c,d).Indeed, for n d = 10 18 cm −3 (n d = 10 19 cm −3 ), a linear fit to the calculated values of the left-hand side of Eq. 10 gives a slope of -664 meV/V (-177 meV/V), whereas the slope associated to the variation of the right-hand side term is -917 meV/V (-549 meV/V).
Following the procedure in the AE method we use the value of n obtained from Fig. 4 to subtract the bias dependence.The result is shown in Fig. 5e,f and it can be seen that this leads to an unphysical increase of Φ AE with V bias .The error becomes more severe as n d is increased.At low (intermediate) doping, Φ AE varies from by 30% (325%) in the range of V bias considered, leading to a change ∆Φ AE = 3.73 k B T (∆Φ AE = 7.31 k B T ).Thus, the intrinsic accuracy of the AE method depends strongly on multiple factors.On the one side, the non- linear increase in Φ AE with V bias suggests that using a single value of V bias is not sufficient to obtain an accurate estimate of Φ.On the other side, the change in ∆Φ AE with n d at a given V bias indicates that the AE method is unsuited for comparative analyses of the variation of Φ with doping.These facts call for a more direct and general strategy for the characterization of M-S interfaces under working conditions.

B. Electronic properties of the interface
A strong advantage of the DFT+NEGF simulations is that they allow the visualization of the electronic strucof the interface and the direct tracking of its changes when n d and V bias are varied.This makes it possible to analyze the calculated I-V bias characteristics in terms of the electronic structure of the interface.Fig. 6 shows the local density of states [48] (LDOS) of the two-probe model at equilibrium (i.e., at V bias = 0 V) along the direction normal to the interface at the different doping densities considered.Increasing the doping has a two-fold effect on the electronic properties of the system: on the one side, W D decreases from ∼200 nm to ∼20 nm when the doping is increased from n d = 10 18 cm −3 to n d = 10 20 cm −3 , as a direct consequence of the increased screening of the n-doped silicon.Increasing n d also shifts the Fermi level towards the silicon conduction bands.In particular, at n d = 10 18 cm −3 the conduction band minimum (CBM) of silicon at Z > W D lies at E − µ L,R = +20 meV, whereas at n d = 10 19 cm −3 and n d = 10 20 cm −3 it lies at E −µ L,R = -40 meV and E −µ L,R = -100 meV, respectively.It is also worth noticing how the macroscopic average [49] of the Hartree potential along the direction normal to the interface, V H (blue lines in Fig. 6), follows the profile of the silicon CBM close to as well as far away from the interface.Similarly to what happens for the electronic bands, V H becomes constant at Z > W D , indicating that the electronic structure starts to resemble that of the infinite periodic bulk.A closer analysis also reveals that a finite density of states extends considerably on the semiconductor side of the interface, due to penetration of the metallic states into the semiconductor side [50][51][52].
Due to the lack of a well-defined electronic separation between the metal and the semiconductor, it is difficult to provide an unambiguous value for Φ based on the electronic structure data only.However, due to the fact that V H closely traces the CBM, it is still possible to estimate the Schottky barrier by defining Φ pot as the difference between µ L and the maximum of V H on the semiconductor side of the interface, V max H (see Fig. 6).We calculate Φ pot = 412 meV and Φ pot = 342 meV for n d = 10 18 cm −3 and n d = 10 19 cm −3 , respectively.The solid line is a guide to the eyes.Filled squares: variation of the slope-dependent term of Eq. 10 (same as in Fig. 5c).The solid line is a guide to the eyes.Filled triangles: φF as a function of V bias .The dashed line shows the bias dependence V bias /n from Eq. 10.The energy on the vertical axis is relative to the semiconductor chemical potential µR.(b) Same as (a), but for intermediate doping.
For n d = 10 20 cm −3 the barrier is considerably lower, Φ pot = 133 meV, reflecting the more pronounced Ohmic behavior observed in the I-V bias curves.Focusing on the low and intermediate doping cases, it can be noticed how the values of Φ pot are considerably larger than those of Φ AE at V bias → 0 V.In particular, at low doping Φ pot − Φ AE = 112 meV, whereas at intermediate doping the difference is even larger, Φ pot − Φ AE = 286 meV.
A consistent physical picture that rationalizes the I-V bias curves can be obtained by studying the doping dependence of the spectral current 7b,d shows the profiles of V H obtained at forward bias in the bias range 0.02 V < V bias < 0.2 V for low and intermediate doping densities.As V bias is increased, V max H shifts towards higher energies due to image-force effects [2], and V H becomes progressively flatter on the semiconductor side.The overall result of these changes is a decrease of the barrier φ F associated with the thermionic emission process from the Si(100) conduction band to Ag(100) (see Fig. 7a): The associated spectral currents I(E) are shown in Fig. 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 − µ L > Φ pot .However, in the present case I(E) is finite also at E − µ L < Φ pot , indicating that electron tunneling has a non-negligible contribution to I.This contribution is much larger for intermediate than for low doping densities.Indeed, at V bias → 0 V, the position of E(I max ) lies very close to Φ pot in the low doping case, as expected in the case of a nearly ideal Schottky diode.Conversely, at intermediate doping E(I max ) lies well below Φ pot , indicating that electron tunneling has become the dominant transport process.
The trend of I(E) with V bias is consistent with these considerations.At low doping, E(I max ) is pinned to V max H , whereas the onset of finite I(E) at E−µ L < Φ pot moves towards energies, following the variation of V H . On the other hand, at intermediate doping the overall shape of I(E) remains the same as V bias is increased, and the variation of E(I max ) follows closely that of V H .We also notice the presence of a narrow resonance at E − µ L = +0.395eV, whose position is independent on n d and V bias .This is due to a localized electronic state at the interface which is pinned to µ L .
The variation of φ F with V bias can be related to the slope-dependent term of Eq. 10 by assuming Φ = Φ AE in Eqs.10-11, thus allowing for a direct comparison with the AE data (see Fig. 8).Independently of the value of V bias , the slope-dependent term lies always below φ F , due to the missing contribution of electron tunneling in the AE method: the latter assumes that the current has a purely thermionic origin, and consequently predicts a value of φ F lower than the actual one.In agreement with the previous analyses, this deviation is considerably larger in the intermediate doping case, due to the much larger contribution of electron tunneling.

C. Comparison of the two-probe with the slab model
The results obtained using the two-probe model can be used as a reference to validate the use of finite-size models to describe the Ag(100)/Si(100) interface.Such models are integral parts of the band alignment method often used to evaluate Φ using conventional DFT [53][54][55][56].The method relies on aligning the electronic band structures of the two bulk materials forming the interface on an absolute energy scale by using a reference quantity, often V H [49].The perturbation of the bulk electronic structure in each material due to the presence of the interface is accounted for by either a slab [37] or a fully periodic [57] model.V H is then used as a common reference to align the electronic structure obtained from independent calculations of the two bulk materials.Despite its widespread use, this strategy relies on two drastic assumptions.Firstly, it is implicitly assumed that the electronic properties of the interface are independent of the doping level of the semiconductor.Moreover, it is assumed that the electronic properties in the central part of each side of the interface model are a good approximation to those of the two bulk materials.
Fig. 9 shows a comparison between V H obtained at the different doping densities considered for the two slab models (short and long, see Section III) and for the twoprobe setup.We notice that introducing an effective dop- ing in the slab model, which was not taken into account in previous slab models for Ag/Si interfaces [37], attempts to better mimic the two-probe simulation in which silver is interfaced with n-doped silicon.The profiles of V H of the three different systems have been aligned according to the value of µ L , the side at which Dirichlet boundary conditions are used for the three systems.Irrespectively of the doping level, the doped short slab model provides a poor description of the variation of V H at the interface.In particular, V max H is always ∼200 meV higher that that obtained for the two-probe model.Furthermore, on the Si(100) side of the interface, V H does not decay correctly with the distance from the interface for the short slab model and, even more importantly, it does not converge to a constant value.The situation improves by increasing the width of the Si(100) layer.For the long slab model, at increasingly larger doping densities the profile of V H resembles more and more that of the twoprobe model.Indeed, in the best case scenario, i.e., at n d =10 20 cm −3 , the profile of V H evaluated using the long slab becomes constant in the center of the Si(100) region, albeit still higher than that of the reference by ∼100 meV.The limitations of the slab model in reproducing the electronic structure at the interface are also evident by looking at the corresponding LDOS plots (see Fig. 10).Similarly to what is observed for V H , the short slab model fails to reproduce the band bending observed at low doping, as well as the correct trend in the decrease of W D as doping is increased.The latter is qualitatively reproduced using the long slab model.However, these modest improvements going from the short to the long slab model come at the expenses of a much higher computational cost.In fact, each DFT calculation for the long slab model takes on average 338.6 s/step.Conversely, each DFT+NEGF calculation using the two-probe model is approximatively one order of magnitude faster, taking on average 46.6 s/step.This suggests that, in addition to computational efficiency, there are also more fundamental reasons for making DFT+NEGF the method of choice for describing M-S interfaces, as in the twoprobe setup the two main assumptions of the band alignment method are naturally lifted.We emphasize that, although the results presented in this paper are specific to the Ag(100)/Si(100) interface only, similar conclusions are likely to hold true for all systems for which the poor screening on the semiconductor side of the interface results in space-charge effects that extend over widths of the order of several nanometers.

V. CONCLUSIONS
In this work, we have presented an approach based on density functional theory (DFT) and nonequilibrium Green's functions (NEGF) for realistic metalsemiconductor (M-S) interfaces modeling.Our approach is designed to deal effectively and correctly with the nonperiodic nature of the interface, with the semiconductor band gap and with the doping on the semiconductor side of the contact, and allows for a direct theory-experiment comparison as it can simulate I-V bias characteristics.Using a Ag/Si interface relevant for photovoltaic applications as a model system, we have shown that our approach is a better alternative to (i) analytical approaches such as the "Activation Energy" (AE) method to analyze the I-V bias characteristics of non-ideal rectifying system with non-negligible tunneling contribution, and (ii) finite-size slab models to describe the interface between metals and doped semiconductors.This DFT+NEGF approach could pave the way for a novel understanding of M-S interfaces beyond the limitations imposed by traditional analytical and atomistic methods.
gies for their input on this work.DS acknowledges sup-port from the H.C. Ørsted-COFUND postdoc program at DTU.

FIG. 1 .
FIG. 1. Scheme showing the two-center "spill-in" terms used 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 φi located close to the L/C boundary.In (a,b) the blue shaded regions indicate the integrals performed in the C region, whereas the red and green regions indicate the "spill-in" terms for the Hamiltonian and for the electronic density, respectively.

FIG. 2 .
FIG. 2. (a) Fitting procedure for the TB09 xc-functional c parameter.Squares (blue): calculated indirect band gap of bulk silicon E TB09 gap obtained for different values of the c parameter.Dashed line (blue): fit to the computed data of E TB09 gap vs. c obtained by linear regression.Dotted line (black): experimentally measured bulk silicon band gap [44].Dasheddotted line (orange): optimal value of the c parameter (optc), obtained as the intersect between E TB09 gap and the fit to the E TB09 gap data.(b) Region around the indirect band gap in the bulk silicon band structure calculated using the optimal c parameter determined from (a).

FIG. 3 .
FIG.3.Geometries employed to simulate the Ag(100)/Si(100) interface using two-probe models (a) or slab models (b).Silver, silicon and hydrogen atoms are shown in grey, beige and white, respectively.

1 / 2 d
[10 18 cm −3 -10 20 cm −3 ].As discussed in more detail in Section IV, the width of the Si(100) layer needed to describe accurately the interface in the two-probe simulations depends on the size of the depletion region (W D ) on the silicon side of the interface.The relation between W D and n d is 1/W D ∝ n , so that progressively narrower C regions can be used as the doping level is increased without any loss in accuracy.Therefore, in the following, the results presented for n d = 10 20 cm −3 , n d = 10 19 cm −3 and n d = 10 18 cm −3 refer to calculations performed with C regions of widths W CC Si(100) = 47.84Å, W CC Si(100) = 197.436Å and W CC

FIG. 4 .
FIG. 4. Calculated I-V bias (a) and forward bias I/(1 − e q|V bias |/k B T )-V bias (b) characteristics at n d = 10 18 cm −3 (blue triangles), n d = 10 19 cm −3 (green squares), n d = 10 20 cm −3 (red dots).In (a), the values of I at n d = 10 18 cm −3and n d = 10 19 cm −3 have been multiplied by a factor 10 and 100, respectively.The solid lines in (b) are fit to the data in the range 0.02 ≤ V bias ≤ 0.08 V using Eq. 7. The ideality factor n extracted from the slope of each fitted curve is reported using the same color as the corresponding curve.

3 n
FIG. 5. (a,b) Empty dots: calculated I-T data at different bias voltages for n d = 10 18 cm −3 (a) and n d = 10 19 cm −3 (b).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 bias .The values of the left-hand side have been extracted from the slope of the fitted I-T curves in (a,b).The solid lines are linear fits to the data.The right-hand side of Eq. 10 has been plotted using the value of Φ AE calculated at V bias = 0.02 V, which approaches the value of Φ at V bias = 0 V. (e,f) Schottky barrier height Φ AE evaluated using Eq. 10 as a function of V bias .

3 FIG. 6 .
FIG. 6. Local density of states (LDOS) of the two-probe setup at equilibrium for n d = 10 18 cm −3 (a), n d = 10 19 cm −3 (b) and n d = 10 20 cm −3 (c).The energy on the vertical is relative to the system chemical potential µL,R.Regions of low (high) LDOS are shown in dark (bright) color.The blue line in each panel indicates the macroscopic average of the Hartree potential VH subtracted the electron affinity of bulk Si and µL,R.The yellow vertical line in each panel indicates the associated Φ pot .
FIG. 7. (a) Scheme of the electronic structure of the Ag/Si interface at forward bias voltage.(b) Profile of VH for different V bias at low doping.The energy on the vertical axis is relative to the electron affinity χ of bulk Si and the metal chemical potential µL.The vertical lines indicate φF at V bias = 0.02 V (blue, solid) and V bias = 0.2 V (red, dashed).(c) Solid curves: spectral current density I(E) for different V bias at low doping.The dashed line indicates the value of Φ pot .VH and I(E) curves calculated at increasingly higher V bias are shown in blue→green→yellow→red color scale.(d,e) Same as (b,c), but for intermediate doping.

3 FIG. 8 .
FIG. 8. (a) Filled circles: energy of maximum spectral current E(I max ) in Fig. 7c as a function of V bias at low doping.The solid line is a guide to the eyes.Filled squares: variation of the slope-dependent term of Eq. 10 (same as in Fig.5c).The solid line is a guide to the eyes.Filled triangles: φF as a function of V bias .The dashed line shows the bias dependence V bias /n from Eq. 10.The energy on the vertical axis is relative to the semiconductor chemical potential µR.(b) Same as (a), but for intermediate doping.

FIG. 9 .
FIG. 9. Profile of VH along the direction Z normal to the interface plane, calculated for the two-probe setup (blue solid line) and for the short (green dotted line) and long (red dashed line) slab models, at doping densities n d = 10 18 cm −3 (a), n d = 10 19 cm −3 (b) and n d = 10 20 cm −3 (c).The vertical black solid line indicated the position of the interface.The vertical green (red) line indicates the position Si(100) layer farthest from the interface in the short (long) slab model.