Modelling charge transport in gold nanogranular films

Cluster-assembled metallic films show interesting electrical properties, both in the near-to-percolation regime, when deposited clusters do not form a complete layer yet, and when the film thickness is well above the electrical percolation threshold. Correctly estimating their electrical conductivity is crucial, but, particularly for the latter regime, standard theoretical tools are not quite adequate. We therefore developed a procedure based on an atomically informed mesoscopic model in which ab-initio estimates of electronic transport at the nanoscale are used to reconstruct the conductivity of nanogranular gold films generated by molecular dynamics. An equivalent resistor network is developed, appropriately accounting for ballistic transport. The method is shown to correctly capture the non-monotonic behavior of the conductivity as a function of the film thickness, namely a signature feature of nanogranular films.


I. INTRODUCTION
Cluster assembled metallic (or, simply, "nanogranular") films may play an important role in the development of emerging technologies. In particular, they show a resistive switching behavior [1] that can be exploited in the fabrication of electrical devices able to process and store data in the same physical unit [2][3][4], as requested by the neuromorphic computing paradigm [5,6]. Such behavior emerges in the near-to-percolation regime [7][8][9], when deposited clusters do not form a complete layer of the film yet, as well as when the film thickness is well above the electrical percolation threshold [10][11][12][13]. In particular, for this latter situation a well-established explanation of the underlying physical mechanisms is still missing.
Atomistic simulations may help to get insights on the microscopic mechanisms responsible for such phenomena. Correctly estimating the electrical conductivity becomes therefore crucial [13,14]. To this end, we developed an atomically informed mesoscopic model which provides accurate conductivity estimates for systems composed by interconnected gold nanoclusters.
The conductivity of nanogranular films is strongly affected by the high degree of porosity and defects in the metallic component. As the first step, we resort to atomistic simulations based on molecular dynamics to create realistic structures that capture the complexity of nanogranular films at the nanoscopic scale, a method that has been successfully used in the past to analyse morphology and mechanical properties of such systems [15]. Provided with a realistic representation of the atomic-scale complexity of the system, we proceed by calculating its conductance by means of an Equivalent Resistor Network (ERN). Within such an approach, a system is typically approximated with a network of interconnected resistors whose impedance are determined by its local values; the overall resistance is therefore cal- * mlopez@dsf.unica.it culated using Kirchhoff's circuit laws applied to the network.
The typical scale of the inhomogeneities of a nanogranular film is, however, comparable to that of the electron mean free path, l e , in the corresponding crystalline phase, i.e. l e = 37.5 nm for Au [16]. Specializing on films with thickness well beyond the percolation threshold, we include in the ERN the ballistic component of electronic transport, which dominates at the length scale at which inhomogeneities occur, according to the following picture. Inhomogeneities in the metallic component of a nanogranular film are mainly due to the cluster landing impacts occurring during the deposition stage. They can be characterized as layers of highly disordered (amorphous) matter either between adjacent clusters or within the clusters themselves. We therefore model the metallic component as a collection of amorphous regions mixed with pristine crystalline ones and assume that electronic transport within each region and between regions of the same phase (whose length scale is typically just of a few nanometers) is mainly ballistic, while between regions of different phase the transport is diffusive. Such a picture is encoded in an ERN by requiring that resistors contribute to the overall resistance either ballistically or diffusively, whether they connect regions with same or different phases, respectively. While the diffusive-like behavior is readily obtained by letting the resistors abide by the classical Kirchhoff's circuit laws, the ballistic-like behavior is enforced by assigning to each resistor a value of resistance that does not simply depend on the local structure of the system but on the entire region it belongs to and reflects the size scaling typical of ballistic transport.
Accurate characterizations of the ballistic component of electronic transport in metal nanostructures can be obtained using ab initio methods [17]. In particular, we use density functional theory (DFT) combined with nonequilibrium Green's function (NEGF) techniques (i) to study ballistic transport in structures mimicking the inhomogeneities found in the nanogranular film and (ii) to determine appropriate values of resistance for the ERN. . The corresponding thickness of the sample is < 6 nm and no percolation path connecting two opposite sides of the sample has been created yet; in the middle panel the creation of the first percolation paths is already achieved, corresponding to a film thickness of 20 nm; right panel shows the film in an advanced growth stage with a film thickness > 60 nm.
Provided with such an input, the ERN can be finally used to get an estimate of the conductance of the entire simulated system.
To demonstrate the robustness of our procedure, we have simulated the growth of a nano-sized sample of a nanogranular film assembled by cluster deposition, close to the experimental conditions of Ref. 10, and calculated its resistance at various stages of growth.
The paper is organized as follows: in Section II we describe the methodology used to simulate the growth of the nanogranular Au film by means of classical Molecular Dynamics (MD). In Section III we discuss the procedure to accurately estimate the conductance of the specific gold micro-structures observed in the simulated film. In Section IV the ballistic ERN used to compute the total resistance of the film is presented. Finally, in Section V we present the results provided by our electrical model and compare them to a set of experimental results.

