Fully resolved currents from quantum transport calculations

We extract local current distributions from interatomic currents calculated using a fully relativistic quantum mechanical scattering formalism by interpolation onto a three-dimensional grid. The method is illustrated with calculations for Pt$|$Ir and Pt$|$Au multilayers as well as for thin films of Pt and Au that include temperature-dependent lattice disorder. The current flow is studied in the"classical"and"Knudsen"limits determined by the sample thickness relative to the mean free path $\lambda$, introducing current streamlines to visualize the results. For periodic multilayers, our results in the classical limit reveal that transport inside a metal can be described using a single value of resistivity $\rho$ combined with a linear variation of $\rho$ at the interface while the Knudsen limit indicates a strong spatial dependence of $\rho$ inside a metal and an anomalous dip of the current density at the interface which is accentuated in a region where transient shunting persists.


I. INTRODUCTION
The standard way to measure a bulk resistivity ρ is the four-point-probe technique [1-3] which assumes isotropic current propagation. In the ongoing pursuit of miniaturization in electronics, such measurements have been extended to study the enhancement of resistivity in thin films [4][5][6][7][8][9][10][11] (and wires [10,12,13]) with a common theme being the estimation of a single effective value of ρ for a given thickness (radius) d [13]. When the Fermi wavelength of the conduction electrons is comparable to d, the wave nature of electrons gives rise to finite size effects [14] that must be taken into consideration. If the mean-free-path λ is comparable to d, the whole concept of a local resistivity becomes moot and as illustrated in Fig. 1(a), specular and diffusive reflection from surfaces play a role in determining the current distribution in such films [4,15]. The corresponding case of a multilayer is illustrated in Fig. 1(b) where interfaces give rise to specular and diffusive scattering and the finite transmission through interfaces leads to shunting of current which is typically addressed in experiments using parallel resistivity models [16,17].
The field of spintronics originated with multilayers comprising alternating thin films of magnetic and nonmagnetic metals [18,19] and these continue to play a pivotal role [20][21][22][23]. The accurate estimation of various spin transport parameters is intimately connected with knowing how much charge current flows in the different layers that are of the order of 1-10 nm thick. A very recent attempt to determine this current distribution combined different thin film resistivity models with fourpoint-probe measurements for a large number of samples where the individual layer thicknesses were varied systematically [24]. This indirect approach was made necessary by the absence of a direct method to observe how current flows in the different layers of multilayer samples. Stejskal et al. concluded their study by emphasizing the need for more detailed structural characterization in order to be able to reduce the non-negligible variations they found in the model-dependent allocation of currents to individual layers.
Although the transport properties of metals are known to be dominated by states close to the Fermi energy [25-FIG. 1. Schematics of (a) electron scattering in a thin film geometry of finite thickness d illustrating the role of surface scattering. An electron wave-packet (depicted by the black arrows) undergoing only specular scattering at the surface (blue arrow) results in a uniform current distribution j flowing in the z direction whereas diffuse scattering from a rough surface (red arrow) results in a non-uniform current profile j(x). (b) Scattering at an interface (thick black horizontal line) between two slabs of finite thicknesses d1 and d2 where the horizontal dashed lines represent a surface or the next interface in a bilayer/multilayer geometry. At an interface, a part of the incident electron wave (I) is transmitted (T) and the rest is reflected (R) specularly or diffusely eventually leading to current equilibration along y. A parallel resistivity model describes x-independent resistivities and, corresponding to these, uniform current distributions. A realistic current distribution j1(x) and j2(x) resulting from a resistivity gradient across the interface and from details of the scattering (specular and diffuse) is sketched on the right. Note that the spatial dependence is only expected to be significant over a length scale determined by the mean free path of the materials. Coordinate axes are shown for reference. 27], there have been few attempts to include the full complexities of the Fermi surfaces associated with partially filled d bands in theoretical studies of transport parallel to the surfaces of thin films [28], or in the context of multilayers, parallel to the interfaces, the so-called currentin-plane (CIP) configuration [29][30][31]. The most sophisticated model used by Stejskal et al. [24] was the phenomenological electron gas model of Fuchs [4] and Sondheimer [15], generalized to treat two different surfaces [32] and to include transmission through interfaces between two different metals [16] as well as the effect of grain boundaries [6,33]. The purpose of this paper is therefore to explore the possibility of calculating the spatial distribution of currents in realistic multilayers and thin films entirely from first principles including temperatureinduced lattice disorder [34]. To do so, we introduce a discrete scheme to interpolate local currents [35] calculated using a fully relativistic DFT based scattering code [36] and apply it to evaluate the full spatial profile of currents in thin films of Pt and Au which are of interest to the spintronics community as well as in Pt|Au and Pt|Ir multilayers.
The paper is arranged as follows. In Sec. II A some aspects of the scattering problem that are relevant for the calculation of interatomic currents are briefly summarized. The planar averaging introduced in Ref. 35 is relaxed in Sec. II B to obtain fully spatially resolved local currents on a three dimensional grid. Inspired by fluid physics we introduce streamlines to visualize the current flow in Sec. II C. Although the same methodology can be straightforwardly applied to study the spatial distribution of spin currents, the present work will for simplicity focus on charge currents. In Sec. III the methodology presented in the previous section is illustrated: in the Knudsen limit for a thin film of Au and a Pt|Au multilayer in Sec. III A; in the classical diffusive limit for a thin Pt film and a Pt|Ir multilayer in Sec. III B. The nonnegligible effect of the choice of lead material is examined in Sec. III C and we conclude with a brief discussion and outlook in Sec. IV.

