A kinetic Monte Carlo Approach for Boolean Logic Functionality in Gold Nanoparticle Networks

Nanoparticles interconnected by insulating organic molecules exhibit nonlinear switching behavior at low temperatures. By assembling these switches into a network and manipulating charge transport dynamics through surrounding electrodes, the network can be reconfigurably functionalized to act as any Boolean logic gate. This work introduces a kinetic Monte Carlo-based simulation tool, applying established principles of single electronics to model charge transport dynamics in nanoparticle networks. We functionalize nanoparticle networks as Boolean logic gates and assess their quality using a fitness function. Based on the definition of fitness, we derive new metrics to quantify essential nonlinear properties of the network, including negative differential resistance and nonlinear separability. These nonlinear properties are crucial not only for functionalizing the network as Boolean logic gates but also when our networks are functionalized for brain-inspired computing applications in the future. We address fundamental questions about the dependence of fitness and nonlinear properties on system size, number of surrounding electrodes, and electrode positioning. We assert the overall benefit of having more electrodes, with proximity to the network's output being pivotal for functionality and nonlinearity. Additionally, we demonstrate a optimal system size and argue for breaking symmetry in electrode positioning to favor nonlinear properties.


I. INTRODUCTION AND MOTIVATION
Brain-inspired computing represents an approach to advancing computation in machine learning applications by emulating the information processing of biological neural networks [1].Compared to von-Neumann computing, it leverages massively parallel operation without separating memory and processing [2].The implementation of these brain-like infrastructures aim to overcome limitations of traditional hardware and ultimately reduce energy consumption of current machine learning applications [3,4].Nanoparticle (NP) networks offer a promising avenue in this field.One current approach uses percolating scalefree NP networks with small-world properties [5][6][7][8].The intrinsic architecture of these networks includes conductivity dynamics which are crucial for computational tasks [9,10].Tunnel gaps forming between adjacent areas or clusters of nanoparticles act as nonlinear units or memristors, with their resistance exhibiting nonlinear responses to alterations in conductance.While the origin of the nonlinearity here arises from the connections of large clusters of NPs and not the particles themselves [6], our approach is different.Our research in this domain centers single electronics within gold NP networks [11,12].Each NP serves as a conductive island, tunnel-coupled with its neighboring NP by an insulating organic molecule.When a NP is charged, its excess charges can tunnel to an adjacent NP only when the potential difference between them surpasses the repelling force exerted by the Coulomb energy.This Coulomb blockade effect imbues the charge tunneling dynamics in between NPs with a nonlinear activation function [13].The network is surrounded by electrodes, which serve as network inputs, outputs or controls (see FIG. 1).The latter manipulate the network's potential landscape and its charge tunneling dynamics, enabling the mapping of various functionalities to the input and output electrodes.Previous research has demonstrated the network's input to output dependence to reconfigurably function as any type of Boolean logic gate when applying suitable control electrode voltages [14,15].FIG. 1 shows a circuit diagram of a small NP network and electrode setup.The integration of theoretical frameworks and simulation methods is of key importance in order to enable a closer understanding and subsequent optimization of these systems.In this paper, we use a kinetic Monte Carlo (KMC) approach to simulate electronic charge tunneling dynamics within the NP network.Prior research has already demonstrated the effectiveness of the KMC approach for studying single electronics devices [16][17][18][19].In particular, Van Damme et al. [20] functionalized small networks of 16 nanoparticles to act as any Boolean logic gate.Our work extends previous analysis in several directions.First, it addresses the influence of system size (number of NPs) as well as number and location of electrodes to optimize the chance to observe Boolean logic functionality.This information is relevant for the future design of such NP networks.Second, from a methodological perspective we introduce quantitative indicators for negative differential resistance and nonlinear separability.These nonlinear properties are imperative for a broad spectrum of brain-inspired computing [21] and machine learning applications [22].Here, we explicitly show that these indicators are quantitatively related to the ability of the network to display Boolean logic functionality.Third, due to the availability of fast array-based computations we highlight the efficiency of the implementation of the KMC approach.
The system of [14] consists of Au NPs which are interconnected by 1-octanethiols as an insulating organic molecule.The network is placed on a highly doped Si substrate with an insulting SiO 2 top-layer and is surrounded by 8 Ti/Au electrodes.From an electrostatic point of view, we just speak of a network of conductive islands interconnected by resistors and capacitors in parallel.Electrodes are understood as voltage sources U E .FIG. 1 shows a circuit diagram of a 9-NP network in a regular grid-like connection topology.We derive the capacitance values of the network using the image charge method (see Supplementary Section VII A).
The capacitance C ij between island i and j depends on the NP radius r i , r j , the distance between both NPs d ij and the permittivity of the insulating material ϵ m .The self capacitance C i,self for the isolated island i close to the SiO 2 layer depends on its radius r i and the permittivity of the insulating environment ϵ SiO2 .Even though the simulation tool allows to setup variable radii or distances, for this work values are set constant with r = 10 nm and d = 1 nm, respectively, following [14].For a given connection topology, we setup the capacitance matrix C where diagonal components (C) ii represent the sum of all capacitance values attached to NP i and its self capacitance, while off-diagonal components (C) ij represent the cross capacitance in between NP i and j.In this work, we will not cover the impact of having variable connection topologies.Therefore, for all results we align NP in a regular grid like two-dimensional topology, such as in FIG. 1.