II. SIMULATED FILM GROWTH AND STRUCTURAL ANALYSIS
The cluster assembled metallic film was obtained by simulating the multiple landing of 210 Au clusters deposited in 6 different steps [15], by classical MD. All the gold clusters were first thermalized at 300 K for 150 ps. The size population of the clusters was constructed with 70% of the clusters of diameters 8.8 nm and 30% of diameters 1.3 nm, thus reproducing the experimental sizedistribution [13].
As for the growth process, the first 35 Au clusters were deposited at random positions and normal impact direction on top of the substrate. The average kinetic energies per atom of the landing clusters were fixed to 0.25 eV/atom consistently with the results obtained from the experiments [18]. Periodic boundary conditions were applied in the in-plane directions normal to the growth one. Finally, the clusters were left free to evolve according to Newtonian dynamics. The snapshots of the film corresponding to three different deposition stages are displayed in Figure 1: after the first deposition steps we observe a film characterized by isolated grains, most of the substrate surface being unoccupied. We can observe in the subsequent deposition steps the formation of cavities giving rise to the expected film porosity, i.e. ∼ 30%.
MD simulations have been performed using the LAMMPS code [19], integrating the equations of motion by the velocity-Verlet algorithm. The Nose-Hoover thermostat with relaxation time equal to 100 fs was used to control the temperature. The Au-Au interactions were sampled using a 12-6 Lennard-Jones potential with a cut off at 0.8 nm. The Lennard-Jones parameters have been optimized in order to reproduce several properties such as surface tension density in good agreement with experiment [20], i.e. = 5.29 eV and σ = 2.62904Å. The Au substrate (grey atoms in Figure 1) with dimensions of 24.5×24.5×5 nm 3 was constructed with the (111) surface exposed to the deposition of the clusters. The four bottom layers were kept fixed in order to mimic a bulk material and a slab region (1.7 nm thick) adjacent to the fixed slab was thermalized at room temperature.  Figure 2 displays a gold grain landed on the substrate: while before landing, by construction, the cluster is perfectly spherical, the collision with the substrate strongly affects its shape and structure. We observe in the deposited film two kinds of atomic arrangements: cubic (fcc) and non-cubic (non-fcc) gold. A Polyhedral Template Matching (PTM) Analysis performed with Ovito [21], a scientific analysis software for molecular simulation models, allows to distinguish between those Au atoms sitting in fcc sites and those which are not (related to planar and bulk defects). In the left panel of Figure 2 we show a section of the PTM analysis performed on the gold cluster after landing. For the sake of clarity we have excluded from the analysis those atoms belonging to the substrate. We observe that, due to the collision with the substrate, the fcc symmetry is broken, thus originating regions with different crystal structure. These defects are local and separate different fcc regions (orange colored atoms) within a single grain. In addition to that, the fcc symmetry is also broken at the surface of the cluster creating a non-fcc shell all around the grain. The shell is found before and after the collision, thus, it is not produced by the exceeding kinetic energy after landing. However, we do observe an increase of the shell thickness after cluster landing. Moreover, the nonfcc shell is responsible of having non-fcc layers separating adjacent deposited clusters in the film.
The PTM analysis provides very accurate atomically resolved structural phase maps of the clusters assembled to form the film. However, its heavy computational cost prevents from using PTM to distinguish the different phases during the evolution of the simulated film which counts with more than 10 6 atoms. Therefore, we rather measure the local atomic density, n l , of the film: we observe that the presence of defects induces a slight increase in the Au-Au bond length which shifts from b m = 2.89 A for fcc-Au regions to b m = 2.99Å for non-fcc ones. This effect can be seen in Figure 3(a) where the radial distribution function calculated on a sub-region of the granule containing fcc atoms (orange line) and non-fcc atoms (red line) is shown. A more evident effect is the broadening of the peaks for the non-fcc region. The upward shift of the first peak for non-fcc regions is directly translated to a decrease of the local density that we define as n l = N r /V r , where N r is the number of gold atoms contained in that particular region of volume V r . Two specific local density threshold values, n v and n c , are used to set the density ranges corresponding to vacuum, non-fcc and fcc gold. We consider a region to be vacuum if its local density falls below n v , while cubic gold is defined as n l > n c . Finally, non-cubic gold corresponds to n v < n l < n c . In order to evaluate the agreement between the two methods we have averaged the former PTM analysis over the regions on which the local density approach is performed, i.e. a regular grid with element size 5Å×5Å×5Å, thus V r = 0.125 nm 3 . By doing so we obtain the color map displayed in the middle panel of Figure 2. In the right panel, the corresponding local density map is shown with n v = 0.0005 atoms/Å 3 and n c = 0.048 atoms/Å 3 . Despite the much lower computational cost a good agreement in the ratios between the fcc/non-fcc/vacuum occupied volume and the total volume is observed. A side effect is a slight over-estimation of the non-fcc shell's thickness for the local density approach.

