Resistivity size effect due to surface steps on ruthenium thin films computed with a realistic tight-binding model

The resistivity size effect, sometimes also called the classical size effect, was first reported in 1901.¹ The effect is observed as an increase in the resistivity $ρ$ of thin metallic films with decreasing film thickness $d$⁠. The implications for interconnect technology have become significant as feature sizes have decreased to nanometer scales.² It is generally accepted that diffusive surface scattering of electrons is responsible for the classical size effect. Theoretical descriptions often are based on the Fuchs–Sondheimer (FS) model.^3,4 The FS model is obtained from a solution of the Boltzmann transport equation (BTE) with inelastic electron scattering occurring in the bulk and partially diffusive scattering at surfaces implemented through boundary conditions. Another scattering mechanism occurs via grain-boundary scattering in polycrystalline thin films. This is often described using the Mayadas–Shatzkes model.⁵ Yet efforts to understand and mitigate the resistivity size effect have been hampered by an inadequacy of theoretical modeling tools. Specifically, while it is accepted that surface roughness due to various defects, including vacancies, adatoms, and surface steps, is connected to diffusive electron scattering, there has been a lack of theoretical approaches to elucidate various mechanisms at the atomic scale and their relative importance.

Recently, atomic-scale modeling of transport in thin films based on first-principles electronic-structure methods has begun to address mechanisms for diffusive surface scattering. The size-dependent resistivity of thin Cu films with frozen-phonon disorder and missing surface atoms has been explored using density-functional theory (DFT) calculations employing the non-equilibrium Green’s function (NEGF) approach.⁶ There, it was determined, in general agreement with previous theoretical expressions, that surface scattering adds a term to the resistivity proportional to $1/d$⁠, where $d$ is the film thickness. The connection between the surface roughness and the specularity parameter in the FS model was theoretically explored for Cu films.⁷ Using the same NEGF approach, the influence of steps was determined for thin Cu films with $(001)$ orientation.⁸ In those calculations, the film thickness was six monolayers (ML), corresponding to $d≈1.1$ nm. By computing electron transmission through a region with steps of various heights, it was demonstrated that scattering from steps results in an additional resistivity $Δρ$ proportional to $ω/(ξd)$⁠, where $ω$ is the root-mean-squared film roughness due to steps and $ξ$ is the lateral correlation length of the roughness. This general relation was experimentally verified for tungsten thin films.⁸ However, while these studies have shed significant light on the resistivity size effect, there remains a need for realistic but more computationally efficient approaches with the objective of elucidating transport over longer length scales and thicker films more consistent with the experimental range.

To address this need, we have developed an approach based on tight-binding (TB) electronic structure calculations and the kernel polynomial method (KPM)⁹ for enhanced numerical efficiency. To approach the level of first-principles accuracy, the TB model is parameterized by fitting to the results of DFT calculations of ruthenium. Transport properties, namely, the resistivity, are determined for Ru thin films with surface steps. The results are analyzed using the approach implemented in Ref. 8, where the effect of scattering from steps was treated in the Landauer formalism. The results demonstrate an additional resistivity $Δρ$⁠, which is proportional to $1/d$⁠, consistent with previous theoretical and experimental results.⁸ However, the model predictions indicate a significantly smaller effect than what is predicted for Cu films and also seen in W films.⁸ This is possibly due to a higher ballistic conductance per unit area for Ru in contrast to Cu and W.

In the tight-binding (TB) approach, a local basis set is used with hopping matrix elements described within the Slater–Koster formalism.¹⁰ The single-particle eigenstates $|ψλ⟩$ can then be specified as a linear combination of the tight-binding basis states,