B. Single-Electron Tunneling
The network consists of conductive NPs interconnected by insulating molecules.The insulating molecules serve as tunnel barriers and functionalize the network.Excess charges residing on NP i tunnel to an adjacent NP j if their associated potential ϕ i exceeds the repelling force resulting from the Coulomb energy E C,j = e 2 (C)jj of NP j.For a NP with 4 next neighbors at maximum we receive a charging energy of about E ≈ 16 meV.This behavior leads to a switch-like nonlinear activation function for the charge tunneling dynamics where charges may or may not tunnel to their nearest neighbors corresponding to the current network's potential landscape.We define the network state ⃗ q(t) as the number of excess charges residing on each NP at a given time step.Then, the network's potential landscape is calculated as with C −1 as the inverse of the capacitance matrix.The internal electrostatic energy is described as When voltages are applied to the surrounding electrodes, charges either enter the network or getting drained from it.The Helmholtz's free energy determines the tunneling process and consists of the difference of total energy stored inside the network E and work W done by external electrodes In its current state, the network will enter its next state associated to a change in free energy, either by a charge tunneling event from NP i to j or by a charge tunneling event from NP i to electrode E with elementary charge e.The change in free energy caused by a tunnel event serves as a measure of probability for this specific event with tunnel rate where R ij denotes the resistance for jump path i → j and T the network's temperature.If we apply constant voltages to all electrodes and execute at maximum about 10, 000 charge tunneling events, the network will eventually settle in an equilibrium with constant electric currents entering or exiting the system via the electrodes.We revisit Eq. ( 6) and see that increasing the network's temperature T leads to drastically larger rates.Then, differences in the potential landscape are negligible and thermal excitation will be the dominant factor for the charge tunneling through the network.Overall, the device looses its functionality due to the linearization of the charge tunneling dynamics, i.e. its nonlinear activation functions (see Supplementary Section VII B and FIG.S2).Therefore, we have to choose the temperature, so that condition E C > k B T is true.In this work we will stick to a temperature value of T = 0.28 K as in [14] which easily satisfies the latter condition with k B T ≈ 2.5 • 10 −2 meV much smaller than the maximum charging energy of E ≈ 16 meV.Additionally we also want to neglect coherent quantum processes, i.e. co-tunneling.For this we have to choose the tunneling resistance R ij in between NPs to be much higher than the quantum resistance of R t = h e 2 ≈ 25.8 kΩ.In our simulations we set the tunneling resistance for all tunneling processes to a constant value of R = 25 MΩ, which is sufficient to assume charges to be confined on our NP islands.