A. Quantum transport
A typical two-terminal transport configuration is sketched in Fig. 2 with a scattering region (S) sandwiched between ideal left (L) and right (R) crystalline leads. In the adiabatic approximation, atoms in the scattering region are displaced from their mean positions with a Gaussian distribution of displacements characterized by a root-mean square displacement ∆(T ) chosen to reproduce the experimental resistivity [37] at a given temperature T [34]. Such disorder would break the translational symmetry completely and make it impossible to solve the Schrödinger equation. To remedy this, we introduce periodic boundary conditions in the xy plane with an N × M "lateral supercell" comprising N and M unit cells in the x and y directions, respectively, whereby the disorder is assumed to be periodic. It turns out that remarkably small supercells are sufficient to eliminate observable effects of the residual periodicity as long as the temperature is not too low [36].
The transport problem now reduces to one of solving the single particle Schrödinger equation inside region S using the propagating Bloch states of the periodic semi-infinite leads as boundary conditions. To do so in practice, we make use of a "wave function matching" (WFM) scheme [38][39][40] formulated [36,41] for a basis of tight binding (TB) muffin tin orbitals (TB-MTO) [42][43][44] and the atomic spheres approximation (ASA) [45]. TB-MTOs form a localized orbital basis |i with i = Rlms where R is an atom site index and lms have their conventional meaning. In terms of the basis |i , the wavefunction Ψ can be expressed as (1) and the Schrödinger equation becomes a matrix equation with matrix elements i|H|j . Ψ is a vector of coefficients with elements ψ i ≡ i|Ψ extending over all sites R and over the orbitals on those sites, for convenience collectively labelled as i R . |Ψ R is a projection of the total wave function |Ψ onto the orbitals on atom R The minimal TB-MTO basis along with the local density approximation (LDA) of density functional theory (DFT) [46,47] makes the scattering problem tractable for scattering regions comtaining 10 4 -10 5 atoms. A detailed description of the TB-MTO-WFM transport scheme can be found in references [41] and [36].