III. ELECTRON TRANSPORT IN GOLD NANO-JUNCTIONS
The ballistic component is supposed to dominate the electronic transport within individual nanoparticles at the considered length scales. Other transport mechanisms such as tunneling and hopping are not included since a strong-coupling regime [42] is expected at the considered temperature and lengths, neither Coulomb interaction and quantum interference effects (see Supplemental Material). We assume the two individuated gold phases, i.e. fcc and non-fcc, to have different electronic transport characteristics. This is justified by the fact that non-fcc Au regions are characterized by a lack of symmetry that effectively reduces the number of opened conduction channels for ballistic transport [22]. Moreover, we also assume that the estimation of the ballistic conductance of the two phases is sufficiently accurate to build up a reliable resistive model. The conductance of fcc and non-fcc Au is estimated by a blended NEGF-DFT approach. In particular, the computation of the conductance of gold nano-sized systems [17,23,24] has been boosted by recent experiments on electronic properties of atomic-sized gold structures [25][26][27][28][29][30][31]. The remarkable agreement between estimates and experimental values for different lengths, cross-sections and crystallographic orientations proves the accuracy of this approach in the study of electron transport in gold systems at the nano-scale.
Using such an approach, we are therefore able to perform a comprehensive study of the conductivity of a twoterminal device containing gold junctions mimicking the structures individuated in the structural analysis of the MD samples, i.e. fcc and non-fcc. More specifically, in order to compute the conductance of gold junctions that might be representative of those found in the simulated film, we proceed as follows: we first set the crystallographic orientation of the device electrodes of the two-terminal device, once for all. Atomic-scale Au junctions are build up with different lengths, L, and crosssections, A, in between the two electrodes. To mimic the fcc phase, we ask the atoms belonging to the central scattering region to keep the crystallographic orientation of the electrodes (as shown in 4(a)), while non-fcc ones are requested to (i) present no specific crystallographic orientation under the PTMA and (ii) to present a broader radial distribution function than fcc as observed in 3(a). This is equivalent to ask to each added atomic layer to change orientation with respect to the previous one. In such a scheme, the first allows electrons to see the symmetry of the lattice along the device, while the latter incorporates the non-homogeneity of the medium found in the deposited clusters. The g(r) for fcc and non-fcc gold junctions are displayed in Figure 3(b). Another possible approach to evaluate the conductance of the defects individuated in the simulated film is to simply carve out from the film those regions we are interested in. The reason to avoid this approach is the fact that the non-homogeneities found in the film extend only for few atomic layers, while the required calculations that allow to specify a unique value for the conductance depending on A and L require the consideration of lengths and cross-sections beyond that limit (see length-scales in Figure 4(b) and (c)).
We describe the gold electronic structure selfconsistently using DFT within the Generalized Gradient Approximation (GGA) as implemented in the SIESTA package [32]. Core electrons are modelled with Troullier-Martins nonlocal pseudopotentials, while the valence electrons are expanded with a double-ζ basis set. The mesh cutoff is 300Ry and a 10×10×10 k-point mesh is used for the 4 atoms unit-cell. We relax all the atomic coordinates till atomic forces are below 0.04 eV/Å and 0.10 eV/Å for fcc and non-fcc Au, respectively.
For the conductance calculations, we have used TRAN-SIESTA [33], which is based on the combination of DFT with the NEGF technique. Therefore, calculations on transport properties are based on the Landauer scheme of elastic scattering probability [34]. Within such a scheme, given a certain bias, V , it is possible to compute the current, I, after self-consistently solving the NEGF and the electrostatic potential to get the electronic density matrix. The conductance of the device is then computed as G = I/V . Semi-infinite 8x8 100-Au electrodes, 4 layers thick, sampled with a converged k-point grid of 3x3x20, are considered.
In the ballistic transport regime the conductance of a material is well known to be independent of the device length. The first step is to compute the conductance, G, for fcc gold junctions for different cross-sections and increasing lengths ranging from tens to hundreds ofÅ.
The computed values of G against the wire length, L, are displayed in Figure 4(b) in units of the quantum of conductance, G 0 = 0.0000775 Ω −1 . We observe fluctuations in the computed G values for short wires (L < 50Å) while G converges to a constant value for longer wires, as expected for ballistic transport. Many experiments [27,35] and computational works [23,24] have reported an increase of G for wires with increasing A, due to the increase in the number of opened conduction channels. A linear relation between G and A is expected with a slope depending on the crystallographic orientation [25]. We obtain values close to 1.2 G 0 and 2.9 G 0 for the two type of gold junctions considered in Figure 4(b). In particular, as the cross-section of a wire oriented along the [1 0 0] direction is increased from 1.7Å (corresponding to 4 unit cells) to 7.0Å (corresponding to 16 unit cells) a increase of ∆G =1.7 G 0 is observed corresponding to the transition from 1.2 G 0 to 2.9 G 0 , close to the expected value provided by the simplified free electron model used in Ref. 25, i.e. 1.8 G 0 . Figure 4(c) shows the conductance for long wires (L = 100Å) for increasing cross-sections for fcc (orange) and non-fcc Au wires (red). The linear trend for fcc junctions is characterized by a slope of 0.04 G 0 /Å 2 while it is less pronounced for non-fcc, i.e. 0.02 G 0 /Å 2 . As expected, non-fcc Au junctions reveal less conductive than fcc ones, for all considered A.