III. MODELLING AND SIMULATION
Based on the findings in the previous section, we explain how to build a KMC simulation tool to model the charge tunneling dynamics in NP networks.We also explain how to functionalize the network to mimic Boolean logic gates.In the last part we define the fitness value as a way to evaluate the quality of the network behaving a a particular Boolean logic gate.Based on this quality factor we derive Q NDR and Q NLS as new quantities to measure the properties of negative differential resistance (NDR) and nonlinear separability (NLS).

A. Kinetic Monte Carlo Method
For a given network topology of N NP nanoparticles, initialized capacitance matrix C and set of electrode volt-NanoNets Simulation Tool Compute Γij(∆Fij)
In the next step we compute the cumulative distribution function (CDF) of all tunnel rates CDF Γ (n).The CDF consists of N T steps, one for each charge tunneling event n with its last step CDF Γ (n = N T ) as the total rate constant k tot .Now, for the KMC procedure we first need to sample two random numbers r 1 and r 2 and find n where As n was selected, we execute its corresponding charge tunneling event from NP i to NP j.This will leave the current state ⃗ q(t 1 ) towards ⃗ q(t 2 ).Afterwards we update the time scale Since only two elements of the state vector have been updated, for the next potential landscape we can just evaluate where ⃗ j is a vector including zeros and elementary charge −e at position i and +e at position j.NanoNets shows an overview of the KMC procedure.We implemented the simulation tool in Python using NumPy [23] and Numba [24] for fast array-based computations (for NanoNets GitHub see [25]).The process of updating potentials and selecting the next event from the sorted CDF is achieved in about O(log(N )).However, since the whole system is coupled via the potential landscape, we have to calculate all tunnel rates after each tunneling event, thus changed potential landscape.This leads to a time complexity of about O(N ).FIG. 2 shows the process time for the major parts in the KMC procedure.
Process time for calculating all tunnel events (blue curve) and advancing into the next step by selecting an event (red curve).

B. Computing Functionalities
In this paper we want to investigate computing functionalities in the form of Boolean logic gates as a systematic property of the NP network.For Boolean logic we need two inputs and one output.Inputs, either set on/1 or off/0, produce an output signal depending on the underlying logic operation.In our case, two electrodes will always serve as input electrodes.These electrodes are set to U I ∈ {0.0, 10.0} mV.Accordingly, there is an output electrode, which is always grounded, thus U O = 0.0 V.At the output we will measure the electric current as a function of time t passed in the KMC procedure and number of elementary charge jumps from network to output N Net→O and from output to network N O→Net .All other electrodes are called control electrodes with voltages set inside a range of U C ∈ [−50.0,50.0] mV.The system's phase space is spanned by N C electrodes with each point in phase representing a set of four different electric currents (I 00 , I 10 , I 01 , I 11 ) corresponding to the four possible input electrode voltage combinations.When increasing the network size, the electrode voltages drop across an increasing amount of NPs.To ensure comparable results, we scale the voltages of control and input electrodes to prevent the electric current at the output electrode from simply converging to zero as the network expands.To achieve this, we conducted experiments to assess the dependencies between input electrode voltage and output electrode electric current (I-V curves) across multiple systems.Through this analysis, we derived a scaling factor that maintains the electric current at a consistent level when multiplied by the initial voltage value, see Supplementary Section VII C and FIG.S3.Summarizing the simulation procedure, we first initialize a network of regular grid-like topology, electrode positions and electrostatics.Then we apply constant voltages to all electrodes and evolve the system into its equilibrium state.Afterwards, we start to track the time scale and the number of charge jumps exiting and entering the system via the grounded output electrode.This allows us to compute the electric current I and its relative uncertainty u I .The simulation ends when the KMC procedure reaches u I = 0.05 = 5% uncertainty or when ten million charge tunneling events have been executed.