where $λ$ represents the eigenstate, $i$ represents an atomic site, and $κ$ represents one of the tight-binding orbitals on the site. The expansion coefficients $cλ,iκ$ can be determined by direct diagonalization of the TB Hamiltonian. The starting point for the TB model was taken from the model for Ru developed in Refs. 11 and 12. However, in the model developed here, we assume an orthogonal basis in contrast to the non-orthogonal basis used in the original model. Consequently, the model needed to be reparameterized using DFT results. As in the original TB model,¹¹ the DFT energies were shifted so that the total energy was equal to the sum of the eigenvalues of the occupied electronic states. This involves a shift in the eigenvalues $ϵλ$ by an energy $V0$⁠. The new eigenvalues are then $ϵλ′=ϵλ+V0$ and the total energy is

where $f(ϵ)$ is the Fermi–Dirac distribution function at zero temperature and $μ′$ is the electronic chemical potential. The shift $V0$ was made separately for each structure simulated with DFT so that the computed cohesive energy of the fit TB model would accurately reproduce the DFT results. Finally, the TB model includes 5s, 5p, and 4d orbital states, for a total of 9 orbitals per Ru site.

Parameterization of the model was accomplished by a detailed fit to electronic band structures for bulk Ru obtained from the DFT calculations. First-principles DFT results were obtained using Quantum Espresso.^13–15 The pseudopotential used was obtained from Pslibrary 1.0.0¹⁶ and was generated via the projector-augmented wave (PAW) method¹⁷ and includes contributions from 4s, 5s, 4p, 5p, and 4d orbitals. The Perdew–Burke–Ernzerhof (PBE) exchange-correlation functional¹⁸ was used in the calculations. Self-consistent calculations were performed using a primitive unit cell with a $20×20×20$ Monkhorst–Pack mesh.¹⁹ The cutoff energies used for the wave function and density were 60 and 480 Ryd, respectively.

While a description of the original TB model has been reported elsewhere,¹¹ some details on the dependence of the on-site energies and hopping integrals are given in the  Appendix. Fitting to the DFT database was accomplished by a random walk in parameter space using the original parameters as a starting point. At each iteration, a small random change was made to each of the TB parameters. This new trial set was used to build the TB Hamiltonian and calculate the band structure for the atomic structure and $k$-points specified in the training data. The training data consisted of 7 hcp structures with different lattice parameters a and c to capture the distance dependence in the hopping integrals and the environmental dependence for the on-site energy terms. For each of the seven structures, there were 40 total $k$-points sampling a path through Brillouin zone both along high-symmetry directions and at high-symmetry points. At each $k$-point, the character corresponding to each eigenvalue was matched between the DFT and TB calculations. The character table corresponds to the representations of the space group $P63$mmc for the hcp lattice. We note that without matching representations, the risk of misfitting the bands is substantial, especially in cases in which the bands are very close together in energy or in the case of band crossings. A cost was assigned to the proposed parameterization by calculating the total sum of squared differences between the DFT and TB energy bands at each $k$-point and for each structure. In computing the cost, each band had an associated weight which controlled how much influence that band had on the final fit. The weights controlled the relative contribution from each band to the total sum. Since this TB model was intended to be used to calculate transport properties, the highest weights were given to the bands in a region centered on the Fermi energy, extending 3 eV above and 3 eV below the Fermi energy, which will be referred to as the conduction region. The low-lying, occupied bands in the region below the conduction region were weighted about 20% less than those in the conduction region. Finally, the higher energy bands in the 10–30 eV range, which stay far away from the conduction region, were weighted the least; up to 90% less than those in the conduction region. In this way, the cost function favors fitting of the electronic states responsible for electron transport. During each iteration of the fit, the cost function was compared to the previous iteration. In cases where the cost increased, the trial parameters were rejected. If the cost decreased, the trial parameters were accepted. The iterative minimization was allowed to proceed until a minimum value of the cost function was obtained. The particular model obtained in this way depends on the specific DFT data points used and the details of the cost function. Moreover, the minimization procedure may not represent a global minimum.