B. Interpolation of interatomic currents onto a three dimensional grid
We begin with expressions [35] for the charge current j P Q c and spin current j P Q sα between atoms P and Q that are given in terms of the block H P Q of Hamiltonian matrix elements and vectors of expansion coefficients i P |Ψ and i Q |Ψ obtained by solving the scattering problem. A summation over lms is implicit. σ is a vector of Pauli spin matrices σ α and α labels the polarization direction of the spin current. The materials whose transport properties we wish to study are crystalline materials or substitutional alloys at finite temperatures whose constituent atoms are displaced at random from the sites of a Bravais lattice. Determining the spatial distribution of currents in the scattering region thus requires interpolating all j P Q onto regular real space meshes as a function of x, y, z. In Fig. 3 we illustrate the discretization of an arbitrary transport geometry. To generalize the interpolation and averaging of currents for an arbitrary geometry generated by translation vectors T 1 , T 2 and T 3 that are not necessarily orthogonal to each other, we first perform an affine transformation T [48][49][50][51] of the translation vectors with a combination of shearing and scaling transformations into a dual space of orthonormal vectors T 1 , T 2 and T 3 that lie along the Cartesian x, y and z coordinate axes. Mathematically, we can express this as where R is the parent coordinate space and O is the dual space. Since the collinearity of points along a given direction in the parent space is mapped into a corresponding collinearity in the affine transformed space ( §23 of [48]), averaging of local quantities in O say x(≡ T 1 ) can be treated as being equivalent to averaging in the direction of the corresponding translation vector (T 1 ) in R. We now apply T to the disordered geometry and map all atomic coordinates from R to O thus transforming the disorder from the parent space. The disordered geometry in O is then divided into boxes whose dimensions D x , D y and D z are determined from the average distance between consecutive atomic layers in the x, y and z directions, respectively, such that the number of boxes is equal to the number of atoms and each box contains exactly one atom; the latter is guaranteed if the temperature-induced atomic displacements are much less than the interatomic separations. The regular lattice of boxes is constructed in such a way that the centres of gravity of the atomic coordinates and boxes coincide. j P Q ≡ (j P Q c , j P Q sx , j P Q sy , j P Q sz ) is imagined as a current through a wire connecting the positions of atoms P and Q with (arbitrary) cross section A P Q , see Fig. 4. Since microscopic details of the spatial distribution of j P Q are unknowable, we assume a homogeneous flux of current between P and Q. The tensor current density of the wire is such that ← → j P Q V P Q = j P Q ⊗d P Q where V P Q = A P Q d P Q is the volume of the wire P Q, d P Q is the vector pointing from P to Q and d P Q its length. The direct product j P Q ⊗ d P Q is estimated in the parent space R as it depends on the components of j P Q and d P Q . Since the affine transformation T preserves ratios of distances between points in a line ( §24 of [48]), the current contribution to each box can be determined in the dual space O using a linear interpolation scheme.
Unlike the planar averaged scheme [35] where a onedimensional interpolation was enough to evaluate the variation in the z direction only, a three dimensional interpolation is implemented here. A general case is shown in Fig. 4 for a pair of atoms P and Q in O with their centres at (x P , y P , z P ) and (x Q , y Q , z Q ) respectively. Note that P and Q can be inside or outside the box "b". Only boxes that are intersected by the wire P Q receive a contribution from the current j P Q which makes it a classic computational problem of "collision detection" where one seeks the point of intersection of the path of an object and a surface of interest. Mathematically, this requires simultaneously solving the equation of line P Q connecting the pair of atoms {P, Q} with equations describing the six faces of each box for all possible {P, Q}. This Illustration of the discrete current scheme where the grid is zoomed to show how the current j P Q between atoms P and Q is interpolated into the box b. The current contribution from "wire" P Q to the box b comes from the segment U V (coloured red) that lies inside the box. amounts to 7 × C(n, 2) × C(n, 1) = 7 2 (n 3 − n 2 ) equations where n is the number of atoms in the geometry with a computational effort that scales as ∼ O(n 9 ). To efficiently perform interpolations for systems with multiple configurations of 10 4 − 10 5 atoms, we instead take advantage of the orthogonality of the dual space O and employ a line clipping algorithm [52] to determine all boxes which are intersected by wire P Q as well as the points of intersection U and V .
The line P Q is first parametrized as for 0 ≤ c ≤ 1. Each box can be described by six boundaries, two for each direction: Equation (7) is only satisfied by the points on the line P Q that lies inside the box. Then c takes a range of continuous values from which the intercepts U and V can be obtained as Note that U = (x P , y P , z P ) when P lies inside b and V = (x Q , y Q , z Q ) when Q lies inside b, which we denote as P ∈ b, Q ∈ b, respectively in the following. We define a parameter β that indicates how much of the wire lies outside the box on either side.
where d AB is the norm of the vector d AB pointing from point B to point A. (In the planar averaged scheme [35], only the projection d z AB of d AB on z, is used to evaluate β P Q .) Since the current changes between Q and P , we make a linear interpolation Thus, for a box of volume V b , the contribution from ← → j P Q V P Q is given by The average current density tensor in the box b is then and we take this value of ↔ j b to be the average current density at the centroid of b. By interpolating all interatomic currents into all boxes b, we obtain the complete spatial variation of the current density. Multiplying the current density by the cross sectional area of the box b perpendicular to z yields the current per unit voltage applied across the leads, a conductance. Summation of the normalized currents for all boxes lying in a given xy plane should then be equal to the Landauer-Buttiker conductance that is calculated independently of the interatomic currents and interpolation scheme. This provides a check of the whole local current formalism. Finally, centroids of the boxes are affine transformed back into the parent space, We observe spatial oscillations in calculated spin currents because of the interference between reflected and incident electron matter waves. Although these oscillations are real, they are not present in semiclassical descriptions of transport. Because they are found to be attenuated away from the leads in parallel with the corresponding decrease in the unscreened particle accumulation, we follow Ref. 35 and use the latter to reduce these quantum fluctuations to facilitate analysis using semiclassical transport formulations. Since lateral supercell sizes do not exceed more than a few hundred atoms in our calculations, we perform such averaging using only the planar averaged unscreened particle accumulation. Details of this averaging can be found in Ref. 35.