C. Analysis Methods
For a single simulation run, we sample 20, 000 different control combinations, resulting in a set of four electric currents (I 00 , I 10 , I 01 , I 11 ) for each combination, totaling 80, 000 currents.To characterize nonlinear properties we introduce the three parameters M l /M r signify the increase in output current upon altering the first/second input voltage, whereas X quantifies the cross-correlation between the two input voltages.In qualitative terms, M l and M r can be interpreted as an effective mobility of the output current concerning one of the two input voltages, while X denotes a measure of nonlinear coupling between both inputs.Our objective is to evaluate the network's ability to be reconfigured into each Boolean logic gate using the fitness function as a quality metric for the accuracy with which a given set of currents (I 00 , I 10 , I 01 , I 11 ) represents a specific gate [14].
The fitness comprises the slope or signal m, the mean squared error or noise M SE, and the absolute offset |c|.
The signal is defined as where ON /OF F represents the set of electric currents corresponding to the gate's on/off currents, and |ON |/|OF F | is the cardinality.For example, for an AND gate, ON = {I 11 } and OF F = {I 00 , I 10 , I 01 }.
The noise is defined as ] and the offset as c = I of f .In total, the fitness consists of a signal-to-noise ratio where logic gates with small offsets are more favourable from the experimental perspective as expressed by δ > 0. Theoretically it is easier to choose δ = 0.In Supplementary Section VII D and FIG.S5 we cover the influence of δ > 0.
For δ = 0 we can exactly rewrite Eq. ( 12) for the individual gates in form of the quantities introduced in Eq. ( 11) Our goal is to estimate the first and second moment of the fitness distribution for each gate in terms of the statistical properties of M l , M r , and X.For this purpose we use two approximations: Firstly, we approximate < F (x) > by F (< x >) and, secondly, we neglect cov(M l , X) and cov(M r , X).Both terms are negligible relative to cov(M l , M r ), as demonstrated in Supplementary Section VII E and FIG.S6.Furthermore, since we only consider regular grid-like networks with a symmetric electrode connection, we naturally have ⟨X⟩ = 0 as well as identical distributions for M l and M r when averaging over a sufficiently larger number of control voltages.Therefore, we may abbreviate < M l >=< M r >=< M > and < M 2 l >=< M 2 r >=< M 2 >.Using these approximations, for the AND, OR, NAND, and NOR gates we can write and for the XOR and XNOR gates Nonlinear behavior is in particular relevant for the NAND, NOR, XOR, and XNOR gates.For these gates we want to define measures which reflect the likelihood to find gates with sufficiently high fitness.For the NAND/NOR gates we consider the ratio < F NAND/NOR > / Var(F NAND/NOR ).Both, if < F NAND/NOR > is not too negative and/or the standard deviation Var(F NAND/NOR ) is sufficiently large, there is a higher probability that positive fitness values occur, i.e. the system starts to display the respective functionality.On this basis we define In the last step we have neglected the Pearson correla-tion between M l and M r .
C I 1 The properties of the tanh-function mathematically restrict Q NDR to values between 0 and 1.Since ⟨M ⟩ is always positive, the actual upper bound is 0.5.For values of Q NDR close to zero, the possibility of NAND and NOR gates is very small.Qualitatively, this measure allows us to assess the nonlinear feature of negative differential resistance (NDR) required to achieve NAND/NOR gates.We mention in passing that the variance in Eq. ( 14) becomes larger for stronger correlations between M l and M r .This increases the probability to find gates (AND, OR, NAND, NOR) with higher fitness values.From Eq. ( 15) we may conclude that XOR and XNOR gates have the same statistical properties.Since in this case the likelihood to find well-performing exclusive gates with high fitness values is directly related to the size of the variance of the fitness distribution, we can introduce as a general measure for nonlinear separability (NLS), strongly reflecting the size of the nonlinear coupling X between both inputs.

IV. RESULTS
In our study of nanoparticle networks manipulated by surrounding electrodes, key design features include the number of nanoparticles (N NP ) and the location and quantity of control electrodes (N C ).In this section, we first analyze the dependence of logic gate fitness (F ) on the location and number of control electrodes.
Given the significant impact of electrode positioning, the subsequent section addresses an increase in system size while considering two distinct electrode positioning setups.All results are contextualized within the framework of nonlinear properties.Additional insights into parameter dependence on N NP and N C can be found in Supplementary Section VII E and FIG.S7.

