Neuromorphic nanocluster networks: Critical role of the substrate in nano-link formation

Atomic cluster-based networks represent a promising architecture for the realization of neuromorphic computing systems, which may overcome some of the limitations of the current computing paradigm. The formation and breakage of synapses between the clusters are of utmost importance for the functioning of these computing systems. This paper reports the results of molecular dynamics simulations of synapse (bridge) formation at elevated temperatures and thermal breaking processes between 2.8 nanometer-sized Au$_{1415}$ clusters deposited on a carbon substrate, a model system. Crucially, we find that the bridge formation process is driven by the diffusion of gold atoms along the substrate, however small the gap between the clusters themselves. The complementary simulations of the bridge-breaking process reveal the existence of a threshold bias voltage to activate bridge rupture via Joule heating. These results provide an atomistic-level understanding of the fundamental dynamical processes occurring in neuromorphic cluster arrays.


Introduction
Existing computing technologies based on the traditional von Neumann paradigm may now be approaching their fundamental limits and require an increasing amount of energy resources [1][2][3]. Partly driven by this problem, a new interdisciplinary research area called neuromorphic computing has appeared recently [4], being motivated by the incredible capability of the human brain to perform specific tasks, like pattern recognition, in a very energy-efficient way [5][6][7]. Neuromorphic computing aims to mimic the function of a human brain by assembling artificial neurons and synapses [8] and is considered promising to overcome some of the limitations of conventional computers [9][10][11].
A proposed architecture for neuromorphic computing is the assembly of interconnected, synapse-like switching devices [11][12][13] such as memristors, which change their resistance depending on the history of the bias voltage applied to the system [14,15]. Switching behaviour in such devices can be achieved due to the dynamic formation and breakage of nanoscale conductive filaments [16]. Metal cluster films fabricated by the gas-phase deposition technique have been reported to show similar non-ohmic switching behaviour [17][18][19][20][21][22], therefore networks of atomic clusters and small nanoparticles can function as neuromorphic network building blocks [23]. Benefiting from the developing cluster deposition technologies [24][25][26], this approach is simpler and possibly cheaper than other techniques, such as lithographic methods [27][28][29], as employed to create robust intra-and inter-device connections in a deterministic way.
The mechanisms of formation and breakage of the conductive connections between clusters have been discussed in the literature. In connection with experimental observations, notably conductance measurements [18,30], the electric field-induced surface diffusion and evaporation and the van der Waals interaction between metal atoms are considered the main drivers of the formation of atomic-scale connections between the clusters [31][32][33]. After the formation of the connection, if no bias is applied, the width of the connection increases with time following a powerlaw function, as shown in experimental observation by scanning transmission electron microscopy (STEM) and simulation results using the kinetic Monte Carlo approach [34,35]. The exponent's value can be attributed to the combined effects of diffusion caused by surface tension and the viscosity of the material [36,37]. On the other hand, the breakage of the connections is explained by electromigration introduced by the electric current running through the connections. Numerous studies have emphasized the contribution of electromigration to the breakage or resistive switching of metal nanowires (NWs) [38][39][40][41].
An alternative mechanism for the breakage is the Joule heating induced by the current running through the connections [17]. High temperatures up to thousands of kelvin could be induced by small voltage drops on the order of ~100 mV when the size of a point contact in a complex nanoscale junction is smaller than the mean free path of electrons [42]. The thermally induced breaking of NWs has been frequently reported [43][44][45], and relevant computer simulations have been carried out [46,47]. Our prior research has demonstrated that Joule heating can break the nanofilaments connecting the clusters [48]. Strong local heating can be observed in a NW network when current flows [49]. The breaking process of the NW by heat generation can be observed in situ using transmission electron microscopy [50]. Meanwhile, multiple quantum conductance models and related experiments and simulations have been employed to evaluate the current through a filament [51][52][53][54][55]. These studies have established a foundation for evaluating the Joule heat generation in a nanofilament based on the bias voltage applied.
Although various experimental studies and computer simulations treat neuromorphic cluster dynamics on various substrates [56][57][58], previous studies on the mechanism of connection formation and breakage have typically focused merely on the clusters themselves [31][32][33]. The role of the substrate in the formation and breakage of conductive filaments remains unresolved and needs further exploration.
In this paper, the formation and breaking of a nanoscale conductive filament between nanometer-sized gold clusters placed on a carbon substrate are studied by means of molecular dynamics (MD) simulations. Our result shows that the graphite substrate serves as a pathway for the clusters to diffuse along, and the formation and breakage of the connection always happen along the substrate. The characteristic times for bridge formation by diffusion of atoms along the substrate are evaluated at elevated temperatures (T = 500-800 K) for different cluster orientations and distances between the clusters and used to estimate the bridge formation time under experimentally relevant conditions. After the formation of the bridge, the heat generated by the electric current for different bias voltage values is estimated using a model approach for determining the conductance of nanoscopic systems, and the nanofilament breaking process on the substrate under this thermal stimulus has been simulated. The presented results provide an atomistic-level understanding of the fundamental processes involving the conductive connections between deposited metal clusters. Our results highlight the important role of the substrate in the bridge formation and breakage processes, which offers new insight into the fundamental processes in cluster-based neuromorphic computing systems.