C. Current streamlines
From here on we only consider the charge current and will therefore drop the subscript c. Because of the assumed superlattice periodicity, the average current in the y direction, j y (x, z) = 0. The visualization of j(r) reduces to a problem in two dimensions if we average over y and, to do so, we introduce a current stream function ψ(x, z) (not to be confused with the wavefunction) by analogy with the velocity stream function in fluid physics [53]. In the steady state, charge conservation requires that ∇.j = 0 and for the xz plane, this reduces to Defining ψ(x, z) such that automatically leads to (14) being satisfied. ψ(x, z) = constant is a path whose tangent at any point gives the direction of the current vector j = j x i+j z k at that point. This defines a streamline and the volume flow per unit width between streamlines connecting the left (L) and right (R) leads is This region can be thought of as a conducting strip carrying a constant flux of current analogous to a streamtube for incompressible fluid flow such that crowding of streamlines at a region in the flow-field indicates a local increase in the magnitude of the current [53].

D. Mean Free Path
In the relaxation time approximation (RTA), the conductivity is given in terms of the k dependent velocities In the low temperature limit − ∂f ∂ε → δ(ε − ε F ) and (17) becomes an integral over the Fermi surface S F . Assuming additionally that τ (k) = τ (ε(k)) then where D(ε) is the density of states. Both D(ε F ) and v 2 F can be evaluated from standard bulk LMTO electronic structure calculations [54] and since σ ≡ 1/ρ is known [37,55], τ can be expressed as 19) and the mean free path can be estimated as λ = τ v 2 F 1/2 .

III. RESULTS
Different regimes of electron transport can be identified depending on the ratio of the electron mean free path λ to the critical dimension d of the scattering geometry that is the Knudsen number (Kn). When λ d, we are in the classical limit where the flow of current is well described by Ohm's law. The other extreme is the Knudsen limit where λ d, size effects and interface or surface scattering dominate and transport deviates from Ohm's law [56]. To illustrate the three-dimensional current scheme presented above, we consider thin films and two-component ...A|B|A|B... multilayers where the thickness of the ith layer is d i . In this paper, only fcc metals are considered as sketched schematically in Fig. 5 with a charge current flowing in the [001] direction, parallel to the A|B interfaces. A k-point sampling of ∼ 160 N × 160 M for an N × M supercell is used throughout the paper. From now on we use the terms current and current density interchangeably. Unless stated otherwise, all currents are averaged over y and are calculated for a temperature T = 300 K ("room temperature"); when averaging currents over the scattering region, z ∈ S, a few layers close to the leads where transient effects are observed are omitted.
A. Knudsen limit

Au thin film
We begin by calculating the current in a free-standing [110] oriented thin film of Au. The thin film is modelled as a Au|vacuum "multilayer" by alternating 60 atomic layers of Au with five layers of "empty spheres" (with nuclear charge Z = 0 to simulate vacuum [57,58]) in the x direction so that N = 60 + 5. A periodicity of M = 3 layers in the y direction is imposed and the scattering region is 90 atomic layers thick in the z direction, see Fig. 5. A value of the root-mean square displacement ∆ of the atoms in the scattering region was chosen to reproduce the room temperature bulk resistivity of Au, ρ 300 Au = 2.3 ± 0.07 µΩ cm [37,59]. The thickness of the slab in the x direction is approximately 1 2 λ 300 Au where λ = 34.5 nm was estimated in the RTA as described above. Ballistic Au leads were used to minimize transient effects at the lead|scattering region interface, between crystalline and thermally disordered Au. In Fig. 6(a) we plot the charge currentj z (x) flowing in the z direction averaged over the y and z directions as a function of x (black symbols). This shows a gradual concentration of the current away from the surfaces and towards the middle of the film. A strong z dependence is also apparent from plots ofj z (x, z 0 ) averaged over y and z = z 0 ± 5 atomic layers for different values of z 0 ∼ 3, 5, 7, 9 nm measured from the left lead. The current streamlines are plotted in Fig. 6(b) where they are superimposed on a colour map showing the magnitude of the charge current. The colour map shows a larger current density at the centre of Au in both x and z directions. On closer examination, streamlines are not parallel to the z axis but exhibit curvature. This demonstrates that the current distribution has not reached its asymptotic form in z which is not surprising as the length of the scattering region is only about half of the mean free path λ [60]. Apart from a rapid decay of the current density within a few layers of the surface that is described by a "specularity coefficient" p in the Fuchs-Sondheimer framework [4,15], an approximately linear variation of the conductivity from the surface to the middle of the film is observed at the centre of the scattering region furthest from the leads. An "effective resistivity" clearly conceals substantial variation in the current density for film thicknesses that are commonly used in spintronics.
x (nm) 6. Current distribution in a thin film of Au at 335K [59]. (a) Top curve, black symbols: current densityjz(x) obtained by averaging over y and z. Colour symbols:jz(x, z0) obtained by averaging over y and z = z0 ± 5 layers in the z direction for four different values of z0 that are offset from the black curve in steps of ∆j = 0.002 for clarity. The error bars, that are smaller than the symbol size, indicate the average deviation over 10 random configurations of disorder. (b) Streamlines of the current vector j(x, z) = jxi + jzk in the xz plane. The colour contour in the background corresponds to the magnitude of the current.

