The system of (Bose et al., 2015) consists of about 100 single Au NPs (20 nm diameter) 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₂ 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 . Figure 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 S1). The capacitance C _ij between island i and j depends on the NP radii 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₂ layer depends on its radius r _i and the permittivity of the insulating environment ϵ SiO 2 . 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 (Bose et al., 2015).

For a given connection topology, we define 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 Figure 1, to isolate the upcoming results from additional disorder effects. Although we measure deviations in our system from the experiment due to the absence of disorder, still the overall statistical trends remain consistent (see Supplementary Section S4). Furthermore, the general statistical concepts to characterize non-linear behavior in the present type of systems are sufficiently general to be applied to the disordered case as well. Finally, it should be stressed that due to the generally different choice of control voltages, the overall system is not fully symmetric even for a symmetric network.

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 ) j j of NP j. For a NP with four next neighbors at maximum we find 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 ϕ ⃗ t = C − 1 ⋅ q ⃗ t , (1) with C ⁻¹ as the inverse of the capacitance matrix. The internal electrostatic energy is described as E = 1 2 q ⃗ t ⋅ ϕ ⃗ t . (2) When voltages are applied to the surrounding electrodes, charges either enter the network or getting drained from it. The Helmholtz 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 F = E − W . (3) 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 Δ F i j = e ϕ j − ϕ i + e 2 2 C − 1 i i + C − 1 j j − 2 C − 1 i j , (4) or by a charge tunneling event from NP i to electrode E Δ F i E = e U E − ϕ i + e 2 2 C − 1 i i , (5) 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 Γ i j = − Δ F i j e 2 R i j ⋅ 1 − exp − Δ F i j k B T − 1 , (6) 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 around 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. In this work we will not distinguish between charge tunneling events in between nanoparticles or in between nanoparticles and electrodes except for the above difference in free energy. We argue that this is a fair approximation, as also in the experimental situation the tunnel barrier with the electrodes is formed by an insulating layer of molecules.

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 S2; Supplementary Figure 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 (Bose et al., 2015) 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.

For a given network topology of N _NP nanoparticles, initialized capacitance matrix C and set of electrode voltages U _E , we firstly initialize the network’s state q ⃗ = 0 ⃗ and potential landscape ϕ ⃗ = 0 ⃗ . Then for all possible tunneling events we calculate the change in free energy by Eqs 4, 5 and tunnel rates by Eq. 6.

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 x ₁ and x ₂ and find n where C D F Γ n − 1 < x 1 ⋅ k tot ≤ C D F Γ n . (7)

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 t 2 = t 1 − ln x 2 k tot . (8)

Since only two elements of the state vector have been updated, for the next potential landscape we can just evaluate ϕ t 2 = ϕ t 1 + C − 1 ⋅ j ⃗ , (9) where j ⃗ is a vector including zeros and elementary charge −e at position i and +e at position j. Figure 2 shows an overview of the KMC procedure. We implemented the simulation tool in Python using NumPy (Harris et al., 2020) and Numba (Lam et al., 2015) for fast array-based computations (for NanoNets GitHub see (Mensing, 2023)). 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 also landscape. This leads to a time complexity of about O(N). Figure 2 shows the process time for the major parts in the KMC procedure.

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 I O = e ⋅ N Net → O − N O → Net t (10) It enters the time t passed in the KMC procedure and the 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 ₀₀, I ₁₀, I ₀₁, I ₁₁) 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 simulations 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 S3; Supplementary Figure 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 times 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 = 5% uncertainty or when ten million charge tunneling events have been executed.

For a single simulation run, we sample 20,000 different control combinations, resulting in a set of four electric currents (I ₀₀, I ₁₀, I ₀₁, I ₁₁) for each combination, totaling 80,000 currents. To characterize nonlinear properties we introduce the three parameters M l = 1 4 I 11 − I 01 + I 10 − I 00 M r = 1 4 I 11 + I 01 − I 10 − I 00 X = 1 4 I 11 − I 01 − I 10 + I 00 . (11)

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 configured into each Boolean logic gate using the fitness function F = m M S E + δ ⋅ | c | (12) as a quality metric for the accuracy with which a given set of currents (I ₀₀, I ₁₀, I ₀₁, I ₁₁) represents a specific gate (Bose et al., 2015). The fitness comprises the slope or signal m, the mean squared error or noise MSE, and the absolute offset |c|. The signal is defined as m = I o n − I off = ∑ i ∈ O N I i | O N | − ∑ j ∈ O F F I j | O F F | where ON/OFF represents the set of electric currents corresponding to the gate’s on/off currents, and |ON|/|OFF| is the cardinality. For example, for an AND gate, ON = {I ₁₁} and OFF = {I ₀₀, I ₁₀, I ₀₁}. The noise is defined as M S E = 1 4 [ ∑ i ∈ O N ( I i − I o n ) 2 + ∑ j ∈ O F F ( I j − I off ) 2 ] and the offset as c = I _off . 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 S6; Supplementary Figure S6 we cover the influence of δ > 0. Specific realisations of the four electric current values for a range of different fitness values are shown in Supplementary Figure S4.