The resulting model showed quite good agreement for the electronic band structure, especially near the Fermi energy $EF$⁠. Sample results corresponding to the TB ground-state structure are shown in Fig. 1 for $a=2.74$ Å and $c=4.357$ Å along a few high-symmetry directions in the Brillouin zone. In addition to excellent fits to the electronic band structure, the model is also able to reproduce the cohesive energy over a range of atomic volumes. In Fig. 2, the cohesive energy extracted from both the TB model and the DFT calculations is shown. The specific points correspond to the cohesive energy for a given volume with $c/a$ optimized. The resulting cohesive energy curves were fitted with the Birch–Murnaghan equation of state^20,21 to obtain the atomic volume $V0$⁠, cohesive energy $Ecoh$⁠, bulk modulus $B$⁠, and pressure derivative $B′$ of the bulk modulus. These values are shown in Table I along with a comparison to experiment^22–24 and the results obtained from the original non-orthogonal TB model.¹¹ The results demonstrate good agreement between TB, DFT, and experiment where results are available.

When simulating film structures, it was found using direct diagonalization that under-coordinated surface sites have excess electron site occupancy in comparison to bulk sites. We define the excess occupancy of a site $Qi$ using the TB expansion coefficients,

where $Zi=8$ for each Ru site corresponding to the number of valence electrons. For bulk Ru sites, $Qi=0$⁠. Results for $Qi$ are shown in Fig. 3 for a surface step in a Ru thin film. While small deviations $Qi≠0$ might be physical, Fig. 3 shows values that correspond to as much as one additional electron per surface site, resulting in a rather large negative surface charge density. Previous TB models of metallic surfaces that have included Coulomb self-consistency tend to show much smaller local site charges.²⁵ For example, when Coulomb interactions along with onsite energy terms are included, deviations from $Qi=0$ for the $(0001)$ surface of hcp Ti were found to be much less than $±0.01$⁠. This reflects the substantial energy cost associated with accumulation of charge on a few sites. To address this issue without including self-consistent Coulomb interactions, we implemented adjustments to surface site energies to enforce $Qi=0$ at all sites. The results obtained with these corrections are shown in Fig. 3. The approach used a local Lagrange multiplier for each surface site which was chosen to exactly enforce $Qi=0$⁠. The Lagrange multipliers were iteratively varied during the calculation until $Qi=0$ to within the desired threshold was obtained at each site. Modifications to the local site energies were retained for subsequent transport calculations.

In determining the site corrections, the structure displayed in Fig. 3 was used with periodic boundary conditions in the plane. Specifically, the structure in Fig. 3 represents one repeat unit used to determine the adjustments to onsite energies. The repeat unit was chosen to be long enough along the $[112¯0]$ so that adjacent steps did not interact. The site corrections obtained for surface sites and sites in the vicinity of a step were then applied to the much large film structures. An identical approach was applied to determine corrections for $(11¯00)$ films.

Transport calculations were performed using KPM.⁹ This was accomplished by expressing the components of the conductivity tensor in terms of a double KPM expansion in the energy domain. The components $σαβ(E)$ of the conductivity tensor are obtained from the Kubo–Greenwood formula,

where the system volume is $Ω$ and $v^α$ is the electron velocity operator defined using the Heisenberg equation-of-motion,

where $x^α$ refer to the Cartesian components $α$ of the site coordinates. To efficiently evaluate the trace, we used the truncated-basis approximation that employs a set of $R$$(9×N)$ dimensional position-space random vectors $|Rr⟩$ to estimate the trace, namely,

The length of the vector $|Rr⟩$⁠, therefore, represents the $N$ sites with nine orbitals on each site. We define vectors

and

which allow for a more compact expression for the conductivity tensor,

To evaluate the vectors recursively for each random vector $|Rr⟩$⁠, the operators are approximated using KPM with $NT$ terms in the expansion,

and

where $Tn(z)$ are Chebyshev polynomials of the first kind and $gn$ are kernel coefficients. In Eqs. (10) and (11), the energy and the Hamiltonian have been rescaled to fit within a unit interval in order to comply with the Chebyshev polynomial’s convergence range: $ε=(E−Em)/ΔE$ and $h^=(H^−Em)/ΔE$⁠, where $Em$ and $ΔE$ are the band middle point and half width, respectively.