Pt|Au multilayer
Modern experiments frequently use heterostructures in which the thicknesses d i of constituent layers are comparable to the corresponding bulk mean free paths λ i . To demonstrate the deviation from bulk behaviour we choose a Pt|Au multilayer. Both Pt and Au are fcc metals with only a 2% lattice mismatch so that such an epitaxial multilayer might be prepared without undue structural disorder. Thermal disorder corresponding to the room temperature bulk resistivities of Au, ρ 300 Au = 2.6±0.07 µΩ cm and Pt, ρ 300 Pt = 10.8±0.5 µΩ cm, was used with different mean square displacements in the Au and Pt layers. The mean free path in Au, λ 300 Au = 34.5 nm is almost ten times that in Pt, λ 300 Pt = 3.74 nm [62]. Choosing d Au < λ 300 Au and d Pt > λ 300 Pt should make any size effect apparent at room temperature. We construct a scattering geometry as shown schematically in atomic layers each of Pt and Au in the x direction with a periodicity of 3 layers in the y direction and 90 layers thick in the z direction corresponding to a total of 32400 atoms in the scattering region. A charge current is injected from ballistic Au leads in the z direction. The resulting charge current distribution in the Pt|Au multilayer is plotted in Fig. 7where the error bars, that are smaller than the symbol sizes, correspond to the uncertainty with which the experimental resistivities are reproduced in our scattering calculations.
The average shuntingj Au z /j Pt z ∼ 2.3 i.e, the ratio of the mean current value in Au to that in Pt, indicated by the solid black horizontal lines in each layer plotted in Fig. 7(a), is much lower than expected from the ratio of the bulk resistivities ρ 300 Pt /ρ 300 Au ∼ 4.2. The charge current is seen to be constant and saturated inside Pt while in Au a rapid increase in the two atomic layers adjacent to the interface followed by a continuous variation to the centre of the layer is observed. We note a small, anomalous dip in the current density in the Pt layers next to the interface. The current density variation at the interface and in the Au slab is clearly not amenable to description using a single resistivity. An almost immediate saturation of the total current carried by the Au and Pt layers is observed when these are plotted as a function of z in Fig. 7(b) making our results independent of the length of the scattering region. The streamlines plotted in Fig. 7 are parallel to the z axes inside Pt indicating that there is no net flow of current across the interface into Au as asymptotic shunting of the total charge current density is reached in the z direction. However, a redistribution of the current inside Au is visible from the color contour and small curvature of the streamlines that is similar to what we saw for the Au thin film. As the total current in each layer is independent of z, one might use a parallel resistance model to estimate the current shunting noting, however, that this would not correctly describe microscopic details of the variation near the interface and in the Au layer just presented.

B. Classical limit
The classical diffusive limit is achieved when the mean free path is much shorter than other critical dimensions of the structures being studied, λ d. Memory requirements currently limit the lateral supercell size to H = N × M ≈ 300 − 400 atoms for which the maximum length of scattering region is about L ≈ 130 atoms when spin-orbit coupling is included. Because the computational effort for metallic systems scales as approximately H 2 L [36,41], considerably longer scattering regions can be studied with smaller lateral supercells. To attain the diffusive limit we need to consider high resistivity materials or study elevated temperatures or both.