IV. EQUIVALENT RESISTOR NETWORK WITH BALLISTIC TRANSPORT
An ERN model is used to evaluate the electrical conductivity of the simulated film. A 3D regular grid is superimposed to the deposited film and each ERN cell is filled with one of the following: fcc Au, non-fcc Au or vacuum depending on the local value of n l as explained in Section II. A xy projection of the grid over a single landed gold cluster is shown in Figure 5(a). We assign to each ERN grid element a resistance value as follows: for vacuum elements this is set to 10 15 Ω, while for Au elements, this is computed sticking to the geometrical dependencies found in Figure 4(b) and (c). The rational behind the value chosen for vacuum is to assure those elements do not contribute to the final resistance value. This can be achieved by setting this value to infinite, which carries numerical issues. Instead, we choose to set this to a sufficiently high value with respect to the fcc and non-fcc ones, both falling in the kΩ range.
From the simulated sample we can roughly distinguish two types of interface: interfaces separating two grains with different crystallographic orientation ("fcc-fcc") and those separating ordered grains from disordered regions ("fcc-non-fcc"). In our modelling, the latter type is always recognised as an interface that disrupts the electronic transport so the total resistance of a slab of fcc gold in contact with a non-fcc region is equal to the sum of their resistances, capturing the decoherence of electrons when reaching the interface and the interruption of the ballistic transport. In other words, the interface limits the regions where the electronic transport is considered ballistic. Instead, fcc-fcc interfaces, on the other hand, affect the transport only if the change of symmetry in going from one grain to the other is high enough so the density analysis detects the intermediation of a non-fcc region in between, otherwise the interface is effectively neglected. It must be remarked that clean fcc-fcc interfaces rarely occur as one can see in Figure 2.
Thus, for a given mesh element containing Au atoms in a fcc (non-fcc) phase, we compute the number, (N x , N y , N z ), of consecutive cells along each cartesian direction containing cubic (non-cubic) gold. We then give a unique resistance value for each transport direction as ] −1 and i = x, y, z. Doing so we assure the total resistance of a given chunk of gold does not depend on its length, and only the cross-section A determines its final value as expected for ballistic transport. Therefore, given a transport direction, if the total resistance value for a single fcc(non-fcc) mesh element is R c (R a ), the total resistance for N consecutive fcc(nonfcc) cells along that transport direction equals R c (R a ). Instead, the alternative stacking of fcc and non-fcc elements result in a total resistance that equals the sum of the consitutive parts, as represented in Figure 5(a).
The ERN grid counts with I×J×K elements the dimensions of which have been chosen so as (i) to minimize the computational cost of solving iteratively the ERN for the considered structures, and (ii) to have enough spatial resolution to well resolve the intra-and inter-granules structure. We set the element size to 5Å, so the ele-ment volume and the minimum resolved area are 0.125 nm 3 and 0.25 nm 2 , respectively. With this the ERN has 50×50×K elements with K increasing at each deposition step in order to include all deposited clusters. Once all the mesh elements count with a resistance value, a finite bias, V bias , across the sample is applied: the voltage is set to V = V bias for those grid elements belonging to one of the electrodes, while it is set to V = 0 V otherwise. We have used V bias = 0.06 V to generate all the data included in Section V. We stress at this point that the R T estimation of the total film resistance produced by the present linear model does not depend on this parameter.
The obtained electrical network is analyzed by solving the Kirchhoff equations. We solve them iteratively updating the node voltages V ijk using the formula where R ijk,x is the resistance of the (i, j, k) grid element in the x direction, etc., and keeping fixed the electrode voltage. Iterations are performed until the variation of the sample total resistance between iteration steps is less than 0.01 Ω.
From the node voltages and the resistances, one can calculate the current flowing through the simulated sample. Each (i, j, k) grid element counts with a three component current vector {I ijk,x , I ijk,y , I ijk,z }. In Figure  5(b) we show the converged current map at z = 15 nm (k = 30) for the simulated sample at a very advanced growth stage along with the corresponding fcc/non-fcc grid at that film height. The total current vector along the bias direction (x direction in 5(b)) is computed as I = ( jk I ijk,x , jk I ijk,y , jk I ijk,z ). The three components of the total current vector are displayed in 5(c). We observe that, after reaching the converged solution of the ERN, the x component of the total current vector, I x , fluctuates around a constant value (1.5 mA with V bias = 0.06V ) all along the transport direction (black line in Right panel of Figure 5(c)). For I y and I z we observe the generation of internal currents that cancel each other (see ∼ 0 values close to the extremes of both curves) so no current gets in or out the system through the directions normal to the transport one. Finally, we compute the total resistance of the film as R T = V bias /Ī whereĪ is the mean value the total current flowing through the film along the bias direction.
To reach the converged solution for the ERN represents the most expensive part of the present model in terms of CPU time (600 seconds for a 50 × 50 × 75 ERN run in a single core).

