Percolation with plasticity for neuromorphic systems

We develop a theory of percolation with plasticity systems (PWPs) rendering properties of interest for neuromorphic computing. Unlike the standard percolation between two large electrodes, they have multiple ($N\gg 1$) interfaces and exponentially large number ($N!$) of conductive pathways between them. These pathways consist of non-ohmic random resistors that can undergo bias induced nonvolatile modifications (plasticity). The neuromorphic properties of PWPs include: multi-valued memory, high dimensionality and nonlinearity capable of transforming input data into spatiotemporal patterns, tunably fading memory ensuring outputs that depend more on recent inputs, and no need for massive interconnects. A few conceptual examples of functionality here are random number generation, matrix-vector multiplication, and associative memory. Understanding PWP topology, statistics, and operations opens a field of its own calling upon further theoretical and experimental insights.


Introduction
Devices for neuromorphic computing remain among the most active areas of research with a variety of models for neurons, synapses and their networks [1][2][3][4][5][6][7]. They are typically built of nonvolatile memory cells and interconnects wired in a certain architecture.
Here we introduce a concept of neuromorphic architectures where neither artificial memory cells nor interconnects are required. They are based on natural disordered materials with percolation conduction [8][9][10] such as amorphous, polycrystalline, and doped semiconductors, or granular compounds. Among percolation conduction materials, we consider those exhibiting plasticity, i.e. the ability to undergo nonvolatile changes in their resistances in response to strong enough electric fields. They include metal oxides and chalcogenide compounds used with resistive random access memory (RRAM) [11] and phase change memory (PCM) [4], granular metals [12] and nano-composites [13].
A commonly recognized mechanism of the above mentioned plasticity assumes a narrow dielectric gap that can be bridged by a conductive filament. The filament is created in response to electric bias of certain polarity and is dissolved by changing the polarity to opposite (bipolar RRAM); alternatively, it can be created by a pulse of certain amplitude and length and dissolved by a pulse with different parameters (unipolar RRAM and PCM). The details of such microscopic mechanisms have been discussed in multiple publications on PCM and RRAM devices; they are of no significance here.
The nonvolatile resistance changes play the role of memory records. In particular, conductive or insulating states in PCM or RRAM device can represent the logic states 0 and 1. The continuing quest for multivalued memory is aimed at developing structures with multiple nonvolatile memory states. Our consideration below renders a unique approach to multivalued memory based on a variety of natural percolation with plasticity (PWP) materials.
A prototype PWP illustration in figure 1 shows several percolation conductive pathways connecting multiple electrodes. It assumes certain voltages E i applied to all electrodes but one used to measure the electric current I. We will see in section 4 below how a measured value of I can be interpreted as a product of vector of voltages and matrix of conductances, thus presenting an example of such a computing operation as the

