Modeling Resistive Switching in Nanogranular Metal Films

Films produced by assembling bare gold clusters well beyond the electrical percolation threshold show a resistive switching behavior whose investigation has started only recently. Here we address the challenge to charaterize the resistance of a nanogranular film starting from limited information on the structure at the microscopic scale by the means of Bruggeman's approach to multicomponent media, within the framework of Effective Medium Approximations. The approach is used to build a model that proves that the observed resistive switching can be explained by thermally regulated local structural rearrangements.


I. INTRODUCTION
The synthesis of materials and devices with new properties by means of controlled manipulation of their structure up to the atomic scale is nowadays paving the way to a plethora of new technological applications. Among these, neuromorphic computing architectures 1 , which are believed to have the potential of overcoming the looming end of Moore's scaling law, are becoming possible thanks to 'resistive switching junctions' 2 , that is, devices that suddenly change their resistance under the action of an electric current in a non-volatile and reversible way.
Recently, cluster-assembled gold films have been reported to present an unexpected form of resistive switching 3,4 . By supersonic cluster beam deposition of bare Au nanoparticles on various substrates, metallic films characterized by a complex microstructure have been grown well beyond their electrical percolation threshold. Nanoparticles deposited on the substrate basically retain their individuality, thus forming films of (nano)granular matter (henceforth "ng-films"), a medium highly porous and rich of interfaces among misoriented crystal grains. Besides increasing the electrical resistance of ordinary (i.e. continuous or, equivalently, atom-assembled) thin films, their nanostructure results in a remarkable dynamical response to a sufficiently strong external electric field. For purely metallic films, such a feature has been observed until then only in highly discontinuous films with thickness near the percolation threshold 5,6 . Even under the action of a direct electric current a ng-film presents a resistance characterized by nearly ohmic regimes (during which only small fluctuations occur) and abrupt jumps to higher as well as to lower values, after which the system either returns to the same ohmic regime (spikes) or reaches a new one (steps). Such a dynamical non-ohmic behavior in a metallic system reflects most likely some sort of microstructural change that needs to be investigated experimentally as well as theoretically. This work represents an attempt to build some appropriate theoretical tools for such an endeavor.
The outcome of a complete theoretical characterization of the resistive switching phenomenon would be a model able to reproduce the abrupt changes in resistance and to provide a reasonable estimate of the spanned values. The challenge is twofold. On the one hand, one needs to identify and characterize the underlying mechanisms responsible for the jumps of resistance. On the other hand, one must be able to estimate the resistance of a ng-film, which, even in the ohmic regime, is not at all an easy task to accomplish, as the granular nature of the system is not captured by the standard models used to estimate the resistance of crystal films. The aim of the present work is therefore to establish a theoretical framework that allows to estimate the resistance of a ng-film starting from ingredients which could be accessible to ab initio and semiempirical methods. In other words, the present attempt to elaborate an analytical formal theory linking the observed resistive switching behavior to some underlying microstructure evolution should pave the way of a multiphysics approach to this intriguing phenomenon, where the theoretical model will be fed by physical information computed from first principles.
The first steps in that direction were made already in Ref. 3, where the experimental results were qualitatively reproduced using a dynamical resistive network model. The inherent simplification in representing the articulated nanoparticle interconnections with a net of resistors and the demanding computational cost of solving the corresponding Kirchhoffs equations make, however, the approach unsuitable to model realistic systems. We therefore adopt an Effective Medium Approximation (EMA), which allows to accurately estimate the electrical resistivity of a ng-film starting just from little knowledge about its microscopic structure. In Section II we present the approach and explain how, in a steady situation, it can be used to account for the high degree of porosity of a structure, as well as for other structural features that define the nanogranular character of the system. In Section III we add a dynamical evolution due to thermally regulated structural changes, and show that they indeed rule over the observed resistive switching behaviour; we also explain how they can be modeled within our framework. Two realizations of this modeling are shown to capture many features of the experimentally reported ngfilm resistance 3,4 .
Such analysis allows us to trace a clear path for future investigations aimed at identifying the actual atomicscale mechanisms occurring in the films. Full extent and limits of the approach are discussed in the last section.