V. RESULTS
The hierarchy of MD, NEGF-DFT and ERN models allows to compute the evolution of the film electrical resistance as gold clusters are deposited on the substrate. We have reproduced the percolation curves for the simulated sample and the corresponding data are displayed along with the experimental measurements from Ref. 10 in Figure 6. We considered two different electrostatic bias conditions by setting the electrodes along the x (x-bias) and y directions (y-bias), as depicted in Figure 6. The determination of the film thickness, t, is performed following the definition used to plot the results of the experiments in Ref. 10: we count the number of cells containing Au atoms in the ERN grid, N c . Given the sample in-plane dimensions we compute t b =0.125nm 3 N c /(24.5nm) 2 corresponding to the thickness of a bulk film with the same amount of matter. Finally, the porosity of the sample is introduced to obtain the final thickness value, t = 1.35t b . We name t p the film thickness at which the first percolation path is created. For t < t p , the film is characterized by isolated clusters and a infinite resistance (> 10 13 Ω for numerical reasons). Since the current implementation does not take into account electron tunneling and hopping effects, our data starts being meaningful after the creation of the first percolation path, which creates a real contact between the the electrodes. We observe the first percolation path occurring at t p,x = 5 nm and t p,y = 10 nm, for the two considered bias conditions. After that, two growth stages are identified from the resistance-thickness curve [36]. At first, few inter-grain electrical contacts exist producing a film characterized by poorly connected aggregates and R T values in the 1-10 kΩ range. In this stage, known as geometrical percolation stage, R T abruptly decreases down to hundreds of Ω due to the increase in the paths becoming available for electron transport as clusters land and interconnect. Next, a transition from insulating to ohmic behavior is observed, i.e R T smoothly decreases as t is increased. This transition is defined as the thickness at which the quantity R T t 2 reaches its minimum value t m [37,38]. The good agreement achieved for the determination of this parameter is shown in Figure 7. For all data sets the transition to ohmic behavior is achieved around t m = 20 nm. However, this parameter can suffer huge fluctuations from one sample to another [10]. The predicted R T values show larger fluctuations in comparison to the experimental curve, these being more important at the first percolation steps. The fact that the simulated film has nano-scale dimensions, i.e. 25×25nm 2 , make granular effects still visible at all deposition stages and the landing of a single cluster produces huge variations in the computed values. Instead, the experimental curve corresponding to a macro-scale sample, presents a smoother behavior.
For t > t m , the total film resistance further decreases (see Figure 6). Quantitatively, our electrical model slightly underestimates the total resistance of the film for every deposition step. For instance, the ratio between the computed and the measured value at t = 40 nm is 0.6 and 0.7 for x-bias and y-bias, respectively. Both experimental and simulated data follow a power law decay (R T ∝ 1/t α ) as clusters are deposited and this represents a qualitative agreement with the experimental data set. The experimental curve provides α exp = 0.9, implying that the resulting film resistivity, calculated as ρ = R T t, should increase instead of remaining constant as expected for bulk materials and also observed for atomically assembled Au films. The power-law exponents for the two simulated curves equal to α x = 0.94 and α y = 0.95 for that thickness range, fulfilling the condition α < 1 to observe an increase of ρ. Indeed a quite remarkable agreement with experimental results. This effect is often attributed to an increase of the sample roughness as suggested in Ref. 39-41. In order to support the plausibility of this hypothesis we have computed the roughness of the film surface, Rq, which we estimate by means of the standard deviation of the film height measured on every (i, j) mesh element as where z ij is the height profile of the sample andz stands for the mean sample height. Differently to experimental realizations, we do have access to Rq at every deposition step so a sound characterization of the film roughness dynamics is achieved.
In Figure 6(b) the evolution of Rq is shown, displaying clearly two alternate regimes: one in which clusters are accumulated in few spots of the (x, y) plane producing an increase of Rq corresponding to t < 40 nm and t > 45 nm, and another in which the deposition of few clusters are enough to drastically reduce Rq generating void regions within the film which results in the expected film porosity. The former, characterized by a non-uniform deposition of clusters, produces a reduction of the slope of the R T curve, due to the fact that the clusters do not contribute to create new conduction paths. The latter, in which new clusters occupy regions connecting separated gold regions, effectively reduces the film resistance. The periods of these two regimes are expected to depend on the ratio between the sample size and the clusters dimensions, being shorter for smaller ones. The correlation between the modeled film resistivity and its roughness for t > 40 nm can be observed in Figure 8. Thus, the present electrical model well captures the role of surface roughness in the evolution of the film growth.