Standard percolation conduction
We start with a recap of the pertinent percolation concepts [8][9][10] describing electric transport through a network of exponentially different random resistors. It is dominated by the percolation cluster formed by conductive bonds connecting the electrodes as illustrated in figure 2. The bonds are composed of the elements with minimum resistances whose total concentration is just sufficient to form the electrode connecting pathways. The characteristic mesh size of the percolation cluster L c is much greater than the linear size a of a microscopic resistor (see figure 2). It is effectively uniform over large distances L L c where L c is the correlation radius determining its characteristic mesh size. A sketch of the percolation cluster between two electrodes. Shown under magnification is a bond domain consisting of microscopic resistances in series. The nature of resistances is not specified. Most work on percolation assumed resistances of electronic nature. However, they can equally represent the mechanism of ionic conduction [28]. Each bond of the cluster consists of a large number of microscopic resistors R i . Their exponential randomness is described as R i = R 0 exp(ξ i ) where quantities ξ i are uniformly distributed in the interval (0, ξ max ) and i = 1, 2, . . ..
The physical meaning of ξ depends on the type of a system. For definiteness, we assume here ξ i = V i /kT corresponding to random barriers V i in noncrystalline materials where k is the Boltzmann's constant and T is the temperature. In reality, the nature of percolation conduction can be more complex including e.g. finite size effects and thermally assisted tunneling between the microscopic sites in nanocomposites [30,31]. These complications will not qualitatively change our consideration below.
The cluster constituting microscopic resistors exhibit non-ohmicity due to the field induced suppression of their barriers V i , according to where J i is the resistor current, J 0 is a constant, U i = E i a is the voltage applied to the barrier, a and E i are, respectively, the barrier width and local electric field, and q is the electron charge. An important conceptual point [16,[32][33][34] is that the applied voltage concentrates first on the strongest resistor of a percolation cluster bond (resistor 1 in figure 3) suppressing it to the level of the next strongest (resistor 2 in figure 3), so the two equally dominate the entire bond voltage drop. It then suppresses the nextnext strongest resistors (3, 4, 5, . . . in figure 3), etc. Based on the analysis of such a percolation process, it was derived that the entire percolation cluster changes its structure as L c = a V 0 /qEa, and ΔV = V 0 qaE (2) resulting in the macroscopic non-ohmic conductivity, Here E is the macroscopic electric field strength, L c and ΔV are, respectively, the field dependent correlation radius and maximum barrier decrease in the percolation cluster, and V 0 is the amplitude of barrier variations. Two assumptions behind the above discussed nonohmicity are: (a) the volatility of bias induced changes where each microscopic resistor adiabatically adjusts its resistance to the instantaneous bias. (b) The quasistatic nature of biasing implying time intervals exceeding the relaxation times of all microscopic resistors, i.e. t τ m = τ 0 exp(ξ max ) where τ 0 ∼ 1 ps is on the order of the characteristic atomic vibration times. The quasistatic bias limitation can be relaxed to describe non-ohmicity in the pulse regime [17]. It was argued that for a cluster bond a bias pulse U of length t generates the current This result applies when t is shorter than the maximum relaxation time τ max = τ 0 exp(ξ max ) and is formally different from that of dc analysis [32] by the substitution [17] ξ max → ξ t ≡ ln(t/τ 0 ); the two results coincide when ξ t = ξ max . For the entire percolation cluster, the modification ξ max → ξ t predicts the current, to the accuracy of an insignificant numerical coefficient in the exponent [16,17].
To avoid any misunderstanding, we note that the field induced change in resistances of percolation clusters described in this section is volatile (i.e. it disappears when the field is removed). It should not be confused with the nonvolatile plasticity introduced in the next section.

Percolation with plasticity systems
PWP phenomenon differs from the standard percolation conduction in both topology and non-ohmicity. The former is such that the proposed PWP has multiple (N 1) electrodes; the latter is due to the nonvolatile nature of bias induced changes. These properties can be somewhat different for PWP systems with large (L L c ) and small (L L c ) geometrical dimensions L as described next.

General
We start with noting some general properties of PWPs relevant to neuromorphic applications.
(1) The variations ΔR ij (with respect to the averages R ij ) in the interelectrode resistances R ij are random quantities that are uncorrelated with any desired accuracy for not-too-close electrodes. Consider exp(ξ α ) for a bond of N ij resistors and ξ α uniformly distributed in the interval (ξ min , ξ max ). It is then straightforward to obtain the correlation coefficient between the resistances of (i, j) and (k, l) bonds, (6) where N ij and N kl represent the numbers of microscopic resistances in those bonds, and N s is the number of resistances shared between them. We describe each bond as a random walk. Then, if the two bonds are not close geometrically, separated by distances L ijkl exceeding a N ij + N kl , then their overlap is exponentially small, and The averages implied by the definition for ΔR ij can be readily measured for an ensemble of geometrically similar pairs, such as (1,9), (2,8), (3,7), (4,12), etc. in figure 1.
(2) Related to the above item (1), there is a characteristic distance l ij , above which the electrodes i and j are electrically independent (no crosstalk between them). To estimate that length, we use it in place of L ijkl in equation (7) setting also N ij ∼ N kl ∼ L/a where L is the sample size and a is the microscopic resistor length. This yields, Assuming as a rough estimate a ∼ 10 nm and L ∼ 1 cm, yields l ij ∼ 10 μm.
(3) High frequency inductive coupling of conductive pathways can be estimated based on a model of a wire and a loop of diameter L. Using the standard electrodynamics, the ratio of the induced current over the primary current then becomes (in the Gaussian system) Here ω is the frequency, c is the speed of light, and R is the resistance of the bond. Assuming as a reference, values L ∼ 1 cm, R ∼ 1 Ω [= 1/(9 × 10 11 ) s cm −1 in Gaussian units], and ω ∼ 1 GHz, the latter ratio is of the order of 1. However, lower frequencies and especially higher resistances (typically above 10 3 Ω with PCM and RRAM applications) make that ratio small and acceptable.
(4) The characteristic RC times related to the writing and reading processes for a single microscopic resistor can be estimated as ∼ 10 ps for R ∼ 1 Ω and C ∼ 10 pF (corresponding to a 1 cm sample with the same macroscopic resistivity as that of 1 MOhm RRAM resistor with 10 nm linear dimension), rather competitive against the background of modern technology.
(5) The property of plasticity takes place when the local electric field exceeds its material dependent threshold value for resistance switching as discussed next.

Large PWPs
The concept of infinite percolation cluster survives if the electrode sizes l and interelectrode distances L ij are much larger than L c , in which case (we call it 'large PWP') resistances R ij are determined, in the ohmic regime, by electrode geometry, similar to the case of steady currents between finite size electrodes in massive conductors, such as grounding electrodes in a soil. Here f ij is a dimensionless function whose shape depends on the electrode locations through the confinement of electric currents by the sample boundaries. The macroscopic conductivity σ in the equation for R ij is taken in the limit of infinitely large percolation systems where it is uniquely defined. While f ij can be numerically modeled with for particular electrode configurations, some general statements can be made based on the available examples and a simplified model [35]. The latter presents an electrode as a metal hemisphere immersed into a macroscopically uniform medium formed by the percolation cluster over large scales as illustrated in figure 4. A rough estimate is given by to the accuracy of a numerical multiplier that depends on a particular electrode geometry. Here O(x) means 'of the order of x' and it is assumed that l/L 1, l/L ij 1. It follows from equations (10) and (11) that resistances R ij do not depend strongly on the interelectrode distances in the ohmic regime. Yet, the differences in the interelectrode resistances will exist due to statistical fluctuations in their connecting bonds. The bond of length L ij in a large PWP will contain L ij /L c quasiindependent cells of the percolation cluster. Each cell has on average resistance R c making the average bond resistance R ij = R c L ij /L c . Taking into account that resistances of individual cells exhibit random variation of the order of R c , one thus arrives at the estimate for the characteristic relative variations of resistances, While relatively small, the fluctuations δR/ R ij are still significant enough to experimentally discriminate between different interelectrode resistances. For example, δR/ R ij ∼ 0.1 assuming L c ∼ 10 nm and L ij ∼ 1 μm. Unlike the standard percolation between two flat electrodes, in large PWP the electric field systematically decays with distance r l from a small electrode due to the current spreading (again, similar to the case of grounding electrodes) as illustrated in figure 4. That geometrical effect will significantly alter the nature of nonohmicity making it most significant in the proximity of r ∼ l around the electrode. One can show that due to that field suppression, the steady state non-ohmic conduction will be limited to r ∼ l ξ max qaU/LkT (13) where U is the total voltage applied to a sample of length L. The pulse regime limitation will be described by a similar formula with the substitution ξ max → ξ t . Figure 4. A large PWP geometry with two hemisphere electrodes and a fragment of percolation cluster resolved in the magnifying glass. Shown in dash are the equipotential (arcs) and electric field (arrows) lines in the proximity of the right electrode illustrating the geometrical effect of electric field spreading with distance from the electrode.

Small PWPs
The concept of infinite cluster fails when l and/or L ij are smaller than L c ('small PWP'), in which case one has to consider multiple conductive paths unrelated to the infinite cluster. Both the cases of large and small PWP are possible across a broad variety of percolation systems. For example, the current L ∼ 10 nm-node technology belongs in small PWP with [16] a 0.3 nm and V 0 /T 100, i.e. L c ∼ 30 nm. With the latter parameter values, increasing L to microns and beyond will result in large PWP networks. The resistances of conductive pathways in small PWP exhibit significant variability. Next, we estimate their statistics. Based on equation (4) a chain resistance is estimated as R = R max exp(−δV/kT) where δV = V max − v max and v max is the maximum barrier in that chain, R max = R 0 exp(V max /kT), and V max is the maximum barrier in the entire system, V max v max . Assuming uniformly distributed barriers, the average number of resistors with the barriers above a given v max in an n-resistor chain is where V min is the minimum barrier. The Poisson probability of a chain having no barriers greater than v max is P V (n) = exp(−n v ), and the probability of finding a chain with a maximum barrier in the interval kT around v max is Multiplying P(v max ) by the probability P n (L) ∼ (L/an) exp(−L 2 /na 2 ) of an n-resistor chain connecting points distance L from each other, we obtain the probability density of n-chain with a given barrier v max (to the accuracy of a numerical multiplier in the exponent). Integrating that product over n by steepest descent and expressing v max through R yields the probabilistic distribution density, It follows that resistance spectrum is a gradual function with a certain characteristic width ΔR. As estimated separately for the cases of small and large L, width ΔR can be approximated for the entire range of L by the following equation: For all practical values, it encompasses multiple orders of magnitude. Because of the dispersion in the values v max between different pathways, the nonohmicity exponents will vary from one R ij to another. More specifically, instead of ξ max in equation (4) the value v max /kT should be used with the probability distribution of equation (14). That additionally broadens the distribution of path resistances in the nonohmic regime; we omit here the obvious formal description of that effect.

Plasticity by switching
A unique non-ohmicity feature of PWP is its nonvolatile nature rendered by the underlying material (say, of PCM or RRAM type). Each microscopic element of a conductive path can exist in either high or low-resistive state whose respective resistances, R > and R < , are markedly different. R > resistances are random, all exceeding R < . The applied bias concentrated on the strongest of R > resistors (in the manner of figure 3) will change them to R < by switching, i.e. by long-lived structural transformation not adaptable to subsequent voltage variations. In the steady state bias regime, the next strongest resistor will be stressed with practically the same voltage as opposed to the above discussed case of volatile non-ohmicity in the standard percolation clusters. In the first approximation, an originally resistive percolation bond will transform into its conductive state by n discrete steps where n is the number of its microscopic resistances, similar to a falling row of dominoes arranged in the order of descending ξ's.
The latter behavior can be more complex in large PWPs due to the geometrical field distortion illustrated in figure 4. There, the field will eliminate large resistances in the region of characteristic length given by equation (13) which then becomes a sort of low resistive protrusion into the bulk material. Such a protrusion will concentrate the electric field similar to the lightning rod effect. As a result, switching will take place in the next high field domain growing that protrusion further, etc, until it reaches the opposite electrode. While the kinetics of such a process can be readily described, we will omit it here.
For switching to occur, the local field on a microscopic resistor must exceed a certain critical value [36,37] E c , typically on the order of 10 5 − 10 6 V cm −1 , and E c decreases logarithmically with the electric pulse (spike) length [37][38][39][40][41]. That temporal dependence opens a venue to the spike timing dependent plasticity (STDP) [29], which is another important property of neural networks. Because almost the entire voltage drops across a microscopic resistance of small linear size a, the microscopic field E is significantly stronger than the apparent macroscopic field and was shown [17] to approach ∼20 MV cm −1 , well above the values of E c sufficient for threshold switching.

Plasticity in the pulse regime
Pulse excitation regime brings in additional physics [17,42]. The random elements of a percolation bond will behave as resistors when their relaxation times τ are shorter than the pulse duration t; however the elements with τ t will act as capacitors causing no significant voltage drop in response to short pulses. The boundary between these fast and slow elements is defined by τ = t, i.e. ξ = ξ t ≡ ln(t/τ 0 ) given that τ = τ 0 exp(ξ).
The resistance R t = R 0 exp(ξ t ) will be the highest of all bond elements resistances and the voltage pulse will concentrate on it for its duration. Should the corresponding local field exceed E c , the resistance R t will switch to R < with other resistors intact, which leads to the interelectrode path resistance decreasing by roughly a factor of [e] (base of natural logarithms) as illustrated in figure 5. It is straightforward to show that the number of such stepwise changes in a bond resistance is estimated as which property can, in principle, be used to create multi-level memory operated by trains of pulses.
Corresponding to equation (17), the characteristic voltage per microscopic resistor is given by [17], where U L is the macroscopic voltage across the bond of length L. According to the numerical estimates in section 7 below, U 1 can create the microscopic field of strength exceeding the switching value. The concept of slow resistors acting in a manner of capacitors has been proposed and verified earlier [42][43][44]. It may be appropriate to additionally explain here that capacitors do not accommodate significant voltages when in series with resistors because they conduct the displacement currents, j D = (ε/4π)(dE/dt), where ε is the dielectric permittivity. The overlay between the capacitor and resistor regimes takes place when j D = j = σE where j is the real (charge transport) current, E is the electric field strength, and σ is the conductivity. The displacement current through a capacitor is due to the rate of field change, unrelated to voltage, rather than the field itself proportional to voltage in a resistor.
The ratio of displacement vs real current can be presented as j D /j = (ε/4πσ)(d ln E/dt). The expression in the first parenthesis represents the Maxwell's relaxation time τ , while the reciprocal of the second parenthesis gives the pulse duration. That takes us again to the criterion τ t for the element of a percolation cluster operating in the capacitive mode.
Relating this understanding with microscopic models, we note that the displacement currents are due to charging/discharging processes in, say, capacitor electrodes, or in certain defect configurations responsible for electric potential distributions in percolation clusters.
To avoid any misunderstanding, we note that only electric currents caused by external stimuli are considered here. In particular, possible purely diffusion currents due to concentration gradients remain beyond our consideration. Such diffusion currents can play a significant role in system transformations including both switching mechanisms and degradation; both remaining beyond the scope of this work.

Reverse plasticity
The above description does not address the important question of reversibility of switched structures back to high resistive states. The feasibility of such reverse process (commonly referred to as RESET) follows from the known practices of RRAM and PCM operations.
One conceivable answer refers to the concept of thermal annealing towards the original high resistive state by Joule heat generated by electric currents over sufficient times. Such a process is qualitatively similar to that of electric fuse where a metal wire melts when too much current flows through it, thereby stopping or interrupting the current. Accelerated by such heating are structural transformations such as thermally activated material diffusion or other type of void nucleation, the rates of which are thermally activated, exp(−w RESET /kT). The activation energy barrier w RESET remains a system parameter at the present level of analysis. The Joule heat elevated temperature T will trigger the process.
The kinetics T(t) of conductive filament temperature will depend on structure details. Skipping the derivation, we show here one relevant example of the adiabatic RESET [45], where heat does not propagate beyond the filament (similar to the known cases of chemical kinetics [46]) Here W ON is the assumed activation energy of resistivity, N is the number of atoms in the filament, and A(t) is the energy deposited during time t. Applying the latter approach to an entire bond where a number M of elements have undergone switching to R < (see section 3.4), we take into account that the latter operate in the current source regime as being in series with much higher (R > ) microscopic resistances in the bond. Power generated in a R < resistance element is then estimated as, where r is the integral resistance of the R > elements, and U L is the voltage across the bond. It follows from equation (20) that triggering a one element RESET will change M → M − 1, which can result in a reverse 'domino' sequence of RESETs through all M low resistance elements or can be conducted sequentially in the pulse RESET regime. Another mechanism of RESET in a partially switched bond takes place in PWP with a degree of ferroelectricity allowing metastable ionic polarizability in response to strong electric fields [47]. Such ferroelectricity is not limited to materials exhibiting macroscopic ferroelectric phase transformations, implying rather the ability to produce metastable ionic displacements commonly described in terms of atomic double well potentials. A variety of materials (such as SiO 2 and HfO 2 ) that are not nominally ferroelectrics possess that property [48,49]. Given that metastable ionic polarizability, the mechanism of RESET evolves through the radial electric field produced by a conductive filament under electric bias [47]. As the bias polarity changes, the earlier created dipole polarization becomes energetically unfavorable, and the energy is gained by eliminating the radial field through the filament rupture.

Examples of functionality
The present paper is aimed to serve as a blueprint for further research and possible implementations of PWP systems where many features remain in their infancy. A cursory description of a few anticipated developments is given next.
Multivalued memory. When high enough voltage is applied between a pair of PWP electrodes, their connecting path will change its resistance depending on the pulse duration as explained in section 3.4 above. The uniquely large number of interelectrode pathways is important here. Given even a modest example of 3 electrodes on each of the faces of a 1 × 1 × 1 cm 3 cube, the total number of electrodes becomes N = 3 × 6 = 18 and the number of perceptive pathways N = N! ∼ 10 16 cm −3 is theoretically greater than that of human cortex resulting in a higher memory capacity than the current 3D crossbar architecture. It will be further enhanced with the functionality of multiple (M 10) records per one cluster bond, so that the total number of records becomes N M 10 17 cm −3 . Various unaccounted factors (crosstalk, leakage, reading errors, material defects) could interfere, and the claim of that superior memory capacity remains to be validated experimentally.
Generation of random numbers. As shown in section 3.1, not-too-spatially-close pairs of electrodes have random resistances uncorrelated with any desired accuracy given in equation (7) when their spacial separation is large enough. Therefore, the resistances of all such pairs of electrodes present a multitude of uncorrelated random numbers. The exponentially great number N is of essence again. The larger the system the smaller degree of correlation between the random numbers can be achieved. For illustration purposes, we consider a single size L a structure (i.e. not very asymmetric), for which we can roughly estimate the typical interelectrode parameters as N ij ∼ N kl ∼ L/a and L ijkl ∼ L in equation (7). With these estimates, equation (7) yields C ∼ exp(−L/a) 1.
More specifically, assuming L ∼ 10 μm sample size with a ∼ 10 nm dimension for a microscopic resistor, our prediction becomes C ∼ exp(−10 3 ) 1. Unfortunately, the latter remarkably low correlation remains more of a theoretical limit than a practical quantity. Indeed, according to equation (6), the coefficient C describes correlation between the resistance deviations ΔR ij = R ij − R ij rather than the resistances themselves. Hence, a prerequisite averaging over multiple electrode pairs in each participating region of the sample must be a part of its practical implementation. Such averaging can be conducted readily by considering conductances G ij = R −1 ij instead of resistances and measuring the integral conductance across the sample. However, it will guarantee only a limited relative accuracy that is not exponentially small and is estimated as C ∼ ΔR 2 ij /( √ N L R ij ) according to the central limit theorem, where N L is the number of electrodes in the area of linear dimension L.
Using the percolation theory estimate from equation (12), ΔR 2 ij / R ij ) ∼ L c /L, and assuming realistic N L ∼ 10-100 yields C on the order of several percent, comparable to but not outstanding against the background of known random generation machines data [51]. We arrive at the conclusion that the estimate in equation (21) represents a theoretical limit and moving closer to that limit remains a challenge of sample optimization.
Matrix-vector multiplication. The measurement based operation of matrix-vector multiplication follows from figure 1. Suppose that A i is the desired product of the vector J j and the matrix F ij . We map these multipliers respectively on the system voltages and conductances. That is achieved by rescaling J j with a certain multiplier (z 1 ) to a convenient interval of electrode voltages E j . Similarly, using a proper multiplier (z 2 ), we rescale F ij so that all its elements fall in the interval of the interelectrode conductances, δG = δR/ R 2 ij with δR from equation (12). The desired product becomes A i = z 1 z 2 I i where I i = j G ij E j is the current through the ith electrode in figure 1. Because the conductance matrix G ij = R −1 ij contains exponentially large number (N 1) of elements covering interval δG, any desired value of G ij can be located at least approximately among the measured conductances with fairly good accuracy without any additional actions. We observe that the role of exponentially large combinatorial number N 1 is most significant in determining the density of interelectrode pathway conductance spectrum allowing approximation of the matrix elements to the accuracy ∼ ΔG/N . After that, applying voltage E j to the electrode j produces a measurable current G ij E j through electrode i, which can be stored e.g. as a partial charge on a certain capacitor C i . Measuring the total of all such partial contributions supplied by electrode i in response to various electrode voltages E j will give the component of vector I i . That procedure can be further improved by choosing different approximate G ij s and using linear regression with all the chosen values to minimize the error.
Brain-like associative learning commonly illustrated with Pavlov's dog salivation experiments (see e.g. reference [52]) is readily implemented utilizing shared portions between bonds of a PWP cluster, such as (1,5) and (1,7) in figure 1. Identifying the 'sight of food' and 'sound' stimuli with signals on the electrodes 5 and 7, predicts that properly and simultaneously triggering both will switch their corresponding pathways (1,5) and (1,7) to a low resistive state making both salivation triggering (through the output on electrode 1). In general, conductive pathways connecting various pairs of the electrodes and sharing the same portion of a PWP cluster will be mutually affected by a single bias-induced change. That demonstrates a single-trial learning model for storage and retrieval of information resembling that of the cortex of the mammalian brain.
Other functionalities based on PCM and RRAM structures for neuromorphic computing appear all attainable with PWP systems. We note that even without utilizing their plasticity, PWPs can serve as high capacity tunable nonlinear reservoirs for reservoir computers. For example, introducing nonohmic (yet volatile) changes in the resistances of pathways (1,5), (11,5), and (5, 7) will produce measurable changes in the resistances of pathways (11,3), (1,7), and (11,7), and the latter will depend on temporal order in which the former changes were introduced. Utilizing the system plasticity will significantly add to that functionality. In fact, our proposed PWPs represent exponentially more powerful reservoir computing systems compared to the one built using a limited number of memristors [50].