II. OHMIC REGIME
A. Synopsis of the residual resistivity contributions The residual resistivity (RR) of ng-films can be orders of magnitude higher than that of films of the same metal in the crystalline phase 7 . In fact, the RR of these latter can already differ quite sensibly from that of the bulk material. Grain boundaries and surface effects are known to play a major role 8 and the models of Mayadas-Shatzkes 9 and Fuchs-Sondheimer 10,11 are routinely used to account for them.
In the case of ng-films, additional structural imperfections have a higher impact on the overall RR. We organize such imperfections in the following way. Foremost we consider the contribution of voids, namely the interstitial vacuum existing among nanoparticles. Although they do not allow for band transport, they still admit hopping and tunneling, whose contribution to the overall current is, however, expected to be much lower than band transport. A second relevant source of RR are the interfaces between nanoparticles. In the best case, they can be assimilated to boundaries between randomly oriented grains belonging to different nanoparticles. Most likely, however, an amorphous layer is present in the contact region, as a result of the impact during the deposition stage 12 . To best of our knowledge, there is no wellestablished method in literature to account for neither imperfections in estimating the resistivity of a film.
Surface effects, which may become relevant for thin films, can be estimated using the Fuchs-Namba model 13 , an extension of Fuchs formalism 10 that also takes into account the roughness of the top layer of nanoparticles. The presence of grains, on the other hand, is not adequately described by the model of Mayadas-Shatzkes, which assumes that grains grow in a columnar structure, with the axes normal to the film plane, and extend from two surfaces of the film. While usually valid for continuous crystal films, these assumptions are not justified for nanogranular ones, where grains exist inside the deposited nanoparticles and are therefore different in shape and orientation from those considered in the Mayades-Shatzkes model. Impurities and further defects inside the grains add the last contribution to the RR.
Such a hierarchy of contributions to the RR of a nanonagranular film is summarised in Fig. 1.