Because the Hamiltonian $H^$ in a tight-binding basis is sparse, evaluating the conductivity tensor in this way scales linearly with system size for a fixed number of random vectors $R$⁠. While the number of random vectors $R$ or terms $NT$ required for convergence of the conductivity tensor may increase with the system size, often this dependence is sublinear, resulting in an overall scaling which is quite advantageous with respect to exact diagonalization.

All KPM conductivity calculations were done using the Jackson kernel and $NT=242$ terms in the Chebyshev expansions. Explicitly, the coefficients $gn$ defined by the Jackson kernel take the functional form,

For ballistic or quasi-ballistic systems, in the absence of a relaxation time scale such as elastic or inelastic scattering times, the number of terms in the Chebyshev expansion affects the computed conductivity, with the conductivity increasing approximately linearly with the number of moments $NT$ used in Eqs. (10) and (11). Because the objective was to compare to experimental results, $NT$ was chosen everywhere to reproduce the bulk room-temperature resistivity for transport in the basal plane. Thus, the truncation of the expansion acts similarly to an inelastic relaxation process, mixing states within an energy window of $ΔE/NT$⁠.²⁶ 

Bulk KPM calculations were done on a rhombohedral supercell containing 327 680 atoms, which measured approximately 17.5 nm on each side. Bulk BTE calculations were done using a two-atom primitive cell over a $64×64×40$ $k$-point mesh. This particular $k$-point mesh was chosen to replicate an equivalent sampling of the Brillouin zone for both the BTE and KPM calculations. The Cartesian $z$-axis was chosen to lie along the crystallographic $[0001]$ direction, and the Cartesian $x$-axis was chosen to lie along the $[112¯0]$ direction. This places the basal plane of the crystal in the $xy$-plane. As a result, there are only two unique, non-zero components in the conductivity tensor. For the $σxx$ and $σyy$ components of the conductivity tensor, 36 random vectors were sufficient for convergence and for the $σzz$ component, 108 random vectors were used. Calculations of the energy-dependent conductivity were tested for convergence relative to the random basis by conducting a t’-test at each energy. Results were considered converged when increasing the number of random vectors used to calculate the stochastic trace yielded statistically indistinguishable results for each energy. The conversion from conductivity to resistivity was straightforwardly done by inverting the (diagonal) conductivity tensor.

For KPM calculations of thin films, orthorhombic supercells were generated with fixed length and width and varying thicknesses. Transport was computed always in the plane of the film along the $[112¯0]$ direction. The length along $[112¯0]$ was approximately 160 nm. The width of the films perpendicular to the transport direction was approximately 5.7 nm for $(0001)$ films and 5.2 nm for $(11¯00)$ films. For BTE calculations on smooth films, orthorhombic unit cells were generated for each film thickness. To ensure the sampling of the Brillouin zone was commensurate with the KPM calculations, transport calculations using BTE were done over a $584×12×0$ $k$-point mesh. For all film calculations, the Cartesian $z$-axis was chosen perpendicular to the surface, and the Cartesian $x$-axis was aligned with the $[112¯0]$ direction. Steps on thin films were generated by removing or adding one layer of atoms from the top surface of a given film, with the step perpendicular to the transport direction. For a section of film with a given thickness, step lengths were pulled randomly from a normal distribution centered on the intended correlation length. For this study, steps with a correlation length of 10 nm were generated with a 5 nm standard deviation. For each film thickness, an ensemble of stepped surfaces were generated for transport calculations. The convergence of calculations for films with steps was determined by comparing the standard error for the ensemble calculation to the standard error for the random-vector truncated basis. When the standard error of the ensemble was of the same order of magnitude or smaller than the standard error introduced by the truncated basis, the results were considered converged. For the surface disorder examined in this study, it was found that a relatively small ensemble consisting of five disorder realizations for each thickness was sufficient to meet the conditions for convergence. For all thin-film KPM calculations, 36 random vectors were used to calculate the stochastic trace. Finally, the error bars shown for smooth film KPM calculations reflect the uncertainty introduced by the truncated random vector basis used for the stochastic trace. For the stepped-film KPM calculations, the error bars shown reflect both the combined uncertainty of the ensemble averaging and that of the truncated random vector basis.