A. Number and Positioning of Control Electrodes
In this section, we explore the impact of the number of control electrodes and their positioning on network functionality, focusing on regular grid-like networks as discussed in Section II A. We use a 7 × 7 grid of NPs (N NP = 49) with one output and two input electrodes.To ensure that the majority of the system resides between input and output, we connect the output to a corner NP, and both inputs to the opposite edges.We start by adding control electrodes near both inputs as we move closer to the output.These simulations are denoted AN C , where N C represents the number of controls.Subsequently, controls are removed from the system, starting with those near both inputs.These simulations are denoted BN C .The upper sketches in FIG. 3 illustrate these simulated systems.This procedure allows us to study the effect on network functionality for varying control numbers and control positioning.
For each network, we sample 20, 000 control voltage combinations, producing 20, 000 sets of electric currents (I 00 , I 10 , I 01 , I 11 ) using KMC.Utilizing Eq. ( 12 for each Boolean logic gate at δ = 0.The fitness distributions for each gate across different networks are depicted as box plots in FIG. 3 a), providing insights into the relationship between functionality and number of control electrodes.First, in agreement with theoretical expectation, see Eq. ( 14) and ( 15), we observe similar distributions for the respective AND/OR, NAND/NOR, and XOR/XNOR pairs.Furthermore, the distributions of the first two pairs are just mirror images of each other.Any remaining variations in < F > and Var(F ) within a pair can be attributed to covariance values between M i and X, which we did not consider in Eq. ( 14) and Eq. ( 15).The impact on fitness when altering the number or position of control electrodes is significantly pronounced.We observe the most dramatic effect for controls in close proximity to the output electrode.The measures of nonlinear behavior, i.e., Q NDR and Q NLS , exhibit much larger values, see FIG. 3 b).Theoretically, one expects a strong correlation between Q NDR and the availability of high-fitness NAND/NOR gates as well as between Q NLS and the number of exclusive gates XOR/XNOR.Indeed, as seen from the box plots, a significant correlation is indeed observed.As a secondary effect, it turns out to be helpful to use the maximum number of control electrodes N C = 9.Only in the case of AND/OR gates, a reduction in the number of controls generally leads to an increase in average fitness, while retaining electrodes near the output results in the most significant second moment.Moreover, we observe a strong impact on the correlation coefficient between M l and M r .The presence of control electrodes adjacent to the output electrode as well as an increase in the number of control electrodes significantly reduces this correlation, going along with a stronger independence of I 01 and I 10 .
As the results in FIG. 3 suggest a substantial impact of the positioning of controls relative to the output or inputs, we performed two additional simulations for a network of size 7 × 7. First, we randomly adjusted the voltages of all 7 electrodes and determined the resulting correlation between individual voltages and the output current.The Pearson correlations are illustrated in FIG. 4. In alignment with previous observations, electrodes in close proximity to the output electrode (specifically E 5 /E 6 ) exhibit a strong impact via notable correlations between electrode voltages and output electric current I.This observation is further supported by the current-voltage dependencies, which are strongest when altering electrode voltages near the output (see FIG. S2).Second, we selected inputs at all possible electrode combinations (E i , E j ) and sampled multiple control voltage combinations for any possible input electrode voltage combination {(E i , E j ), (E i + ∆, E j ), (E i , E j + ∆), (E i + ∆, E j + ∆)} with ∆ = 10 mV.The last two plots in FIG. 4 depict the nonlinear features of NDR and NLS for each input position combination.If a single input is positioned adjacent to the output at E 5 or E 6 , NLS is lost entirely due to the strong correlation with the electric current.For either NDR or NLS, optimal results are achieved by breaking symmetry and attaching one electrode at the corner opposite to the output electrode at position E 0 , with the second input placed elsewhere.In fact, there is not a significant difference in attaching the second input at E 1 , E 2 , E 3 , or E 4 .The advantage of having an asymmetric electrode connection may indicate a potential functionality gain for disordered networks in general.