VI. CONCLUSIONS
We have presented an electrical model able to predict the resistance of cluster-assembled gold films based upon a well resolved description of their nanostructure and related charge transport at that length scales. By reproducing the deposition of tens of Au clusters on a substrate by means of MD simulations we are able to reproduce the growth of a Au nanogranular film. Next, by using our model at different growth stages we compute the resistance percolation curve and compare the model predictions against experimental values obtaining a good agreement. In particular, the power-law dependence of the total resistance with the film thickness is well reproduced as well as the non-monotonic behavior of the film resistivity and its dependence with the film surface roughness. The modeled film is a ∼30% more conductive respect to the experimental one, likely due to the simplified synopsis of transport mechanisms it is based on. Further improvement of the model could be achieved by including tunneling and hopping transport mechanisms so to extend the applicability of the model to the t < t p region of the percolation curve.

TUNNELING AND SINGLE GRAIN CONDUCTANCE
In modeling electron transport in granular materials, one can distinguish different transport regimes using two key quantities [1]: the intergranular conductance g and the intragrain conductance g 0 , both expressed in units of the quantum conductance e 2 /h. In the following we will provide an estimate of those quantities and prove that our system is in a conducting ( g 0 1) and strong grain-coupling (g 1) regime, in which effects such as Coulomb blockade can be excluded.
The intragrain conductance can be estimated as g 0 = E th /δ [1], where E th is the Thouless energy and δ the mean energy level spacing for a single grain. Under the simplifying assumption that our system is made of Au clusters of diameter d = 8.8 nm (instead of two distinct populations with d = 8.8 nm and d = 1.3 nm), the mean energy level spacing can be estimated to be δ ∼ 20 meV [2], while E th = l e v F /d = 800 meV , where we have considered an electron mean free path l e = 40 nm and a Fermi velocity v F 1.4 · 106 m/s, both at T = 300K. Combining those values one obtains g 0 = 40, which indeed is much greater than the unity and proves that are Au clusters are indeed good conductors.
The intergranular (also called "tunneling", even when the regime is not of quantum tunneling) conductance g is harder to estimate from the available information on the experiment only. We therefore resort to our simulation and ab initio calculations to provide an estimate of the conductance of junctions as the one represented in the following picture: where the outer rectangles identify the electrodes and the inner hemispheres the interface between two nanoparticles. In a granular sample composed by mono-sized deformable but incompressible (i.e. volume conserving) spheres, the average surface of contact between two grains S avg is only function of the occupied volume fraction Φ. Considering Φ = 0.675, which is extracted from our simulated sample, we can estimate S avg = 15 nm 2 , based on systems with similar geometrical characteristics [3]. Given S avg the conductance of the corresponding junction can be extrapolated from our ab initio calculations for similar junctions. In particular, using the data for non-fcc interfaces shown in Figure 4 of the manuscript, we compute a tunneling conductance of g ∼ 30, from which we can deduce that our system is in the strong coupling regime. As the ratio between grain and tunneling conductance is g 0 /g ∼ 1.3, we can say that our system is so strongly coupled that can be considered as a porous metal.