Short representative segments of stepped films are shown in Figs. 4 and 5 for the $(0001)$ and $(11¯00)$ surfaces, respectively. The figures include labels for the thickness $d$ and width $w$⁠. The distance between the steps is also indicated. The distance between the steps is comparable to the average step correlation length $ξ=10$ nm in both cases. Finally, periodic boundary conditions were applied in the plane, with the width $w$ the periodic repeat distance shown for the direction perpendicular to the $[112¯0]$ transport direction.

The results for bulk resistivity determined using KPM are presented in Fig. 6. Past experiments^27–29 report the room temperature, basal plane resistivity for single crystal ruthenium as $7.6μΩ$ cm, and the c-axis resistivity to be $5.8μΩ$ cm. As mentioned previously, the number of moments $NT$ in the KPM calculation was chosen to exactly reproduce the basal plane resistivity. The c-axis resistivity of $6.52μΩ$ cm obtained from KPM calculations is larger than experiment yet still reasonable. For comparison, bulk transport was also determined using the Boltzmann-Transport Equation (BTE) within the single-relaxation time approximation. The elements of the BTE conductivity tensor are computed using the expression,

where $∂f/∂ϵ$ is the derivative of the Fermi function. The relaxation time $τ=9.45$ fs was chosen in the BTE calculations to also reproduce the experimental room-temperature basal-plane resistivity $ρRT=7.6μΩ$ cm. The BTE calculations were done with an electronic temperature $Te=2200$ K to provide sufficient smearing for the electron occupations. Figure 6 also shows the results of BTE calculations within the single relaxation-time approximation, which exhibits excellent agreement with the KPM calculation although resistivity along the $c$-axis is slightly overestimated.

Thin film KPM resistivity calculations are shown in Fig. 7. Somewhat surprisingly, even films without steps show increased resistivity with decreasing thickness, despite a lack of diffusive surface scattering. This behavior is especially pronounced for $(11¯00)$ films. Resistivity averaged over an ensemble of stepped surfaces is also shown in Fig. 7, demonstrating an increased resistivity due to scattering from steps.

The FS model is given by the expression,

in which $ρ0$ is the bulk resistivity, $d$ is the film thickness, $λ$ is the inelastic mean-free path, $κ=d/λ$⁠, and $p$ is the specularity parameter. The data for stepped surfaces in Fig. 7 were fit to the simplified Fuchs–Sondheimer model. Specifically, the data in Fig. 7 for stepped films are fit using the expression,

When surface scattering is entirely specular, $p=1$⁠. By contrast, the diffusive scattering limit corresponds to $p=0$⁠. The simplified FS model fit shown in Fig. 7 was the result of a two parameter fit, with $ρ0$ as one fitting parameter and the value of the product $λ(1−p)$ as the second fitting parameter. For the $(0001)$ films, $ρ0=7.62$$μΩ$ cm and $λ(1−p)=0.208$ nm were obtained from the fit. The much larger resistivity size effect computed for the $(11¯00)$ films resulted in fit parameters $λ(1−p)=1.44$ nm and $ρ0=7.44$$μΩ$ cm. Interestingly, however, much of the thickness dependence is seen even in the absence of steps and hence without any diffusive scattering. This is most notable for the $(11¯00)$ films. Increased resistivity is observed for films with both smooth and stepped surfaces, with the stepped surfaces seeing a small additional enhancement in resistivity. Of the studied surface configurations, the resistivity–thickness relationship is most pronounced in $(11¯00)$ oriented surfaces. The coordination of atoms at the film surfaces is the primary difference between the $(0001)$ and $(11¯00)$ oriented smooth films. The $(11¯00)$ surface is populated by sites with two different coordination configurations, with some having eight nearest neighbors and the others having ten nearest neighbors. In contrast, the $(0001)$ oriented surface is populated by sites with only one coordination configuration, with each site having nine nearest neighbors. The higher-coordinated $(0001)$ surface shows a weaker 1/d trend.