B. EMA approach
In estimating the electrical resistivity of inhomogeneous media, EMAs have proven to be particularly   [14][15][16] . The underlying principle is that the inhomogeneous medium is effectively treated as a homogeneous one whose properties are calculated in a meanfield approximation. For regimes well above the percolation threshold, as those we are here considering 3 , such simplification allows to accurately describe very complex systems with simple analytical formulas.
When speaking of EMAs one usually refers to the Maxwell Garnett formalism 17 , the Bruggeman formalism or one of their extensions. We here adopt Bruggeman's approach 18 , which, contrary to Maxwell Garnet one, allows to deal with a mixture of several components treating them all on the same footing (see Fig. 2). In such framework, the resistivity of a multicomponent system is estimated as the inverse of the "effective conductivity" On the left, the real system, which is composed by a mixture of two components characterized by their electrical conductivity. On the right, the model, in which a single region (of either components) is assumed to be a sphere embedded in a single medium of uniform conductivity. As an external field is turned on, the presence of the sphere causes a certain polarization. By imposing that the average polarization due to all regions cancels out, one deduces the relation Eq. (1), which allows to determine the conductivity of the medium σe in a self-consistent manner. σ e , the positive root to the equation 14 where i runs over all the components, Φ i is the relative volume fraction occupied by the component i, for which i Φ i = 1, and σ i is the corresponding conductivity. A first, crude estimate of the resistivity of a ng-film can be obtained by considering the film as a mixture of bulk gold, for which σ Au = 0.041 nΩ −1 m −1 at room temperature 19 , and insulating vacuum, σ V = 0. Assuming that the nanoparticle density in the film is close to that of a set of randomly packed spheres, we consider Φ Au ≈ 0.63 20 , and hence Φ V = 1 − Φ Au ≈ 0.37, in Eq.
(1) to estimate the resistivity of a ng-film to be ρ ≈ 55 m nΩ, which is twice as much that of bulk gold. A more accurate estimate can be obtained if we improve the description of the gold component. For instance, we can include the effects of grain boundaries by taking its conductivity to be that of a polycrystalline film with similar grain size. Using σ Au ≈ 0.01 nΩ −1 m −1 , which is the conductivity of a film with grains of approximately 10 nm (see 21 table 1, sample S1), the resistivity of our ngfilm is estimated to be ρ ≈ 225 m nΩ. This value, which is about ten times that of a crystalline film, is comparable with the values reported by Mirigliano et al. of ρ ∼ 100 ÷ 1000 m nΩ 3 , which suggests that our approach is sound.
The flexibility of the EMA approach, in fact, allows us to be even more accurate. Suppose that within the highly inhomogeneous metal component we can identify regions of relatively uniform conductivity, which we say belong to the same "phase": Eq. (1) can be accordingly reinterpreted as referring not to a mixture of just two components (vacuum and metal), but to a mixture of different phases, which can also belong to the same component (in our case, the metal one). The classification of the various inhomogeneities of a ng-film exposed in the previous section helps us in identifying phases in our systems. For instance, if we believe that amorphous layers at the interfaces between nanoparticles account for an important contribution to the total resistivity, then one can single that contribution out by introducing three phases: a vacuum, an amorphous and a polycrystalline phase. If, on the other hand, one believes that nanoparticle interfaces do not contribute more than simple grain boundaries, but the size of nanoparticles is more important, one can introduce different phases to describe nanoparticles of different size. For instance, in Ref. 3 it is reported that two distinct nanoparticle populations characterized by a different radius (0.7 and 4.4 nm) are clearly discernible, at least in the initial part of the deposition. One can therefore introduce a phase for the population of smaller nanoparticles and one for the larger ones. A resistivity of about ρ ∼ 1000 m nΩ is obtained, for instance, by considering four phases: one for the vacuum, for which σ 0 = 0 and Φ 0 = 0.35, a first metallic phase, which might represent large nanoparticles, with σ 1 = 0.005 nΩ −1 m −1 and Φ 1 = 0.26, a second, less conducting metallic phase, representing small nanoparticles, with σ 1 = 0.001 nΩ −1 m −1 and Φ 1 = 0.26, and a third, even less conducting metallic phase, representing amorphous layers between the nanoparticles, with σ 1 = 0.0005 nΩ −1 m −1 and Φ 1 = 0.13.
The EMA approach here adopted is expected to be accurate for low concentrations of the insulating component, while problems can arise near the percolation threshold 22 . Even though we are interested in films with thickness well above the percolation threshold, we would like to test the performance of the approach in the nearpercolation regime to assess its range of validity. To do that, we use Eq. (1) to estimate the conductance of a ng-film in the growing stage, before the onset of resistive switching, and compare it with the data reported in Ref. 3. During growth, three stages occur: i) nanoparticles start to be deposited and, until no percolation paths are formed, the sample conductance is zero; ii) first percolation paths are formed and a whole first layer of nanoparticles is gradually filled; iii) the deposition continues on top of that layer creating a film of higher and higher thickness but effectively constant density. As the concept of "film thickness" makes sense only in the third stage, in Ref. 3 a percolation curve was obtained by monitoring the mass of the sample and plotting its conductance as function of the "nominal thickness", i.e. the thickness of a film characterized by the sample mass and the final film density, length, and width. This choice causes a "shoulder" (point of non-differentiability) of the curve predicted by the EMA occurring when the nominal thickness hits the height of the first layer, which marks the beginning of the third stage and the end of the near-percolation regime. In the third stage, the EMA predicts indeed a constant conductivity, as one would expect for thick enough films. In a simple scenario in which only two phases, vacuum and metal, are considered, the specific EMA percolation curve for the ng-film considered in the experiment is obtained by fixing the value of four parameters (phase conductivities and relative volume fractions). As before, we assume that the vacuum is insulating and the nanoparticles fill around 65% of the volume of the completed film. The remaining two parameters are then tuned in such a way that the EMA percolation curve matches the experimental data at the end of the near-percolation regime, namely at the shoulder. In Fig. 3 we report a comparison between the measured and the EMA percolation curve in the nearpercolation regime, in which the shoulder has been identified, and matched, with the data at x = 30 nm. As we can see, the predicted behavior in the near-percolation regime is quite close to the observed one, with a very accurate estimate for the percolation threshold. Improving such an agreement has resulted particularly hard. Surface effects, which we estimated using Fuchs model 8,10 using a wide range of values for the unknown electron mean free path (from 40 nm, the value in Au bulk 23 , down to 4 nm, the average size of the smallest nanoparticles deposited), were found to be negligible. Also introducing hopping and tunneling processes via a more sophisticated EMA approach 24 did not improve the agreement. In fact, we argue that, as the data refer to a single sample and not to a statistical average, features like the two bumps around 17 and 20 nm of thickness may pertain to the specific sample considered and not be general features of ng-films. We are therefore prone to consider Eq. (1) sufficiently accurate for the ng-films considered even in the near-percolation regime.