B. System Size
In this section, we investigate the influence of system size on the network functionality of regular grid-like NP networks.Following the experimental device proposed in [14], each network is surrounded by a configuration Simulation results for negative differential resistance (QNDR), nonlinear separability (QNLS), and input mobility correlation (corr(M l , Mr)) in networks with varying numbers of NPs (NNP).These networks are surrounded by five control electrodes (C), two input electrodes (I1,I2), and one output electrode (O).Two setups are considered: Setup A maintains same relative distances between all electrodes, while Setup B keeps the electrodes at fixed positions. of 8 electrodes.We connect the output to one corner of the network, and both inputs to the opposite edges to ensure that most NPs are between the input and output.This leaves five control electrodes to be attached at the remaining corners or edges.Considering the strong relation to electrode positioning indicated in the previous section, we distinguish between two setups, as displayed in FIG. 5 a).In Setup A, we maintain consistent relative distances between all electrodes, connecting them either to a network corner or to the center of a network edge.In Setup B, the two electrodes, adjacent to the output, keep this minimum distance independent of system size.Each simulation covers a specific system size, ranging from a 3 × 3 grid of NPs (N NP = 9) to a 16 × 16 grid of NPs (N NP = 256).We conduct 20, 000 simulations for each system size, sampling control voltage combinations and producing sets of electric currents (I 00 , I 10 , I 01 , I 11 ).
In FIG. 5 b), we break down the fitness into its key nonlinear components: Q NDR , Q NLS , and corr(M l , M r ).There is a noticeable difference between both setups.In Setup B, a saturation behavior is observed for each nonlinear component.For N NP > 100, maximum NDR and NLS are achieved, and both inputs start to become fully independent.In fact, Q NDR even reaches its theoretical upper bound Q NDR = 0.5 as described in Section III C. In contrast, Setup A shows a less pronounced dependence.For systems seeking nonlinear separability, those with more than 50 NPs are not preferred.There is a decline in performance for increasing numbers of NPs, suggesting that smaller systems tend to perform better in terms of nonlinear features, provided there is a minimum of about 16 NPs.Regarding the correlation between both mobility values M l and M r , we observe that independence between input electrodes is never achieved, with corr(M l , M r ) ≈ 0.5 for all networks.As shown in Supplementary Section VII D and FIG.S4 the fitness distributions for AND, NAND, and XOR follow the same trend.Similar to the findings in the previous section, we observe that control electrodes in close proximity to the output electrode play a pivotal role in overall network functionalities, particularly in terms of nonlinearity.As the system size increases in Setup A, the distance between the output and its closest control electrodes also increases, leading to a negative impact on performance.Despite achieving better results in larger systems with Setup B, it is crucial to note that in practical experiments, the network configuration may not resemble Setup B. Therefore, we argue that in experiments, one should be cautious about making the network too large, even though theoretically, saturation in nonlinearity should occur.Optimal results are obtained when aiming for a maximum of about N NP = 100, especially if electrodes can be densely packed.Otherwise, even a smaller system may prove to be more beneficial.It is noteworthy that the saturation observed in FIG. 5 occurs only for δ = 0, as Q NDR and Q NLS have been derived under this assumption (see Section III C).For larger values of δ, indicating a contribution of the logic gates offsets to the fitness value, we deviate from the saturation behavior, showing a clear maximum in fitness at about N NP = 100.Considering the preference for logic gates without an offset in the experimental setup, this observation further supports the identification of an optimal system size at or below 100 NPs.