PWP metrics and implementations
We briefly mention several metrics of the proposed PWP devices following the nomenclature used for other neuromorphic systems [52]. (1) Dimensions and architecture. The above estimate of a superior information density may be reduced to account for larger physical dimension of a single microscopic resistor. However, even assuming a microscopic resistor of PWP dimension a in the range of tens of microns yields the density ∼ 10 9 cm −3 . (2) Energy consumption. Assuming a PWP structure made of the same materials as the existing PCM and RRAM, we expect its energy efficiency to be superior because of the lack of interconnects requiring costly energy [53]. (3) Operating speed/programming time. Generally, PWP devices RC times are greater than those of nano RRAM and PCM. Like other brain-inspired systems, their computational efficiency will be achieved through the high degree of parallelism. We recall in this connection that the combinatorally huge number of memory units N! is exponentially higher than the number of electrodes N. (4) Multi-level states. Assuming a ∼ 10 nm microscopic resistor and 1 cm device, each bond in a PWP cluster will contain hundreds of microresistors; hence, hundreds of multi-level states per typical bond, at the level allowing robust analog operations. (6) Retention and endurance. PWP systems can be superior to the existing PCM/RRAM based devices because of the lack of multiple interconnects triggering degradation [25][26][27].

Similarity with biological systems
The PWP architecture and functionality have similarities with that of biological neural networks. We recall that the latter operate through individual neurons integrating input signals and firing pulses upon exceeding a certain threshold. The neurons communicate with network by means of synapses that inject ions through their membrane channels then electrically altering adjacent objects. At longer distances, the neuro-synaptic entities are connected with each other through axons providing pathways for electric pulses. A significant degree of randomness and stochasticity is introduced by variations of material parameters, particularly, ion channels whose concentrations and characteristics vary even between nominally similar biological membranes [54][55][56][57].
The property of firing pulses upon accumulating enough electric charge is found in PWP's slow elements (τ > t) that operate as capacitors. Such elements turn into resistors when, through their multiple connections in PWP, they acquire voltages sufficient to make real currents larger than the displacement ones [17]. Such accumulation becomes possible when the signals arrive within a certain time interval thus resembling STDP in biological neural networks. The inherent randomness of PWP correlates with that in biological neural networks.
Note however that the neuron analogy described here pertains to the system functionality rather than its structure. In PWP, slow elements can be associated with any of the system components, depending on the pulse duration. Therefore, the PWP neuron-like elements are in a sense distributed throughout the system and the same element can play the role of a synapse, or axon, or neuron depending on the excitation conditions. Such ability of one and the same element to perform information computation, communication, and storage is characteristic of neuromorphic systems [15].
One aspect of PWP operations ignored in the present work is the inevitable structural degradation limiting the system reliability. All types of degradation known for RRAM and PCM devices, such as recrystallization and diffusion processes can exhibit themselves in PWP compounds. A question of their effect on PWP parameters remains open because the synergy of multiple random elements constituting PWP will not necessarily be affected by possible redistribution of their parameters. With regards to that highly nontrivial question, we can however point at a certain analogy with biological systems where physical degradation can or cannot lead to detrimental consequences depending on the system architecture and functions.