Suggested Explanation
Breaking/recreation of junctions An active spot is identified with the interface between two nanoparticles (top left). Such system makes a transition to a higher conductivity state via defect migration (top right); from that, a lower conducting state is reached because of local melting (bottom); finally, the original interface is rebuilt (top left) by action of the original crystals in the nanoparticles that do not participate to the transition (inactive spots, which surround the active ones).

III. MODELING DYNAMICAL PROCESSES
Inspired by the mechanism proposed in Ref. 4, here recalled in Fig. 4, we consider small regions of the ngfilm (say a few nanoparticles) that undergo some cyclic structural change. We call these regions "active spots" and assume they are randomly distributed within the ngfilm. Starting from a given initial atomic configuration, an active spot experiences over time a rearrangement of its atoms, as effect of temperature variations possibly due to fluctuations of the current flowing through it. A change can either lead to a configuration with higher conductivity, if, for example, two nanoparticles coalesce or a defect migrates, or to a lower conductivity, as effect, for instance, of local melting. After several changes, the active spot can be found again in its initial configuration.
The presence of an onset time of a few seconds for the jumps to start suggests that in this transient regime the spots are not active yet. We argue that a current is needed for the Joule heating to provide the missing energy to activate the spots. More specifically, we conceptualize the process sequence as follows: the cycle is initiated as soon as the temperature of the spot exceeds a critical temperatureT : when this situation is reached, the active spots have sufficient energy to overcome the activation barrier that prevented the first rearrangement to occur. In fact, it is reasonable to expect that each of the steps of the cycle is affected in one way or another (i.e. activated or deactivated, at least to some degree) by the local temperature. We therefore suppose that a few critical temperaturesT 1 ,T 2 , ..., which mark different energy thresholds that allow or forbid the relevant transitions, exist. The actual evolution of an active spot is determined by current-related fluctuations of the local temperature. In order to describe the overall dynamical evolution, we adopt a stochastic approach. More specifically, we consider a discrete time variable t a , characterized by a step δt = t a+1 −t a that can be tuned to match the experimental sampling frequency (e.g. δt = 0.1 s in Ref. 3). The dynamics of the system is determined by a rule that dictates how each single active spot changes its state from one timestep to the next. We assume a Markovian evolution, according to which at each timestep an active spot has a given probability of changing its state that only depends on its current state and the ng-film temperature. A complete characterization of the system is therefore achieved by specifying the corresponding transition probability matrix P . If, for instance, we assume a thermally regulated three-step cycle, the transition probability matrix can be written as where the only three degrees of freedom of the matrix are parametrized by functions of the ng-film temperature p nm (T ) for which 0 ≤ p nm (T ) ≤ 1 and are taken to suddenly change around some critical temperaturesT 1 , T 2 , ... Such parametrization takes into account the fact that transitions can only happen in one direction (so, say, P 1→2 > 0 but P 1→3 = 0) and ensures that total probability P n→1 + P n→2 + P n→3 sums up to 1. The functions p nm (T ) are characteristic of the microscopic mechanism one wants to model and their exact form can in principle be determined by means of numerical simulations of a single active spot. As they are function of the temperature, we must also specify how that quantity evolves with time. Temperature variations of the ng-film happen for two reasons, namely i) Joule heating, for which V being the applied voltage, R the resistance of the ngfilm, m its mass and c its specific heat, and ii) heat dissipation due to contact with the environment, for which with τ a characteristic time constant that describes how fast the film heat is transferred to the substrate and T env the environment temperature. Combining the two contributions into a law for discrete time propagation we can write where c 1 and c 2 are two suitable constants and the film conductivity σ(t a ) is calculated via a time-dependent version of Eq. (1), namely where the index i runs over all the phases identified in the system while a represents the time index. If the temperature of the sample spans a wide range of values, the conductivity of the single phases may vary due to the increasing contribution of electron-phonon scattering events, in which case we must consider instead of Eq. (7). Around room temperature, the resistivity of gold, from simple bulk 25 to very complex systems like nanocrystalline films 26 , increases linearly with the temperature and therefore we consider where α i a positive parameter usually called Temperature Coefficient of Resistivity, TRC. A signature of a nonvanishing TRC is a slight drift in the ohmic regime towards higher values of resistivity, which can be recognized in the resistance measurements of a ng-film under the bias of 0.5 V in Ref.