For δ = 0 we can exactly rewrite Eq. 12 for the individual gates in form of the quantities introduced in Eq. 11 F AND = − F NAND = 8 3 M l + M r + X M l 2 + M r 2 + X 2 − M l M r − M l X − M r X F OR = − F NOR = 8 3 M l + M r − X M l 2 + M r 2 + X 2 − M l M r + M l X + M r X F XOR = − F XNOR = − 2 X M l 2 + M r 2 (13)

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 S7; Supplementary Figure S7. 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 l 2 > = < M r 2 > = < M 2 > .

Using these approximations, for the AND, OR, NAND, and NOR gates we can write < F AND/OR > ∝ < M > < X 2 > + < M 2 > + Var M ⋅ 1 − corr M l , M r ∝ − < F NAND/NOR > < F AND/OR 2 > ∝ < X 2 > + < M 2 > + Var M ⋅ 1 + corr M l , M r < X 2 > + < M 2 > + Var M ⋅ 1 − corr M l , M r ∝ < F NAND/NOR 2 > Var F AND/OR ∝ < X 2 > + Var M ⋅ 2 + corr M l , M r < X 2 > + < M 2 > + Var M ⋅ 1 − corr M l , M r ∝ Var F NAND/NOR (14) and for the XOR and XNOR gates < F XOR/XNOR > = 0 < F XOR/XNOR 2 > ∝ < X 2 > < M 2 > . (15)

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 Q NDR = 1 2 ⋅ 1 − t a n h < F AND/NAND > Var F AND/NAND / 2 = 1 2 ⋅ 1 − t a n h < M > < X 2 > / 2 + Var M ⋅ 1 + corr M l , M r / 2 ≈ Var (M) ≫ < X 2 > / 2 1 2 ⋅ 1 − t a n h < M > σ M . (16)

In the last step we have neglected the Pearson correlation between M _l and M _r .

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 NDR required to achieve NAND/NOR gates.

A second statistical access to Eq. 16 is based on the fact that exactly for negative M NDR becomes visible in the input-output dependence. Therefore, when sampling many M values with corresponding probability density p(M), the integral ∫ − ∞ 0 d M p ( M ) denotes the probability to find a realization with NDR. As in general p(M) may be complicated, we aimed to find a function Q _NDR, which should recover essential properties of this integral and thus can also serve as a measure for the ability of the system to show NDR. The following conditions should be fulfilled: (1) Q _NDR only depends on the first moment ⟨M⟩ and the standard deviation σ(M), (2) for ⟨M⟩ ≠ 0 and σ(M) → 0 the function Q _NDR approaches one or 0, depending on the sign of ⟨M⟩, (3) for σ(M) → ∞ the function Q _NDR approaches 0.5, (4) for fixed σ(M) the function Q _NDR decreases for increasing ⟨M⟩ and vice versa. Beyond these limiting cases, one may additionally require that for a Gaussian distribution with |⟨M⟩|≪ σ(M) one reproduces the exact result ∫ − ∞ 0 d M p ( M ) ≈ 0.5 ⋅ [ 1 − c ⟨ M ⟩ / σ ( M ) ] with c = 2 / π . In order to avoid the use of transcendent numbers in our definition of Q _NDR we approximate this result by choosing c = 1. Naturally, there is not a unique choice of Q _NDR , fulfilling all conditions, but our choice in Eq 16 is probably one of the most simplest ones.

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 Q NLS = < X 2 > < M 2 > (17) as a general measure for nonlinear separability (NLS), strongly reflecting the size of the nonlinear coupling X between both inputs.

In Supplementary Section S5, we generalize the approach described in this section towards systems of arbitrary numbers of inputs. However, for the sake of simplicity and as this work only covers two-input devices we will stick to the quantities of Eqs 11, 16, 17.