Pt thin film
We first study charge transport through a thin film of Pt, modelling a free-standing [110] oriented Pt layer as a Pt|vacuum multilayer with N layers of Pt and 5 layers of "empty spheres" repeated periodically in the x direction. The effect of increasing the default periodicity of three atomic layers in the y direction to five atomic layers is studied. To minimize lead|scattering-region interface effects, Pt leads are used. We calculate the average resistivity for different values of N and identify the thickness at which it saturates to the bulk value. At T = 300 K, the bulk resistivity of Pt is ρ 300 Pt = 10.8 ± 0.5 µΩ cm and λ 300 Pt ∼ 3.74 nm.
The thickness dependence of the resistivity of a thin film (or wire) is often studied using the "FS" model formulated some 70-80 years ago by Fuchs [4] and Sondheimer [15] and named after them, in which surface scattering is treated phenomenologically using Boltzmann direction for transport in the z(001) direction and λ is the mean free path of bulk Pt at RT. The separation of periodically repeated thin films by vacuum ("vacuum thickness") is modelled using 5 layers of "empty" spheres in the x direction. Calculations for each thickness are done for two cases with 5 (blue) and 3 (red) atomic layers repeated periodically in the y direction. The dotted lines are calculated using the Fuchs-Sondheimer model [4,15] for three different values of the specularity coefficient p that describes the amount of completely diffusive surface scattering: completely specular (p = 1), partially specular (p = 0.5) and completely diffusive (p = 0). According to (20), choosing p = 1 yields the bulk resistivity ρ b irrespective of d/λ. transport theory. In the case of a monocrystalline "freestanding slab" with only bulk and surface scattering of charge carriers, the thickness dependent resistivity, ρ is given by where d is the thickness of the slab and λ and ρ b are, respectively, the mean free path and resistivity of the bulk material. The "specularity coefficient" p is the fraction of electrons scattered elastically from the surface independent of their velocity and takes values ranging from 1 (specular) to 0 (diffusive). When all the electrons are reflected specularly from the surfaces, the resistivity is identical to that of the bulk; finite-size effects are not considered in the phenomenological FS electron gas model. Values of ρ calculated for different thicknesses of Pt at T = 300 K are plotted in Fig. 8 as a function of d/λ where the thin-film enhancement of the resistivity is clear. The resistivity decays to within a few percent of its bulk value when d ∼ 4λ and follows the FS model with p ∼ 0.5 quite well even though the only surface roughness that is included is what results from thermal disorder. The number of layers M in the y direction has been limited to three throughout this paper in order to be able to study films with a reasonably large thickness given by the number of layers N in the x direction. In Fig. 8 this is shown to be a reasonable compromise for all but the thinnest of films where the results of calculations with M = 3 and 5 are compared. The thickness dependent resistivity analysis confirms that d/λ is the appropriate length scale for transport in slabs with finite thickness.
For N = 60 atomic layers corresponding to a film thickness of ∼ 4.4λ 300 Pt and "bulk-like" behaviour, the charge currentj z (x) resulting from averaging over y and z is 10. Distribution of the charge currentjz in a Pt|Ir multilayer. Forjz(x) in (a) jz(x, y, z) is averaged over y and over z.
The horizontal black lines indicate the asymptotic values obtained by averaging separately over x ∈(Pt,Ir) omitting the region of rapid variation close to the Pt|Ir interface in either layer. Forjz(z) in (b) jz(x, y, z) is averaged over x and y with x ∈ Pt and x ∈ Ir separately. (c) Streamlines of the current vector j(x, z) = jxi + jzk obtained from averaging jz(x, y, z) over y and superimposed on a colour contour of j = j 2 x + j 2 z in the xz plane. The scarcely visible error bars that are smaller than the symbol sizes represent the mean deviation over 10 configurations of thermal disorder.
plotted as a function of x in Fig. 9(a) (black symbols). The current density at different z coordinates (coloured symbols) shows that j z is only weakly dependent on z unlike what we saw in Au. The reason is that λ 300 Pt is much shorter than the length of the scattering region. In Fig. 9(b), the current streamlines in the xz plane are superimposed on a colour map of the magnitude j(x, z) ≡ j 2 x + j 2 z of j(x, z) . Within the uncertainties of the calculation the current is constant except very close to the surface where a rapid decay in the current occurs over a length of 2 nm ∼ 1 2 λ 300 Pt . A more detailed study of transport through thin Pt films will examine the effect of film orientation and surface roughness [63].

Pt|Ir multilayer
We study a Pt|Ir multilayer at an elevated temperature of T = 800 K chosen to make it possible to realise "bulk" like behaviour with scattering region sizes that are computationally tractable. Lattice disorder was introduced as described in Sec. II to reproduce the experimental bulk resistivities of Pt, ρ 800 Pt = 28.1 ± 0.4 µΩ cm and of Ir, ρ 800 Ir = 16.1 ± 0.6 µΩ cm at 800 K [37] for which λ 800 Pt ∼ 1.4 nm and λ 800 Ir ∼ 2.4 nm. Ir and Pt are both fcc metals with an equilibrium lattice mismatch of only 2%. We neglect this mismatch and use a common lattice constant of a Pt = 0.392 nm in the following. The scattering geometry is constructed as sketched in Fig. 5 with 60 (110) planes each of Pt and Ir stacked in the [110] x direction. It consists of 90 (001) atomic layers sandwiched between ballistic Ir leads in the z direction chosen to be the crystal [001] direction with periodicity of three atomic layers in the y direction. A current is injected in the z direction from the Ir leads, parallel to the Pt|Ir interface.
We average j z (x, y, z) over y and z and plot the re-sultingj z (x) as a function of x across the interface in Fig. 10(a). The horizontal black lines indicate averaged asymptotic values of current densities calculated separately for x ∈ Pt and x ∈ Ir. A transition is clearly visible over a length scale of λ i about the interface. The asymmetry of the transition region with respect to the atomic interface simply reflects the difference between the mean free paths of the two materials. Within the error bars of our calculation the ratio of the equilibrium values of the currents averaged over the Ir and Pt volumes,j Ir z /j Pt z = 1.71, mirrors the ratio ρ 800 Pt /ρ 800 Ir = 1.75 confirming that bulk behaviour is recovered inside the slabs. As shown in Fig. 10(b) where the total charge currents in Pt and in Ir are plotted as a function of z, the current injected from the leads attains its asymptotic distribution essentially immediately.
In Fig. 10(c) we plot streamlines calculated for the charge current in a plane perpendicular to the interface in the Pt|Ir bilayer. Streamlines are parallel to the z axis everywhere suggesting no current flow across the interface in the x direction so each material can be treated as an independent transport channel. A colour map corresponding to the magnitude of the charge current j = j 2 x + j 2 z is shown in the background for reference.