3.
Provided with numerical values for conductivity σ i(0) and the TRC α i of the different phases, the parameters c 1 , c 2 of the temperature equation (6), a specific behavior for the functions p nm (T ), the total number of spots and the phase initial densities, it is possible to propagate the model and simulate the time evolution of a ng-film. In Fig. 5 we present a specific realization of such model, whose parameters, reported in Table I and II, were chosen to give rise to the relevant features of the resistive switching behavior observed in the experiments 3 . We can indeed clearly distinguish a big spike and jumps between ohmic regimes characterized by different duration, fluctuation amplitude and resistance value, which can be higher as well as lower than the initial value. The ohmic regimes are characterized by fluctuations due to the occurrence of transitions of the active spots. The amplitude of the fluctuations reduces if the number of active spots is increased. A comparison between the temperature and resistance evolution shows that steps occur every time a critical temperature is crossed. In other words, whenever the transition probabilities are stable, only fluctuations occur, while jumps arise when transition probabilities suddenly change. The duration of a ohmic regime is therefore linked to the time the sample takes to reach a critical temperature. As in the experiments, a voltage threshold, under which no jumps but only fluctuations occur, is also observed, as result of the fact that Joule heating, which is proportional to V 2 , is not sufficient to take the ng-film to cross the lowest critical temperature. Rising the voltage well above the threshold can significantly accelerate the dynamical evolution, as critical temperatures are reached in a shorter time.
A statistical analysis allows to exactly calculate the values of resistivity spanned during the ohmic regimes starting from the transition probability matrix defined in Eq. (2) and the values of conductivity of the single phases. During such regimes, the temperature of the ngfilm is between two critical temperatures and the transition probability matrix is constant, P (T ) = P . We can then say 27 that the average concentration of the different phases Φ i is given by Φ i = P ∞ i1 where P ∞ is the limiting transition matrix P ∞ = P P P...., for which P ∞ 1m = P ∞ 2m = P ∞ 3m . Provided with Φ i for each ohmic regime, one can use Eq. (1) to calculate the effective conductivity, and hence the resistance of the film, for the such regimes without performing a propagation of the stochastic equations. Such an analysis is particularly important for cases in which the conditions of the simulation, or, for that matter, of an experiment, cannot ensure that the entire landscape of allowed values is spanned.
Finally, to prove that the approach can quantitatively, and not only qualitatively, explain the experimental data, in Fig. 6 we show another realization of the model, whose parameters were tuned to match the experimental results for the ng-film with thickness 30 nm under the bias of 0.5 V, reported in Ref. 3,Fig. 4(c). In fact, since the same graph R(t) can be obtained with various choices of the parameters, we assumed certain features of the spots that one would expect from the specific mechanisms represented in Fig. 4. More specifically, we assumed a certain ratio between the conductivity of the different phases and an overall behavior of the functions p nm (T ) characterized by two critical temperatures, one corresponding to the defect migration activation temperature and the other to the local melting temperature. Given then the specifics of the experiment (sample size and density, and frequency of the current measurements), we tuned the two parameters of the conductivity of the first phase σ 1(0) and α 1 , which set a scale for the resistance, the constant c 2 , which determines the time scale of the problem, and the number of active spots, which sets the amplitude of the fluctuations in the ohmic regimes, to match the corresponding features of the experimental data (cf Tables  I and III). As a result, the experimental data are entirely reproduced with a high level of accuracy, thus proving that our approach is indeed suitable for such kind of experiments.  Fig. 5 and Fig. 6. While for the latter fully dimensional values were used, for the former the parameters are expressed in terms of three quantities (one for the conductivity, one for the temperature and the timestep) that set the scale of the problem.   Fig. (6).