Computational methodology
The MD simulations were performed with the MBN Explorer software package [59]. The creation of the systems and analysis of the results was performed with its accompanying multi-task toolkit, MBN Studio [60].
The studied system is a cuboctahedral Au1415 cluster placed on top of a four-layer thick graphite substrate, see Figure 1. The cuboctahedral structural motif was chosen based on the presence of both (100) and (111) facets of comparable total area. Au1415 cuboctahedral nanoparticles (NPs) comprise six rectangular 8x8-atom (100) facets and eight equilateral triangular 6-atom per side (111) facets. This facilitates the comparison between different orientations using the same model as a reference. Carbon substrates (either graphitic or amorphous carbon) are easily available and have been widely used for atomic cluster deposition. Gold has also been a model material in nanocluster science for decades. Therefore, the cluster-support combination considered in this paper is a system that can be easily replicated in future experimental observations.
The system was simulated using periodic boundary conditions so that connecting bridges were formed between the cluster and its neighbouring periodic images. Four different Au1415 geometry arrangements were simulated to demonstrate the impact of different cluster orientations: (i) the (100) facet of the gold cluster was on the graphite layer (labelled as Au (100)); (ii) rotation of the cluster by 45° along the z-axis based on (i) (Au(100)rot); (iii) the (111) facet of the cluster was on the graphite layer (Au(111)); and (iv) rotation of the cluster by 60° along the z-axis based on (iii) (Au(111)rot). These geometrical arrangements are shown in Figure 1 (b-e) respectively. The distance between the cluster and its periodic image was varied by adjusting the simulation box size. The atoms of bottommost layer of the substrate were fixed in place to avoid translational motion of the whole system. The many-body Gupta potential [61,62] was used to model the interaction between gold atoms. The bond-order Brenner potential [63] was used to describe the interaction between covalently bonded carbon atoms within each graphite layer, while a Lennard-Jones potential was employed to account for the van der Waals interaction between the graphite layers and also the Au-C interactions. The parameters for these interatomic potentials are listed in Ref. [57].
In the simulation of bridge formation, the system's temperature was controlled by means of a Langevin thermostat with a damping time of 1 ps; the integration time step was set equal to 2 fs. After the initial energy minimization and relaxation of each system at 300 K, a series of constanttemperature MD simulations were carried out for each geometry and distance between the cluster and its adjacent images within the temperature range from 500 to 800 K. The total duration of each simulation was 10 ns. Importantly, we have not observed any rearrangements in the structure of cuboctahedral clusters during the initial relaxation and equilibration stages.
The initial geometries for the simulations of bridge breaking were taken from the outputs of the bridge formation simulations. Prior to the simulations of bridge breaking, the systems were equilibrated at 300 K for 100 ps using a Langevin thermostat. Gold atoms that appeared in the gap between the neighboring clusters in the initial Au(100) geometry (see Figure 1 (a)) were assigned to the bridge region. The model proposed by Wexler [52,55] was used to estimate the conductance of the bridge and the current running through the system for a specific bias voltage. The heating rate for the system due to the Joule heating mechanism was then evaluated according to the estimated current and conductance. The velocities of gold atoms in the bridge region were rescaled according to the evaluated heating rate. The system then evolved over 2 ps without a thermostat, and the last frame of the simulation was used to re-evaluate the conductance and current for the chosen bias voltage. Multiple subsequent simulations were performed in this way for a total simulation time of 100 ps. The simulations were carried out for the values of bias voltage in the range V = 10 ~ 100 mV applied between the cluster and its neighboring image, which corresponds to the typical cluster sizes and distances between the electrodes in experiments [18,34]. Although the realistic arrangement of clusters in a randomly assembled cluster film is, of course, more complex than the presented model, the chosen cluster geometries, as well as the temperature and bias voltage ranges, are representative of experiments in this research field. This permits us to explore at the atomistic scale the mechanisms which drive the function of such a device.