ELECTRON LOCALIZATION
The inverse of the Thouless energy 1/E th = 1.25 eV −1 gives an estimate of the time that it takes for an electron to traverse the grain. By comparing that with the time takes an electron to pass from a grain to the next one, which can be estimated as (gδ) −1 = 1.7 eV −1 [1], we can say that electrons are effectively delocalized (they time they spend in a grain is roughly the time they need to go from one edge to the other) and effects such as Coulomb blockade can be excluded.

SAMPLE CONDUCTIVITY FROM GRANULAR MODEL
The conductivity of dense nanoparticle systems can be estimated using a model of ordered arrays of conducting, strongly coupled grains, which can be solved analytically [4,5]. Using the parameters discussed above, one obtains σ ≈ 5.246 × 10 5 Ω −1 m −1 . To see how this compares with our estimate, we superimpose the corresponding resistance arXiv:2110.02572v2 [cond-mat.mes-hall] 16 Oct 2021 to Fig. 6 of the manuscript, as follows:

Experiment
Our calculation (X-bias) Our calculation (Y-bias) σ from model in Ref. [6] We observe an underestimation of the resistance that is slightly stronger than our own estimate. We would therefore argue that this can be regarded as a further validation of our approach.
The model includes corrections to the classical Drude conductivity for a granular metal due to Coulomb interaction, which here are found to be negligible. More specifically, in writing the conductivity as σ = σ 0 + δσ 1 + δσ 2 , where σ 0 is the Drude conductivity, δσ 1 and δσ 2 the correction due to the contribution from the large and low energy scales, respectively, we found that σ 0 = 5.25×10 5 Ω −1 m −1 , δσ 1 = −474.2 Ω −1 m −1 and δσ 2 = 56.3 Ω −1 m −1 , which shows that these corrections are both smaller than σ 0 by orders of magnitude. We can therefore deduce that Coulomb interaction and quantum interference can be neglected. Finally, we notice that the Coulomb energy (needed to compute σ 1 ) E c = e 2 /(2C), with C = 4π 0 r d and r = 6.9 for Au, gives an estimate of E c ∼ 12 meV , which is lower than the thermal energy 26 meV .