C. Effect of different lead materials
Perhaps surprisingly, neither Fig. 7 nor Fig. 10 provides any indication of current redistribution from the high resistivity to the low resistivity metal, the phenomenon know as shunting. It transpires that this is because these calculations were carried out using the lower resistivity material of the multilayer pair as the lead material so that there is essentially no contact resistance for the low resistivity channel. Because of the mismatch between their electronic structures, there is a substantial 11. Distribution of charge currentjz injected from ballistic Pt leads into a 300 K Pt|Au multilayer (top) and into a 800 K Pt|Ir multilayer (bottom). To obtainjz(x) in (a) and (d), jz(x, y, z) is averaged over y and over the ten central layers in the z direction furthest from the leads. The horizontal black lines in (d) indicate the asymptotic value ofjz obtained by averaging jz(x, y, z) over y and z and x ∈ Pt or x ∈ Ir omitting atomic layers near the Pt|Ir interface where the current varies continuously. In (b) and (e),jz(z) is obtained by averaging jz(x, y, z) over x and y with x ∈ Pt and x ∈ (Au, Ir) separately. The grey curves indicate the corresponding profiles on injecting from ballistic Au (Fig. 7) and Ir (Fig. 10) leads into Pt|Au and Pt|Ir multilayers, respectively. (c) and (f) show streamlines of j(x, z) = jxi + jzk averaged over y and superimposed on a colour contour of j(x, z) = j 2 x + j 2 z in the xz plane. Error bars represent the mean deviation over 10 configurations of thermal disorder.
contact resistance between the high conductance lead and the high resistivity material. To lowest order, this contact resistance is the same as the interface resistance between the two materials making up the multilayer. While it is desirable that the results of calculations with the scattering formalism should be independent of the materials used for the leads, this is not always possible in practice because of present memory and computational time constraints. We illustrate this by considering what happens if we use the high resistivity material of the multilayer pair as lead material.

Au|Pt
Instead of using Au leads to inject a current into the Au|Pt multilayer, we now use Pt. In the Knudsen limit illustrated by the upper panels of Fig. 11, the current density does not saturate in the x direction in either Pt or Au (Fig. 11(a)) nor does it saturate in the z direction in either Pt or Au for the lengths of scattering region considered here (Fig. 11(b)). Compared to the results with Au leads, the current density in Pt is higher whereas that in Au is lower. Because the current does not saturate in the z direction in Fig. 11(b), thej z (z) plotted in Fig. 11(a) was obtained by averaging j z (x, y, z) over y and over the ten central layers in the z direction (lay-ers 41 to 50 out of a total of 90). As a function of x, the anomalous dip at the interface in the current profile plotted in Fig. 11(a) is much larger than in Fig. 7(a) with Au leads.
In Fig. 11(b), the injection of charge carriers from Pt leads into Au (red symbols) is much lower than the injection from Au leads into Au (upper grey line). This can be attributed to the existence of an interface resistance between the ballistic Pt lead and diffusive Au. The converse applies for the injection into diffusive Pt from a ballistic Pt lead which is now higher (blue symbols) than when Au leads were used (lower grey line). Shunting tries to achieve the asymptotic situation where ρ Ptj Pt z = ρ Auj Au z by diverting current from Pt to Au. We see this happening in Fig. 11(b) withj Au z (z) increasing andj Pt z (z) decreasing towards the centre of the scattering region. With our present computational resources, the number of layers in the z direction required to achieve asymptotic behaviour is not tractable for the supercell sizes considered here.
The streamlines plotted in Fig. 11(c) now curve towards the Au layers of the multilayer indicating the flow of current across the interface from Pt into Au. Charge transport in the non-asymptotic case studied here is not amenable to description using a simple parallel resistance model following Ohm's law. It is worthwhile noting that for nanoscale experiments using heterostruc-tures composed of metals with long mean free paths like Cu, asymptotic current distributions cannot be guaranteed for short lengths and the application of semiclassical models may be contentious.