Bridge formation process
The process of bridge formation between two gold Au1415 clusters deposited on a carbon substrate is illustrated in Figure 2. The process starts with the diffusion of gold atoms from the clusters along the substrate, leading to a narrowing of the gap between the clusters. Then, a thin bridge with the width of one or two atoms is formed so that the clusters become connected. Figure 2(a) shows a schematic illustration of this process, and Figure 2(d) shows the corresponding snapshot of an MD simulation performed at 700 K. The width of the bridge then increases to about half of the cluster's diameter (equal to ~2.8 nm) within the first ~200 ps of the simulation. The height of such a bridge corresponds to one monolayer of gold atoms. The corresponding schematic image and simulation snapshot are shown in Fig. 2(b) and Fig. 2(e), respectively. The increase of the width then becomes slower, and a second layer of atoms appears on the bridge. The system's geometries obtained at the end of 1 ns-long simulations performed at 700 K are shown in Fig. 2(c) and Fig. 2(f). It is notable that in the simulations carried out in this study, the bridge was always formed along the substrate for any initial distance between the clusters above 5 Å and for any cluster orientations. Figure 2(g) shows the initial Au(100)rot cluster orientation, where the shortest distance between atoms of two neighboring clusters is equal to 7 Å, a value smaller than the chosen cutoff distance for Au-Au interactions. Despite this closeness, the atoms still diffuse along the substrate and form a bridge on the substrate, as shown in Figure 2(h).
Due to the limited simulation time (10 ns), the bridge formation was not observed in the simulations performed at room temperature. Therefore, higher temperatures in the range 500-800 K have been considered, and a relation between the bridge formation time (defined as the time needed for the formation of the monolayer bridge shown in Fig. 2(a)) and the thermostat temperature was established. Figure 3 shows the dependence of bridge formation time on the distance between the clusters. The figure shows the simulation results for temperatures above 550 K at which the bridge was formed in every simulation run. Dots indicate the results of each run, and crosses illustrate the average bridge formation time for a specific temperature. The formation time is shown on a logarithmic scale.  The formation of nanofilaments has been shown to be highly affected by surface atomic selfdiffusion [64]. According to surface diffusion theory, the rate at which a single particle (e.g. an atom) diffuses on a surface is given by: where is the diffusion activation energy defined as the value of the potential barrier that should be overcome by an atom to occupy a neighboring position on the surface, is the attempt rate, is the temperature of the system, is the Boltzmann constant, and ∆ is the characteristic time for particle diffusion on the said surface [65][66][67]. The calculation of the diffusion activation energy is a separate research task that goes beyond the scope of this study. However, we note that the diffusion activation energy for gold atoms diffusing on graphite was calculated in several earlier studies by means of density functional theory and molecular dynamics simulations. The values of reported in these studies are in the range 0.05 -0.10 eV [68,69]. An important outcome from Eq. (1) is that ∆ depends exponentially on 1⁄ .
If one thinks of the bridge formation process as an accumulation of the diffusive steps of all the atoms, it is reasonable to consider the bridge formation time to be governed by an exponential relation with 1⁄ . Indeed, our simulations performed at different temperatures result in similar systems' geometries, with the only variable being the bridge formation time. A fit of our data with the exponential function given by Eq. 1 predicts a bridge formation time of ~10 µs at T = 300 K (see Figure 4). Therefore, we believe that surface diffusion can be the dominant mechanism of the bridge formation process for the systems and the temperature range considered in this study.
An alternative to the single adatom diffusion process for the explanation of bridge formation can be the coordinated diffusion of a group of atoms. In Refs. [70] and [71], the growth and coalescence mechanisms of Au NPs of sizes similar to the ones used in this work were found to be the result of "concerted displacements of many atoms".  Figure 4 shows the dependence of bridge formation time on temperature for different cluster orientations. The results of MD simulations (symbols) were fitted for each cluster orientation with an exponential function (lines). The distance between neighbouring clusters was set to 7 Å for all configurations. Interestingly, we did not observe a bridge formation of Au(111) and Au(111)rot orientation (blue triangles and green diamonds, respectively) at temperatures below 750 K over a simulation time of 10 nanoseconds. The predicted formation time at 300 K is ~10 µs for Au(100) orientation, ~1 s for Au(100)rot, and over one year for Au(111) and Au(111)rot orientations. A reason for the significantly faster formation process for the Au(100) orientation may be that the (100) facets contain more atoms compared with the (111) facets of these clusters (64 and 36 atoms, respectively, see Figure 1). Additionally, the (111) facet of gold is a good match to the hexagonal graphitic structure [72], a factor tending to the stabilization of the cluster. In general, the predicted formation times for all orientations at 300 K show a large variation, from ~µs to years, which reflects the experimental complexity of the switching behaviour in randomly assembled cluster films [18,22].
After the bridge formation, the neck width (defined as the width at the middle of the bridge) continues growing as shown in Figure 2. The "neck width -time" relation for Au(100) cluster orientation at different temperatures is shown in Figure 5(a). Eight to ten independent simulations have been carried out at each temperature, and each symbol shows the average time when the bridge reaches a certain width. The data agree well with a power law ~, where is the neck width of the bridge and is time [34,35]. The relation between the power and the temperature is shown in the inset. It follows from the performed simulations that = 0.34 ± 0.03 and it is practically independent of the temperature. A similar analysis was done as a function of the distance between the neighbouring clusters. Figure 5(b) shows the neck width -time relation for different neighbouring distances at 750 K. The power law ~ is also used for fitting the data and the relation between the power and the distance is shown in the inset. In this case, a linear increase of with the distance can be observed. The growth of the neck width has been experimentally observed previously [35]. Monte Carlo simulations have attempted to shed light onto the underlying mechanism [34,35]. It was shown that the neck radius between two coalescing microscale solid particles follows a power law ∝ , where is specific to the physical coalescence process and lies between 1/6 (surface diffusion) [36] and 1/2 (viscous flow) [37]. Our results shown in Figure 5(a) agree well with the power law behaviour, as discussed, and the power = 0.34 ± 0.03 is in the range of these previous simulation results and very close to the in-situ observation result in Ref. [35] where the power ∼ 0.32. Our finding that the exponent is independent of temperature is also compatible with the previous simulation predictions. Note that in Refs. [34,73] the clusters were spherical, and no substrate was considered, so the cross-section of the neck was circular. In the present simulations, the neck has a bow-tie shape with a thickness of one or two atomic layers (see Figure 2(e,f,h)). However, it can be seen that despite the geometry difference caused by the surface interaction between the gold and the substrate, the neck growth mechanism remains unchanged. On the other hand, our result in Figure 6 shows that the mechanism is affected by the initial distance between the clusters, but the value of is still within the range from 1/6 to 1/2. A longer initial distance makes the diffused cluster atoms deviate more from the circular shape when the bridge is formed. is approaching 1/2 for a longer initial distance, which shows that the viscosity flow mechanism becomes dominant when the geometry of the system deviates from the ideal circular shape.

