A coarse-grained Stochastic Resistor Network model
Aiming to describe the electrical behavior of macroscopic cluster-assembled metallic films at a coarse-grained level, our Stochastic Resistor Network model is essentially based on a large three-dimensional (3D) regular lattice of resistors, each one capable of a discrete number of conductive states, which in our work has been fixed to four (see Fig. 2 and Methods). This choice is a trade-off between descriptive capability and the need to limit the complexity of our simulation framework. Resistors are represented as links ij joining a pair of nodes (i, j), which are coarse-grained representations of large sub-regions of the original sample; the links are organized in groups termed layers (which are three in Fig. 2), connected to each other by vertical (z-axis) links. The conductance values were established based upon experimental data on nanostructured Au films. As shown in Fig. 2 our resistor network is provided with an input and an output node, to which an overall voltage \(\Delta V\) is applied.
Left: schematization of a node, a link and viem of the layer 1 of the resistor network. Each link can assume one of a discrete number of conductance values (namely, in this work, \(\sigma _\alpha =10^-11 \frac1\Omega \), \(\sigma _\beta =10^-3 \frac1\Omega \), \(\sigma _\gamma =2 \times 10^-3 \frac1\Omega \) and \(\sigma _\delta =4 \times 10^-3 \frac1\Omega \)). In the picture, links with \(\sigma _ij=\sigma _\alpha , \sigma _\beta , \sigma _\gamma , \sigma _\delta \) are respectively white, orange, blue and purple. Right: 3D view of the full network, where only edges with \(\sigma _ij=\sigma _\gamma \) are colored. \(N_x, N_y, N_z\) correspond to the number of nodes along each network direction. Picture realized with NetworkX44.
For each link connecting nodes i and j, the applied voltage induces the presence of electrical current \(I_ij\) flowing through it, and also a potential difference \(\Delta V_ij\). Applying the spectral theory to the Laplacian matrix of the weighted undirected graph associated to our network31,32,45,46,47,48, we are able to retrieve \(I_ij\) and \(\Delta V_ij\) for each link of the network. The set of these \(I_ij\) and \(\Delta V_ij\) values constitutes the input for Monte Carlo update moves which make the system stochastically evolve. These MC moves are conceived to reproduce the thermal stability of the connections and also other non-linear electron conduction mechanisms, due to inter-cluster and intra-cluster atomic rearrangements, which result in the dynamical creation and destruction of conduction pathways and trigger the switching events24,33,34,35. Such physical effects are reflected in our model through the possibility for each link ij to either increase its own conductance or decrease it, via stochastic jumps across the available discrete conductance levels, due to the heat released by its neighbors (\(\sigma _ij\) grows with a given probability) or to its thermal dissipation (\(\sigma _ij\) lowers with a given probability), respectively—see Methods. Besides, \(\Delta V_ij\) is nonlinearly used either to stochastically degrade or improve \(\sigma _ij\), by comparing it with a suitable threshold \(\Delta V^\mathrmth\) (see Methods for details). The computational cost of the model is highly demanding, since the application of the aforementioned stochastic update moves produces a new configuration of the network, which then requires a new complete solution in terms of currents \(I_ij\) and voltages \(\Delta V_ij\) at each link, to provide the input for the subsequent MC step. In view of the computational effort required, our aim has been to find out and focus on a well-founded set of parameters capable of reproducing the experimental phenomenology.
Qualitatively reproducing experimental electrical transport properties via the SRN model
Our SRN model is capable of generating a rich and complex electrical transport phenomenology which is impossible to forecast in advance; nonetheless, the key ingredients were conceived in the attempt to qualitatively retrieve the experimental features of Fig. 1. Remarkably, we see the progressive emergence of all the peculiarities observed in the electrical conduction properties of nanostructured gold films, as long as the network size is gradually increased. In particular, the size of our regular 3D lattice has been progressively enlarged, up to the size of \(N_x=27 \times N_y=42 \times N_z=3\), which corresponds to 3404 nodes and 8919 links. The simulation of such a large network, endowed with a nontrivial stochastic dynamics, represents an unprecedented attempt to study the complexity required to describe electrical conduction phenomena in cluster–assembled nanostructured metallic systems. Our results thus indicate the first identification of a minimum complexity limit to be considered in order to start achieving such an experimental phenomenology. The simulations of our SRN model have required an extensive use of high performance computating facilities. Despite this size and the resulting intricacy, our SRN model cannot capture the full complexity of the experimental system; for this reason we can only expect a qualitative reproduction of the phenomenology observed in gold nanostructured films. We also point out that, keeping the potential difference \(\Delta V\) fixed, the effect of rescaling by a constant factor, \(\Gamma\), all the conductances, \(\sigma _ij\), and correspondingly, adequately rescaling all the model activation thresholds, has as its only effect an equivalent stochastic dynamics of the model, characterised by the same \(\Delta V_ij\) with currents \(I^\prime _ij = \Gamma \times I_ij\). We show in Fig. 3 the analogue quantities displayed in Fig.1, simulated via the SRN model.
(a) typical evolution of simulated electrical resistance under the application of \(\Delta V =1\ \hbox V\) (red) or 15 V (lavender); measurements are taken every MC step. Black and dark blue dots correspond to resistance switches. (b) Inter-Switch-Interval probability density, computed as for the experimental data, joining data from 10 statistically independent simulations. Also here, at 15 V the distribution is approximately monotonically decreasing and the majority of the RS events are close in time (within \(\approx 1000\) MC steps). The data taken at 1 V feature a longer tail of the distribution, as in the experiments at low voltage. (c) Simulated PSD for measurements taken at 1 V (blue-blue) and at 15 V (dark purple), as a function of frequency f, in log-log scale. Yellow and blue dashed lines represent \(1/f^\alpha \) fitted curves. d: simulated \(I(\Delta V)\) curve, where each point is averaged over 20 independent simulations lasting 8000 MC steps each (the dashed line is only a guide for the eyes).
Figure 3a shows the evolution of \(R_\mathrmtot\) during a typical simulation of 20000 MC steps (where a MC step is the simulation basic time unit) at \(\Delta V = 15 V\) for our SRN model. The analysis of the resistance series data for two-electrode devices (see Fig. 7, Methods section) is carried out following the same protocol used for the experimental data. The qualitative similarity between the experimental and simulated patterns suggest that the SRN model is able to capture the main features evidenced in the experiment: a number of distinguished resistance levels are repeatedly visited by the system with sudden variations of \(R_\mathrmtot\) and subsequent fluctuations around a given resistance value. In the stochastic evolution of \(R_\mathrmtot\) we observe more intense fluctuations with respect to the experimental case; the reason for this lies, as discussed below, in the lower complexity of the SRN model compared to the real system and the consequent greater fragility and susceptibility of \(R_\mathrmtot\) to the stochastic dynamics.
We evaluated the Inter-Switch-Interval distribution of \(R_\mathrmtot\) and the the Power Spectral Density of \(R_\mathrmtot\), i.e. the square modulus of the Fourier transform of the signal. The Power Spectral Density of the measured electrical signal for our SRN model is characterized by a \(1/f^\alpha \) trend24, with typical \(\alpha\) values between 1.5 and 2. ISI distribution is obtained like in the experimental data analysis. The data are plotted in form of probability density distribution function in Fig. 3b for devices polarized at 15 V. Notably, the ISI distributions coming from the SRN model are characterized by the same trend obtained from the experimental data. Even for the SRN model, at short times the data decrease like an exponential with similar timescales for both voltages, while extended tails of outliers are present, as also found in the experiments (see Fig. 1b). The majority of the RS events are most frequently separated by less than 70 MC steps, with very long tails of the distributions for larger MC steps intervals. The comparison of the experimental and simulated ISI distributions allows us to roughly associate to a single MC step an equivalent time, for the observed phenomenology, of the order of 0.1 seconds.
In Fig. 3c, a typical PSD of the SRN total resistance is reported in log-log scale, for simulations performed at 1 V and at 15 V (lavender and dark read points, respectively). Dashed dark blue and pale red lines represent power-law fits of the simulated data, with exponent \(\alpha\). Surprisingly, the simulated PSD associated to \(R_\mathrmtot\) shows a \(1/f^\alpha \) behavior with \(1<\alpha <2\) (in the particular example shown, \(\alpha =1.74\) and \(\alpha =1.72\) at 1 V and and at 15 V respectively). This behavior evidences a pink-noise memory25, qualitatively similar to the noise generated by the cluster-assembled films studied in the experiments24,35.
Our SRN model presents a peculiar \(I(\Delta V)\) relation, featuring very small current values at small voltages, a markedly different regime with a steep increase in I at intermediate voltages, and a saturation regime at higher voltages (see Fig. 3d). The regions at low and intermediate \(\Delta V\) qualitatively resemble the experimental behavior displayed in Fig. 1d. Despite its considerable number of links and nodes, the SRN model unveils finite size effects at high voltages: the limited size of the network, jointly with its collective and cooperative dynamics, yields a saturation in \(|I(\Delta V)|\) values. Given the remarkable qualitative reproduction of the experimental phenomenology, our simulations could provide insights about the evolution and modification of the structure of cluster-assembled Au films induced by current flow and responsible for the observed switching behavior. We concentrate on the effect of high and low voltage application on the modification of the electrical behavior of different junctions24,35.
In Fig. 4, we report the evolution of the fractions of links having \(\sigma _ij=\sigma _\alpha , \sigma _\beta , \sigma _\gamma , \sigma _\delta \) as a function of the MC steps, during a specific simulation where the effect of the application of a high voltage bias (\(\Delta V=\)15 V) to our resistor network has been studied. Initially, a low voltage is applied to a network featuring a purely random conductance distribution. First of all, we notice that, almost independently of the actual initial spatial distribution of the conductances, the system reaches a dynamical equilibrium as a consequence of the stochastic evolution of the model (see Fig. 4, slightly before 20000 MC steps). The application of a high voltage bias significantly alters this steady state, leading the system to a markedly different condition, characterized by a further decrease in the amount of the high conductances and a concurrent increase of the fraction of links having \(\sigma _ij=\sigma _\alpha \).
Evolution of the fractions of links with: \(\sigma _ij=\sigma _\alpha =10^-11 \Omega ^-1\) (grey), \(\sigma _ij=\sigma _\beta =10^-3 \Omega ^-1\) (orange), \(\sigma _ij=\sigma _\gamma =2 \times 10^-3 \Omega ^-1\) (navy) and \(\sigma _ij=\sigma _\delta =4 \times 10^-3 \Omega ^-1\) (purple) as a function of the MC steps. Black dashed vertical lines separate the three simulation phases. At the start of the simulation, one can see the initial probability for the abundance of links having conductivity equal to \(\sigma _\alpha , \sigma _\beta , \sigma _\gamma , \sigma _\delta \). The initial amount of each conductance level is chosen in a range of values which allow for a sufficient degree of current percolation.
In the subsequent phase, where \(\Delta V= 1\hbox V\) again, the stochastic evolution leads it to the onset of a novel dynamical equilibrium (nonetheless characterized by the persistence of the resistive switching activity, see Fig. 3a), in close analogy with what is observed for the physical substrate: the number of highly resistive links almost doubles, corresponding to a picture where a significant reduction of the available paths for the current occurs. This is confirmed by the analysis of the shortest paths (measured weighting each link with the inverse of the current it is traversed by) available for the current from the input to the output, whose number is strongly reduced after that the system experiences a high voltage bias (see Supplementary Information). This result shows that the SRN model, whose stochastic dynamics is strongly correlated, is highly sensitive to the electrical conditioning history of the system, and supports the model based on local rearrangements of grain boundaries to explain the non-linear and non-local conduction properties of cluster-assembled Au films24,35.
Information-theoretic analysis of correlations in the experimental device and the SRN model
In neuroscience, the role of dynamical correlations between different brain areas is recognized as fundamental for cognitive and behavioral integration21,39. Statistical measures derived from information theory have been proposed to characterize the integration of information among functionally segregated groups of neurons, in particular entropy, Mutual Information (MI) and Integrated Information (IN) have been considered to reveal the degree of interconnection/segregation for different regions of the brain or biological neural networks in response to external stimuli39.
Aiming to get a deeper insight into the presence and role of spatial correlations in cluster-assembled Au films, we performed an experimental analysis exploiting the same MI and IN tools21,39. To investigate the information content of the measured and simulated resistance series data, Mutual Information and Integrated Information have been computed as detailed below and in the Methods section. The nanostructured cluster films can indeed be idealized like networks of smaller units with a proper activity with a stochastic behavior. In literature different approaches are used to define the state of such a kind of system22,49. In a general approach, we can consider that the intensity of the electrical activity of each unit corresponds to different states of the elemental units.
Let us consider a generic set X of N Random Variables (RVs), whose subsets of size \(k \le N\) are indicated as \(X_j^k\), where j runs over all the \(\dfracN!k!(N-k)!\) possible choices of a subset of size k. Therefore, the elementary units of X can be referred to as \(X_j^1, j=1 \dots N\). In the following equation, we give the definition of entropy \(H(X_j^k)\) of a generic subset \(X_j^k\), possibly even of \(X=X_1^N\):
$$\beginaligned H (X_j^k) = \sum _i p_i (X_j^k) log_2(p_i (X_j^k)) \endaligned$$
(1)
being \(p_i (X_j^k)\) the probability to find the subset \(X_j^k\) in its i-th state. Exploiting the above general definition of entropy, two useful quantities can be introduced: the Mutual Information and the Integrated Information. MI takes into account the relationship between a subset \(X_j^k\) and the complementary subset \(X-X_j^k\), whereas IN includes the correlations among the basic units which are part of \(X_j^k\)21,39:
$$\beginaligned&\text MI(X_j^k,X-X_j^k)=\text MI(X-X_j^k,X_j^k)=H(X_j^k)+H(X-X_j^k)-H(X) \endaligned$$
(2)
$$\beginaligned&\text IN(X_j^k)=\sum _l \, : X_l^1 \in X_j^k H(X_l^1)-H(X_j^k) \endaligned$$
(3)
In Eq. (3), the index l runs over all the elements \(X_l^1\) belonging to the subset \(X_j^k\). MI describes the statistical dependence between the RVs represented by the subset \(X_j^k\) of X and the complementary \(X-X_j^k\). Note that MI\((X_j^k,X-X_j^k)=0\) if the two partitions of X are composed by independent RVs, while positive values of MI mark the presence of general statistical correlations between two data sets. Integration is, on the other hand, a tentative to give a measure of the segregation (independence) or integration (dependence) among the elementary constituents of the system. IN \((X_j^k)=0\) if the RVs composing the chosen subsystem, \(X_j^k\), are statistically independent one from each other, while IN\((X_j^k)>0\) in presence of correlations among them.
Note that, in principle, one can compute also the average value of MI and IN:
$$\beginaligned&\overline\text MI_k :=\dfrack!(N-k)!N!\, \sum _j \text MI(X_j^k,X-X_j^k) \endaligned$$
(4)
$$\beginaligned&\overline\text IN_k :=\dfrack!(N-k)!N!\, \sum _j \text IN(X_j^k) \endaligned$$
(5)
In the case of our multi-electrode cluster-assembled films37, the system X can be thought of as made by elementary units \(X_j^1\), which are single electrode couples (see Methods, and especially Fig. 7b for a system with \(3 \times 3=9\) possible electrode pairs), where an electrode couple includes one of the three input terminals (1,2,3) together with one of the three output ones (A, B, C). The probability that the two electrodes selected have a given electrical resistance represents instead the basic ingredient from which we aim to infer entropy-related properties. Computing entropy for a set of k electrode couples, \(H(X_j^k)\), using Eq. (1) requires having at disposal an approximated probability distribution \(\p_i\\), where a probability is associated to a state i, corresponding to a list of k resistances simultaneously taken by all the involved electrode couples. More in detail, the time series of the resistance between each pair of electrodes is measured for a set of times \(\t\\); each resistance time series \(R_\t\ (X_j^1)\) is normalized dividing by its mean value \(\barR(X_j^1)\) and then discretized (see Methods). The entropy related to a set \(X_j^k\) of electrodes (and MI and IN as a consequence) can be then obtained computing the joint probability of the resistance states of the k electrode couples belonging to the set.
Inspired by Refs.21,39, in order to understand how correlation between different regions of cluster-assembled films rises in response to the applied external voltage pulses, we experimentally measured MI and IN using a device shown in Fig. 7 and starting from the resistance time series of all the different electrode couples. In the following, we focus on data collected at \(\Delta V=\)1 V, after the conditioning stage (i.e. the application of a high voltage bias of 15 V). As an example, in Fig. 5a, \(\overline\text MI\) is plotted as function of the number k of electrode couples, while applying the voltage bias either to the electrode couples 1-A (azure curves) or to 1-C (orange curves); Fig. 5b shows the equivalent result for \(\overline\text IN\). \(\overline\text MI\) and \(\overline\text IN\) both show an increasing trend as a function of the subset size k. In particular, the Mutual Information monotonically increases, reflecting the correlations of a subset with its complementary subset. The application of \(\Delta V\) to the electrode couple 1-C produces higher MI values: intuitively, in this configuration, the current flows throughout a larger region of the sample, probably enhancing reciprocal correlations among the sub-regions. Simultaneously, \(\overline\text IN\) is larger when the voltage is applied to 1-C, mirroring the growing internal correlations among the elements of the chosen subset. In both cases, the growth of \(\overline\text IN\) is sublinear, whereas a linear growth (represented as a grey dashed line in the picture) would indicate a fully correlated system. We observe a remarkable similarity with the trends usually observed in neuroscience works, where measurements are taken on animal brains21,39. Note that, in that field, the area comprised between the linear growth and the effective \(\overline\text IN\) curve is a measure of the so called neural complexity.
(a) experimental \(\overline\text MI(k)\) for two different electrode couples to which \(\Delta V=1\) V is applied (1-A: azure and 1-C: orange). Data are collected over 275 distinct configurations. (b) \(\overline\text IN(k)\) for the same experimental data. Data are collected over 275 distinct configurations. In both cases, errorbars are given by the evaluation of these observables over all the different ways of choosing a given k. \(\overline\text MI\) is only measured up to \(k=5\), since \(\text MI(X_j^k,X-X_j^k)=\text MI(X-X_j^k,X_j^k)\), while \(\overline\text IN\) is measured for any \(k \in [0,9]\). Grey dashed lines represent the \(\overline\text IN\) growth in a fully correlated system.
(a) \(\overline\text MI(k)\) computed from two statistically independent simulations at \(\Delta V=1\) V (blue and red), as a function of the number of subregions, k, considered. (b) \(\overline\text IN(k)\) computed for the same simulations. In both cases, data are collected over 200 distinct MC configurations and errorbars are given by the evaluation of these observables over all the different ways of choosing a given k. Grey dashed lines represent the \(\overline\text IN\) growth in a fully correlated system.
In the complementary analysis of the simulations, the investigation of spatial correlations plays a crucial role as well. The topology of the nanostructured film just analyzed (see Fig. 7b) is more structured than the simulated network and it has more electrodes; therefore, the Information-theoretic tools based analysis in the SRN model necessarily requires a modification of the notion of \(X_j^k\), which we choose to correspond to a well-defined sub-region of the network (see Fig. S4). In fact the network can be easily divided into N sub-regions (here we set \(N=7\)), each of them gathering the electrical properties of a number neighboring of links. The resistances of each sub-region (averaged over the links belonging to it) are periodically recorded in time, and their distribution is discretized with the same procedure employed for the experimental data (see Methods). In this way, it becomes possible to evaluate \(\overline\text MI(k)\) and \(\overline\text IN(k)\) even in the simulated system. In Fig. 6 we show \(\overline\text MI(k)\) (panel a) and \(\overline\text IN(k)\) (panel b) for two statistically independent simulations at \(\Delta V=1\) V, after the application of a bias of 15 V. Statistical independence is guaranteed by the different initial configuration of conductances of the two SRNs used in simulations 1 and 2, which triggers a different electrical evolution for each of them. Remarkably, also \(\overline\text MI(k)\) and \(\overline\text IN(k)\) extracted from the SRN model simulations present some remarkable correlations among the various sub-regions of the whole simulated network. The growth trend of \(\overline\text MI(k)\) and \(\overline\text IN(k)\) as a function of the subsystem size k is exemplified in Fig. 6, and it is a general feature observed in all our simulations after the application of high voltage bias (see Supplementary Information for further discussion). Conversely, without the previous application of a high voltage, in experiments and as well in simulations we observe a substantial lack of dependence of \(\overline\text MI(k)\) and \(\overline\text IN(k)\) on k. Notably, the voltage which activates such behavior corresponds to the one which triggers the switching events, and, in the SRN model, promotes the reduction of the number of available shortest paths. Therefore, it turns out that a high \(\Delta V\) bias is crucial to allow the system to visit neatly distinguished resistance states in the subsequent phase with a small applied tension. All these results suggest that the experimental and simulated systems have a complex and collective response to external stimuli, with an emerging electrical behavior determined by a subtle balance between the random dynamics and the reciprocal influence among different regions.