Pt|Ir
We now look at what happens when Pt leads are used in the classical limit for the 800 K Pt|Ir multilayer illustrated in the lower panes of Fig. 11. Comparing Fig. 11 (a) and (d), we see a striking difference between the classical and Knudsen limits near the interface. In Fig. 11(d) j z (x) varies gradually across the interface essentially interpolating between the saturated values calculated previously using Ir leads (indicated in grey) with the ratio of the mean values of the current density in each slab (black lines)j Ir z /j Pt z = 1.52 falling short of the ratio ρ 800 Pt /ρ 800 Ir = 1.71. As we saw for Pt|Au, shunting of the charge current is seen in Fig. 11(e) to vary along the transport direction z but the crossover seen near the leads for Pt|Au is now absent. And just as we found for Pt|Au, the reduction in the current injected into Ir from Pt leads can be attributed to the interface resistance between ballistic Pt and diffusive Ir being larger than the negligible resistance between ballistic Pt and diffusive Pt. Streamlines of the current vector in Fig. 11(f) now curve into the Ir layer of the multilayer indicating a net flow from Pt into Ir.
The above examples show that longer scattering regions are needed in order to realize the situation where the currents reach their asymptotic distributions that are independent of the choice of lead material. This will become increasingly difficult as the temperature is lowered. The layers considered in this study are comparable in thickness to those use in many experiments in the field of spintronics where Pt layers are typically 10-20 nm thick and current distributions are expected to be asymptotic because of the longer lengths of samples. However, by using the high conductivity material in a multilayer as lead material, we showed that it is possible to probe the asymptotic state.

IV. DISCUSSION
We have presented a scheme to study spin and charge currents in nontrivial nanostructures containing surfaces and interfaces that builds upon an extremely efficient fully relativistic quantum mechanical scattering formalism [36,64] and illustrated it with a study of charge transport in thin films and multilayers of nonmagnetic materials. The specific examples that we considered viz.
Pt|Ir and Pt|Au multilayers as well as free-standing thin films of Pt and Au are illustrative of transport regimes where the mean free path λ is either much larger than the thickness d of individual layers or much smaller. The ratio λ/d is the Knudsen number (Kn) that is well known from fluid physics. As pointed out early on by Fuchs [4], Sondheimer [15] and others [32], it plays a crucial role in determining how currents are distributed near a surface where diffusive scattering leads to a suppression of the current in a thin film as confirmed by numerous experimental studies as well as by our calculations for Au and Pt. The same current suppression is apparent in just the Au layer of a Au|Pt multilayer for which Kn > 1 but is much smaller in the Pt layer leading to a current density that varies nonmonotonically when we pass from the centre of the Pt layer through the interface to the centre of the Au layer. For a Pt|Ir multilayer at 800 K for which Kn < 1 there is a smooth and continuous variation of the current density through the interface.
Although the need to more accurately describe current distributions in metallic multilayers, including transient shunting effects, is widely recognized in order to interpret spin-transport experiments [17,[65][66][67], there has been little progress in devising improved methods of doing so [24]. At the same time, in the semiconductor world, there is a growing need to describe electron transport in wires whose size is being constantly reduced [13,68] in order to identify improved interconnect materials. While we performed calculations for systems of O(10 4 ) atoms, accessing the asymptotic limit at room temperature given by Kn > 1 for high conducting materials such as Cu with a long mean free path requires calculations on larger systems ∼ O(10 5 ) atoms. This is currently only limited by computer memory and computational time which will be met by the next generation of computers. In conclusion, our fully resolved current scheme makes it possible to accurately predict electronic transport in the complex geometries frequently encountered in modern microelectronics without introducing empirical parameters.

V. ACKNOWLEDGEMENTS
This work was financially supported by the "Nederlandse Organisatie voor Wetenschappelijk Onderzoek" (NWO) through the research programme of the former "Stichting voor Fundamenteel Onderzoek der Materie," (NWO-I, formerly FOM) and through the use of supercomputer facilities of NWO "Exacte Wetenschappen" (Physical Sciences). R.S.N acknowledges funding from the Shell-NWO/FOM "Computational Sciences for Energy Research (CSER)" PhD program, project number 15CSER12 and is grateful to Max Rang for help with testing the code.