Bridge breaking process
The bridge breaking process under different values of bias voltage, from 10 to 100 mV, at T = 300 K was also studied. The geometry of the bridge with a width of one or two atoms was chosen to be the initial state because the current starts to run through the bridge as soon as the bridge is formed, as described in previous experimental works [17][18][19][20][21][22]. The conductance of the bridge is evaluated using the theoretical model suggested by Wexler [52]: where is the Fermi wavevector (in our case, gold, = 1.204 Å −1 , see Ref. [74]), 0 = 2 2 ℎ ⁄ is the quantum conductance unit, is the radius of the cross section of the bridge (in our case is the van der Waals radius of gold atom 1.44 Å), is the electron mean free path, is the cross-sectional area of the bridge (in our case = 2 ), Γ is a slowly varying function of order ~1 which can be seen as constant in our case, and = ⁄ is the so-called Knudsen number (the ratio of mean free path and the representative physical length scale, in our case, the radius of the cross section of the bridge). Experimental results show a "conductance -surface area" relation between the gold link formed by mechanical stretch [55]. The data was fitted well with Wexler's expression. The fitted values Γ( ) = 0.7 and = 38 Å. The conductance of the bridge was calculated with the information above, and then the current and power were evaluated for different values of bias voltage.
At the bias voltage values above 46 mV, gold atoms in the bridge region will have intense thermally-induced atomic motion and eventually break the bridge at t~25 ps. Local heat is generated due to different heat transfer properties between the cluster and the substrate. The temperature of the cluster reaches ~700 K when the breakage happens. With a bias lower than 30 mV, the atomic movement in the bridge part is not intense enough to break the bridge. The temperature of the cluster keeps increasing, and the width of the bridge increases, similar to the bridge formation process. Figure 6 illustrates the different results in the case of high and low bias. The estimated conductance, current and power of heat generation are shown in Table 1. The conductance is given in the units of quantum conductance 0 = 2 2 ℎ ⁄ .  Our results indicate that, at the applied bias voltage higher than 46 mV, the bridge between the neighbouring clusters breaks during the first 20-30 ps and then reforms again on the time scale of about 100 ps. After the configuration shown in Fig. 6(b) is formed, there is a dynamic interplay between the bridge breaking and formation processes over the time scale of about 10 ps. After that, a stable bridge is formed again, becoming thicker as the simulation time increases. In contrast, the bridge does not break and continues to grow when the bias is lower than 30 mV. A bias between these two values will result in random growth or breaking behaviour. Therefore, our results indicate that a certain bias level is required to break the bridge, but the threshold is obscure. Several experimental reports mentioned the existence of a threshold [17,18] and the complexity of the switching behaviour near the threshold has been reported via conductance measurements [18][19][20]. Our simulation results reflect these features and can act as a reference for future studies.