To verify that the resistivity size effect found for perfect films without steps is due simply to differences in the Fermi surface, smooth films were also studied using the BTE. Since the BTE is only dependent upon the nature of the conducting energy bands, any observed increases in resistivity as film thickness decreases can be attributed to band structure effects. BTE calculations of the resistivity are presented in Fig. 8. The calculated resistivity for the smooth $(11¯00)$ oriented surface reveals a fairly strong dependence upon film thickness, whereas the results for the $(0001)$ oriented film show a very weak dependence upon film thickness. This follows the same general trend seen in Fig. 7 for KPM calculated thin film resistivity. This establishes the validity of the KPM calculations and demonstrates that the resistivity size effect can occur in the absence of diffusive electron scattering simply due to changes in the electronic structure due to surfaces. It can be seen, however, that the BTE calculations produce somewhat higher resistivity values in the thin films and appear to extrapolate to a higher resistivity in the limit of a bulk film $d→∞$⁠. This is most likely due to the finite width of the computed films, and differences in how “smearing” occur across the Fermi surface in the KPM and BTE calculations. In the BTE calculations, smearing is accomplished by a Fermi function [see Eq. (13)] and a finite electronic temperature $Te=2200$ K. However, despite the offset in resistivity, the dependence on thickness $d$ is quite similar between KPM and BTE calculations.

The increase in resistivity from monolayer steps cannot account for the experimentally observed resistivity size effect. This is illustrated in Fig. 9, where the KPM calculated resistivity for smooth films is compared to the FS fits to experimental data obtained for single-crystal Ru $(0001)$ films reported in Ref. 28. In the $(0001)$ Ru films used in the experiment, atomic-force micrographs were consistent with monolayer-high steps. However, surface mounds were observed with widths in the range 100–300 nm, which might indicate that the step density of the structures used in KPM calculations was actually greater than that of the films studied in the experiment.²⁸ The parameters used in the FS fits to experiment were $ρ0=7.6μΩ$ cm, $λ=6.71$ nm, and $p=0$ corresponding to completely diffuse scattering.²⁸ The comparison of the KPM results to experimental results for $(0001)$ films suggests that monolayer-high steps do not result in enough diffuse scattering to account for the classical size effect observed in the experiment. In fact, it was noted by the authors of the experimental study, using theoretical insight obtained in a previous study,⁸ that steps would likely contribute less than 1$%$ of the observed resistivity size effect.²⁸ The results computed using KPM agree with this assessment. This indicates the likelihood of other diffuse scattering mechanisms, or potentially a combination of different mechanisms. We will return to this point in Sec. V.

Despite the small effect, the increases in the resistivity due to monoatomic steps can be quantified and interpreted based on the analysis approach recently applied to describe experimental thin film resistivity.⁸ This approach is based on the Landauer formalism applied to electron scattering from steps. The details of this approach are developed in Sec. IV. By comparing the resistivity of smooth films of a given thickness to a stepped film of the same thickness, the rather small classical size effect seen in the KPM calculations can be understood.

We analyze the transport results following the approach described previously in Ref. 8. The basic idea is to treat scattering from step edges within the Landauer formalism such that each step has a transmission probability $η(s)$⁠, which depends on the step height $s$⁠. The resistance $ρ$ of a film with cross-sectional area $A=wd$⁠, where $w$ is the width, is comprised of a “flat film” resistivity $ρff$ and an additional resistivity term due to scattering from $N$ steps which occur along a length $L$ of material,

The term $ρff$ is due to scattering from point defects and phonons and represents the resistivity of a “smooth” film (i.e., without steps). The specific ballistic conductance $g0$ is defined as the ballistic conductance per cross-sectional area $A=wd$⁠,

where $M$ is the number of conductance channels that scales linearly with the cross-sectional area.