V. CONCLUSION
We have developed a versatile and efficient kinetic Monte Carlo simulation tool designed to model charge tunneling dynamics in NP networks.Building upon established concepts of single electronics, the model has been specifically tailored for NP networks.The tool allows for the investigation of larger networks, comprising hundreds of NPs, within reasonable processing times.Within our simulation framework, electrodes can be connected to the network at arbitrary positions and designated as input, output, or control electrodes.We have successfully mapped Boolean logic gate functionality to the dependence between input and output electrodes.Adjusting the control electrode voltages, the charge tunneling dynamics is manipulated, enabling the same network to be configured into any desired logic gate.We assessed the quality of the network in functioning as a logic gate through a fitness function.Increasing the number of control electrodes consistently enhances fitness, with controls in close proximity to the output electrode proving most pivotal.When altering the positions of inputs and controls, our results demonstrated a loss of functionality when inputs are positioned in close proximity to the output.Additionally, we observed that an asymmetric connection of inputs could be advantageous.We further related the fitness of each logic gate to a novel set of variables, enabling the derivation of measurement tools for key qualitative nonlinear features such as negative differential resistance and nonlinear separability.These features, while essential for Boolean logic, also hold relevance for other classification or machine learning functionalities.Upon increasing the network's size, we unveil optimal nonlinear properties for a minimum of about 100 nanoparticles.However, this saturation needs electrodes to remain close to the network's output.If an increase in size is accompanied by an enlargement of the area between the output and its neighboring electrodes, a clear maximum is achieved below 100 nanoparticles.The same holds true when the logic gate offset is also considered in the fitness.In this context, we argue that smaller networks, 100 nanoparticles at maximum, are advantageous in experimental setup where ensuring the correct electrode connection distance (i.e.nanoparticles in between) may be challenging.The results in this work describe the impact of constant voltage signals on the network response.An important next step involves time-varying signals, aligning the input time scales with the intrinsic time scales of the network states.This should allow us to explore memory effects, interpreting the network as a dynamical system due to its recurrent connection topology.Subsequently, the potential application of reservoir computing allow us to utilize NP networks for temporal signal processing, encompassing tasks such as time series forecasting and approximating not only functions, as Boolean logic, but also dynamical systems.Given the findings of this work on the influence of disorder and asymmetry, the study of disordered network topologies or variable NP sizes, resistances or materials may also provide new interesting facets of brain-inspired applications.