Conclusions
In summary, aiming to provide an understanding of the formation and breakage of the conductive filamentary bridges between neighbouring clusters in neuromorphic devices, we built a series of model systems and performed MD simulations on the formation and breaking process of the bridges between clusters on a graphite substrate.
In the bridge forming simulations, a high temperature was used to accelerate the process and provide results within the MD simulation time, and the time scale at room temperature was extrapolated from fittings of a series of these higher temperature results. The diffusion of atoms on graphite substrate contributes greatly to the formation mechanism in our simulations. In the bridge breaking simulations, the conductance of the bridge was estimated by the quantum conductance model, and Joule heating under particular bias conditions was calculated. The temperature variation was applied by re-scaling the speed of atoms on the bridge, and the breaking process was revealed by the simulation. Our results indicate the validity of the Joule heating mechanism in the breaking process of the bridge. We also show the existence of a bias threshold to activate switching behaviour and demonstrate the complexity of the switching process in experiments.
These atomistic simulations provide valuable insight into the mechanisms that govern the behaviour of cluster-based neuromorphic systems. Our models took into account the clusters and the support. Parameters of the model have been chosen according to the experimental conditions, making our results representative and a good reference for further studies. New mechanisms, such as diffusion and different formation and breakage process, were revealed by our model. However, since more complicated and asymmetric situations can exist in a real cluster network, more research is needed to present a full picture of the mechanism in the forming and breaking behavior. We hope this work will contribute to a better understanding of nanocluster networks.
Follow-up works must include multiple nanoclusters, of varying orientations, deposited on a series of substrates. Au binds weakly on graphene, with a small binding energy and large equilibrium separations [75]. A future work might investigate substrates that show strong binding with Au. The comparison of the results with the present work would be of great interest for the realization of cluster-based neuromorphic devices.