Using NEGF methods,⁸ it was found that the transmission coefficient for a step of height $s$ (including both “up” and “down” steps) in a film of thickness $d$ is given by

For a film described by steps all of the same height $s$⁠, the resistivity is then given by

where the factor of $2$ in the denominator accounts for the fact that the total number of steps $N$ is twice the number of “up” and “down” step pairs, and $ξ=L/N$ is the average step-correlation length. In the case where $d≫s$⁠, the resistivity can be approximated by

The resistivity difference $Δρ=ρ−ρff$ can be predicted using the equations above and the computational details for the simulated Ru thin films. For both film orientations simulated, transport was along $[112¯0]$⁠. Consequently, the relevant number of channels for both films is given by

The specific ballistic conductance for this transport direction is

Using the conductance quantum and the lattice parameters for Ru results in $g0=2.998×1015Ω−1m−2$⁠. We use $ξ=10$ nm as the step correlation length. However, the step height for the two films is slightly different. First, we consider $(0001)$ films. In this case, the step height is $s=c/2$⁠. In this case, we obtain for $(0001)$ films,

Second, for the $(11¯00)$ films, the step height is $s=a/3$ and hence, in this case,

In short, the theory predicts a resistivity size effect, which is very similar for the two films.

In Fig. 10, the results of this analysis are shown, with $Δρ$ plotted for both the $(0001)$ and $(11¯00)$ films. Additionally, the theoretically expected $Δρ$ as determined in Eqs. (22) and (23) via the Landauer formalism have been overlaid on the data. It can be seen from this plot that the KPM calculated $Δρ$ follows the same trend and is consistent with the theoretically predicted $Δρ$ for the examined geometries.

This paper presents a practical and realistic approach for the computation of the transport properties of perfect-crystal and defective metallic thin films. Using a realistic TB model parameterized by fitting to DFT electron energy bands, coupled with the efficient KPM numerical approach, systems with $106$ or more sites are easily within reach. The model for Ru developed here was used to describe the resistivity size effect due to monolayer-high steps on Ru thin films with $(0001)$ and $(11¯00)$ orientations.

It was found that resistivity increases with decreasing thickness $d$ even in the absence of steps. This is especially evident for the high-energy $(11¯00)$ surfaces. The $(11¯00)$ includes surface sites with eight and ten nearest neighbors. By contrast, smooth, low-energy $(0001)$ surfaces, which are comprised of surface sites coordinated with nine nearest neighbors, display a weaker resistivity size effect. This potentially suggests a connection between size effects and surface coordination. Specifically, films with higher surface energy and lower-coordinated surface sites may tend to exhibit a measurable size effect even in the absence of surface roughness. Calculations of the BTE show that these differences in orientation and surface coordination determine the transport properties of films without steps. Thus, the results show, for the first time, a resistivity size effect that does not depend on the presence of surface defects.

Steps are found to further increase resistivity by an amount $Δρ$⁠, which also tends to increase with decreasing thickness as $Δρ∝1/d$⁠. This is generally consistent with observations of the resistivity size effect. However, in comparison with experimental results, the effect is quite small. This indicates that other scattering mechanisms are required to explain the observed resistivity size effect. Some possibilities include surface vacancies, adatoms, impurities, size-dependent electron–phonon scattering, and scattering from the substrate. Previous studies of vacancies/adatoms have shown strong scattering and resistivity increases for a surface with half of the atoms removed.⁶ However, thermodynamically it would be expected that vacancy concentrations are quite low at room temperature, and it is more likely that adatoms and vacancies diffuse to form step structures. It is, therefore, still uncertain what scattering mechanisms are most likely responsible for the substantial resistivity size effect seen in metallic thin films.