IV. CONCLUSIONS
Nanogranular gold films with thickness well beyond the percolation threshold present an intriguing, unexpected dynamical response to external electrostatic potentials. Understanding and, ultimately, harnessing such phenomenon can lead to interesting technological applications. Their modeling, however, poses some theoretical challenges.
Even before introducing dynamical effects, estimating the resistivity of a nanogranular film eludes the standard models used for crystal films which do not take into account the complex structure arising from the presence of nanoparticles. In this work we suggest to approach the problem using Bruggeman's formalism for multicomponent media, which belongs to the framework of the Effective Medium Approximations. The approach was designed for macroscopic scales and its use to nanometric scales might not be fully justified 15 ; moreover, it is expected to be accurate only when a sufficiently large number of percolation paths have been established 22 , a regime that can not always be easy to establish by only looking at the overall resistance of a discontinuous system 28 . Nonetheless, apart from those restrictions, the method is typically rather accurate, computationally very inexpensive and, most importantly for our case, quite flexible. We have indeed explained how it can be adapted to describe not only the porosity, but also other features of ng-films, such as amorphous layers at the nanoparticle interfaces and the presence of nanoparticles of different size, complementing other models designed to account for the remaining sources of RR. Although more sophisticated approaches within the EMAs framework may be considered 24 , we provided evidence that the one we adopted can be accurate enough even near the percolation threshold.
Provided with such a tool it is possible to build dynamical models. To give a concrete example of the flexibility and extent of the approach, we built a stochastic model connecting the macroscopic film resistance to a generic class of microscopic mechanisms, characterized by local variations of resistivity due to structural rearrangements. Two specific realizations of the model, with tailored choice of the parameters, were shown to reproduce, qualitatively as well as quantitatively, the characteristic jumps of resistivity, as well as other, less obvious, features of the experimental data. As local structural rearrangements happening in nanogranular films can be investigated by means of atomistic simulations that only have to involve a few nanoparticles, such model can be considered as a starting point for further theoretical investigations aimed at shading light upon such mechanisms.
In summary, we presented a modeling based on EMA that can be used to study the the resistive switching phenomena that seems to characterize nanogranular film also well beyond the percolation regime and we showed how simple thermally regulated local structural changes can explain the most relevant features of the experiments, suggesting a clear path for future investigations.

ACKNOWLEDGMENTS
This work was fully funded by Fondazione CON IL SUD (Grant No: 2018-PDR-01004).