VI. ACKNOWLEDGEMENT
First we want to find the potential ϕ(r) inside a grounded conducting sphere (NP) with radius R using the image charge method (see.FIG.S1).We introduce an image charge q ′ for the charge q positioned at r q relative to the origin of the sphere with distance s q to an arbitrary surface point on the sphere.The image charge is positioned along the line between the sphere center and charge q at r q ′ and distance s q ′ .We have to meet the boundary condition of zero potential on the spherical surface, i.e.
For the two specific surface points aligned on the same axis as q and q ′ , we get the vanishing potentials for the following relationships: Now, for an arbitrary surface points on the sphere with distances to the first charge s q and second charge s q ′ , we get vanishing potentials when inputting into Eq.( 18) with Eq. ( 20).With these results, one can now calculate the total capacitance of two equally shaped spheres (see.FIG.S1).Again using the image charge method, both spheres will be replaced by charges while maintaining the surfaces as equipotentials.The charge q makes the left sphere an equipotential, disturbing the right sphere's potential.The charge q ′ leads to a zero equipotential of the right sphere but disturbing the left.The image of the image q ′′ compensates for q ′ and so on.This procedure leads to the series r q r q' q q' s q s q' θ q q' q'' q''' r q r q'' r q''' r q' FIG.S1.Sketch illustrating the image charge method applied to a single grounded conducting sphere and two adjacent conducting spheres.
and with the potential ϕ = q 4πϵR to the total capacitance Finally, separating the total capacitance in mutual capacitance C ij between NPs i and j and self-capacitance C i,self of NP i, and substituting the distance between both spheres into Eq.( 22), the final equations are achieved, with r q → 2r + d as NP radius r and spacing in between NPs d.

B. Current voltage dependence
For a fixed N NP = 49 network with 8 electrodes, we measured the output electrode electric current (I) with respect to the input electrode voltage (U ).The output electrode, along with all other electrodes except the input, is grounded.The middle plot in FIG.S2 shows current-voltage dependencies obtained by varying the position of the input electrode.We observe a significant increase in current after surpassing the blockade regime when attaching the input close to the output at E 5 .Accordingly, repositioning the input farther away at E 1 or E 2 results in a substantial decreases in current, indicating that most of the electric current might be drained at E 5 .For input at E 0 , the slope is the shallowest.The right plot in FIG.S2 maintains the input at E 0 and uses different colors to represent various network temperatures.Here, it becomes evident that increasing the temperature above approximately 20 K leads to a linearization of the current-voltage dependence, indicating the loss of functionality due to the loss of the network's nonlinear activation functions.Both results are in well agreement with the measurement results from [14].

C. Voltage re-scaling based on System Size
As there is a voltage drop across the entire network, we have to scale the magnitude of electrode voltages according to the number of NPs.Without such scaling, voltage ranges that are above the Coulomb blockade regime for a  small network may not correspondingly surpass blockade for a larger network.This would result in decreasing electric currents with increasing numbers of NPs, not allowing to assess the true impact of system size but merely emphasizing the impact of vanishing network potentials.To address this, we conducted electric current measurements (I) of the output electrode for networks with 8 electrodes.In all network setups, two input electrodes were positioned at opposite edges of the output electrode, both set to the same voltage (U ), while all other electrodes remained grounded.The right plot in FIG.S3 illustrates the I to U dependence for various numbers of NPs (N NP ).Simultaneously, the left plot presents the factor by which U should be multiplied to maintain electric currents of the same magnitude.The current values were normalized to those of the 7 × 7 network.For all simulations, we then applied the N NP -specific factor to the sampling ranges of either U I ∈ {0.0, 10.0} mV or U C ∈ [−50.0,50.0] mV.In FIG.S3 we also clearly see the phenomenon of negative differential resistance, as increasing the voltage reduces the output electric current, at least for a network of about 100 nanoparticles at minimum.

D. System Size and δ-dependence
Here, we present additional insights into the dependence of fitness (F ) on system size.Box plots in FIG.S4 illustrate fitness values for AND, NAND, and XOR logic gates across networks with varying numbers of NPs (N NP ).We compare the results for both setups, A and B, as depicted in FIG. 5 a).Consistently, we observe similar behavior on average between F NAND and Q NDR , as well as between F XOR and Q NLS .The behavior of F AND is complementary to that of F NAND .

5 FIG. 1 .
FIG. 1. Circuit diagram of a gold nanoparticle network.Nanoparticles (golden dots) are tunnel-coupled with each other by organic insulators.The network represent a regulargrid like connection topology.A NP-NP connection is represented by a capacitor and resistor in parallel.The network is placed on a highly doped Si substrate with an insulating SiO2 top-layer (not shown).The nanoparticles are surrounded by 8 Ti/Au electrodes applying voltages U1, . . .U8.

FIG. 3 .
FIG.3.a) Box plots of fitness values for all logic gates for networks with varying numbers of control electrodes.The length of a bar is equal to two standard deviations.Each box corresponds to one of the network diagrams shown above.The networks are projected along the x-axis in first ascending and then descending order of control electrode numbers.b) Measure of negative differential resistance (QNDR), nonlinear separability (QNLS), and Pearson correlation corr(M l , Mr) across variable numbers of control electrodes.
FIG. 5.Simulation results for negative differential resistance (QNDR), nonlinear separability (QNLS), and input mobility correlation (corr(M l , Mr)) in networks with varying numbers of NPs (NNP).These networks are surrounded by five control electrodes (C), two input electrodes (I1,I2), and one output electrode (O).Two setups are considered: Setup A maintains same relative distances between all electrodes, while Setup B keeps the electrodes at fixed positions.
FIG. S2.Current-voltage (I-V) dependence for a network of 49 nanoparticles with 8 electrodes.A voltage is applied to the input electrode, while all other electrodes are grounded.The left figure shows I-V curves with the input electrode positioned at various locations.The right figure shows I-V curves for a fixed input position under different operating temperatures (T ).
FIG. S3.Voltage scaling for different numbers of NPs (NNP).The factor determines how the control sampling range and input states are adjusted depending on the network size.The factor results from the current-voltage (I-V) dependence.When a voltage (U ) is multiplied by the factor each system produces approximately the same electric current (I).