The small resistivity increases computed for the Ru films appear generally consistent with the theoretical analysis based on previous work.⁸ Essentially, the very small increase $Δρ$ with decreasing volume can be traced to the relatively large specific ballistic conductance $g0=2.998×1015Ω−1m−2$ predicted for Ru. By contrast, the values of $g0=9.5×1014Ω−1m−2$ for $[100]$ transport and $g0=8.7×1014Ω−1m−2$ for $[110]$ transport were reported for first-principles calculations of W films.³⁰ Smaller values of $g0$ have also been reported for Cu. Using NEGF and DFT, $g0=0.996×1015Ω−1m−2$ was obtained for Cu.⁸ Similarly, DFT calculations using the Fermi velocity and the density of states of bulk Cu resulted in $g0=1.10×1015Ω−1m−2$⁠.³¹ Consequently, Cu and W films would be expected to exhibit a more substantial resistivity size effect. Indeed, experimental results do apparently demonstrate much larger increases in W films.⁸ However, it should be noted that the increases seen in experiments require a much smaller value $g0=1.1−1.4×1014Ω−1m−2$ than the value obtained from DFT calculations. Moreover, experimental results indicate a temperature dependence which is not explained by Eq. (19), although this effect was possibly due to small errors in the low-temperature value of $ρff$⁠.⁸ 

In summary, computational results for Ru films show that surface steps produce a very small resistivity size effect, and experimental results must be explained by some other mechanism. The behavior found from the results here is generally consistent with previous theoretical studies.⁸ One important result is the validation that larger values of the ballistic specific conductivity $g0$ tend to suppress the resistivity size effect. This observation points to $g0$ as a critical material parameter for thin films designed to mitigate the resistivity size effect. In this respect, Ru appears to be a very promising material, although further work is needed to understand the significantly larger size effect seen in the experiment.

Funding support by the National Science Foundation (NSF) under Grant No. ECCS-1740228 and the E2CDA-NRI Program of the Semiconductor Research Corporation under Task No. 2764.003, as well as support from the Air Force Office of Scientific Research (AFOSR) under Nos. FA9550-19-1-0156 and FA9550-18-1-0063 is gratefully acknowledged. The authors thank K. Coffey for illuminating discussions.

The data that support the findings of this study are available from the corresponding author upon reasonable request.

The authors have no conflicts to disclose.

The tight-binding model described here is based on the general approach outlined elsewhere.^11,12 For completeness, we summarize some of the relevant equations along with the fit TB parameters for Ru. The starting point for the parameterization was taken from the original model reported for Ru.^11,12 However, in the model developed here, we assume an orthogonal basis in contrast to the non-orthogonal basis used in the original model. Consequently, refitting to DFT results was required. The TB model included 5s, 5p, and 4d orbital states, for a total of nine orbitals per Ru site.

The model defines a density function for each atomic site, which depends upon the neighboring sites $j$ within a cutoff radius $Rc$⁠. For the $ith$ atomic site, the density $ρi$ is

in which $Fc(R)$ is a smooth cutoff function given by

with the values $l=0.5a0$ and $R0=14a0$ (⁠$a0$ is equal to the Bohr radius). The values for $R0$ and $l$ were kept fixed and not treated as parameters during the fitting. The on-site energies for the $κ=$ s, p, d orbitals are functions of the density at that site,

where $aκ,bκ,cκ$⁠, and $dκ$ are fitting parameters. The two center Slater–Koster hopping integrals are assumed to take the form of the product between a polynomial and an exponential function,

where $eγ$⁠, $fγ$⁠, and $gγ$ are fitting parameters and $γ$ indicates the type of orbital interaction. The ten interactions considered for $γ$ are $ssσ$⁠, $ppσ$⁠, $spσ$⁠, $ddσ$⁠, $sdσ$⁠, $pdσ$⁠, $ppπ$⁠, $ddπ$⁠, $pdπ$⁠, and $ddδ$⁠. For our TB calculations, a cutoff radius of $Rc=8.731$ Å was used, which corresponds to $16.5a0$⁠. With this cutoff radius, the on-site density and hopping integrals for each atomic site included contributions from over 300 nearest neighbors.

In Tables II and III, we list the model parameters along with the original parameterization for Ru obtained in Ref. 11. The parameters corresponding to the choice of a non-orthogonal basis in Ref. 11 are not shown in Tables II and III.