Numerical estimates of PWP parameters
Assuming a ∼ 1 nm and V max /kT ∼ 100 yields L c ∼ aV max /kT ∼ 100 nm. Empirically, the field reliably leading to switching is 1 MV cm −1 . On the other hand, it is estimated by our theory as U 1 /a, with U 1 from equation (18), which yields U L 0.1 V. The latter corresponds to a fairly attainable electric potential drop of 10 kV across 1 cm thick samples. We conclude that multivalued memory by switching in PWP corresponds to rather moderate voltages at least for large PWP systems. Because the earlier estimated ∼ 10 records corresponds to the bond of correlation radius L c ∼ 100 nm, the number of records per 1 cm thick sample is estimated as ∼ 10 6 . The question of their temporal overlaps remains to be addressed based on statistical properties of percolation clusters.
For small PWP, using the same parameters and sample size L L c one obtains ∼ 10L/L c possible transitions, which however will be all distinctly different. Note that for 10 nm samples that estimate predicts one switching event, which on average is what is observed in the current RRAM and PCM devices.

The role of randomness
A characteristic feature of our proposed PWPs is that they are inherently random and possess exponentially broad spectra of electric parameters. As such, they are suggestive of probabilistic algorithms that employ a degree of randomness as part of its logic [58][59][60].
The uniqueness of our approach is that it proposes an inherently 'random structure', with which even a logically deterministic algorithm becomes randomized. The general power of such approach remains to be tested yet, although some examples do point at its potential: superiority of randomly-wired networks [23,24] and intentionally randomized reservoir computers [18]. Biological neural networks provide of course an ultimate example of neuromorphic computing leveraging randomness (we are not built of well controlled silicon nano-chips).
From that perspective, our approach proposes a practically implementable machinery for testing the potential of 'physically randomized neuromorphic systems'. That goal can be achieved both through computer modeling of PWP that benefits from the earlier developed percolation modeling algorithms, and through direct material implementations, such as e.g. chalcogenide films with multiple electrodes.

Conclusions
We have introduced the PWP systems that, while being essentially random, exhibit various neuromorphic functionalities and similarities to the cortex of mammalian brain. These systems demonstrate rich physics, understanding of which remains in its infancy.
The outstanding challenges in theory go far beyond the standard percolation paradigm including statistical aspects, microscopic phase transformations, heat transfer, AC propagation and signal cross-coupling in random systems with multiple co-existing interfaces. Both analytical research and numerical modeling will help to better understand PWP operations and functions. Multiple material bases can be used to fabricate PWP devices, ranging from the known PCM and RRAM materials to the nanocomposites which beg to be further explored for memory and other neuromorphic applications. The above noted functional similarities with biological neural networks suggests the need for further investigations.
To avoid any misunderstanding we should emphasize that the present work is limited to the basic theory of PWP systems. It only tangentially addresses the related neuromorphic applications that can be tried at both the levels of software implementations and real material devices. Obvious experimental projects potentially triggered by this work include pulse regime nonohmicity in disordered systems, statistical properties of percolation clusters with multiple electrodes, and reservoir computing based on percolation systems.

Data availability statement
All data that support the findings of this study are included within the article (and any supplementary information files).