Ns-Au/ZrO_x bilayer films are deposited on a commercial glass coverslip (Zeus 2 cm × 2 cm x 1 mm) cleaned in an ethanol ultrasonic bath for 15 min. The two gold electrodes have been deposited by thermal evaporation and they are separated by a 4 mm long gap. Electrode’s dimensions are 3 mm × 7 mm x 80 nm. A 30 nm thick gold nanostructured film (ns-Au) is deposited using a Supersonic Cluster Beam Deposition (SCBD) apparatus (Milani and Iannotta, 1999; Borghi et al., 2022), filling the gap between the two electrodes with a strip 1 mm large and 6 mm long. Subsequently, a zirconia (ZrO_x) nanostructured film of 90 nm thickness is deposited with the same aforementioned apparatus, on the top of the gold thin film (Borghi et al., 2018). The SCBD apparatus is equipped with a Pulsed Microplasma Cluster Source (PMCS) that allows the production of neutral clusters in the gas phase, as described in detail in (Wegner et al., 2006). For further details about fabrication method, we refer to our previously published paper (Profumo et al., 2023). In Figure 1, an Environmental Scanning Electron Microscopy (ESEM) image of the nanocomposite is displayed, showing the granular and porous structure of the ns-Au/ZrO_x thin film.

The as-deposited layered ns-Au/ZrO_x films shows ohmic electrical behavior as previously reported for the single layer of ns-Au film (Mirigliano et al., 2019). After the application of constant voltage higher than roughly 16 V, the device switches to the high resistance state ascribed to the partial breaking and reorganization of the underlying ns-Au film, a so-called forming process, well described for ns-Au film in (Nadalini et al., 2023). Figure 2 shows representative snapshots of the electrical response during the forming process of the non-linear electrical behavior of the bilayered samples: this process is related to a joule heating effect at low voltages A), a non-trivial reorganization of the network at higher voltages due to both heating effects and electroforming mechanisms B) and an abrupt increase in sample resistance at a critical voltage C).

After the forming process, the device shows a non-linear and hysteretic behavior reported in Figure 3. I-V curves are obtained applying DC voltages starting from −20 V to +20 V with velocity sweep 1V/s and 2.3V/s. The measurements show that for faster voltage sweeps higher currents are obtained. We observe peculiar electrical behaviors consisting in the opening of the curves at non-zero voltage and non-zero current and a saturation effect almost independent of the velocity of the voltage sweep. As first approximation, we discuss in the next Section 2.4 the saturation effect that is ascribable to the saturation of the PF conduction mechanism at a certain critical applied field.

We verified the stability of the IV curves acquired at the different sweep velocities for ten cycles, as it is reported in Supplementary Material (Supplementary Figure S1). Furthermore, we evaluated the variability of the conduction responses for different devices, nominally identical, i.e., with the same ns-Au and ns- ZrO_x thicknesses. In the Supplementary Material (Supplementary Figure S2) the box plots of the main electrical parameters of eight devices (initial resistance before the forming process, resistance after forming process, current value at 20 V and maximum percentual difference between forward and backward currents at fixed voltage) are reported. We noticed the same hysteretic and non-linear behavior, with the highest variability in the current values, principally due to low reproducibility of the contact area between the bilayered nanogranular structure and the bulk electrodes, which affects the current flowing in the devices.

In situ electrical characterization with temperature spanning from ambient temperature (300 K) to 200 K, has been performed in vacuum (p < 10–4 mbar) on films mounted on a copper cold finger of a helium mechanical cryocooler. I-V measurements are performed with a SPD3303X-E Programmable DC Power Supply and an Agilent 34410A multimeter at 2V/s and 1V/s.

Before describing the memristive-like model developed for our device, we discuss the Poole-Frenkel (PF) effect both in the classical formulation (Frenkel, 1938) and in the modern formulation (Ongaro and Pillonnet, 1989), which incorporates a saturation effect that must be included in a complete model for an emptying-trap process.

The PF effect can be illustrated (see Figure 4A) as the reduction of the potential barrier trapping the electrons to the Coulombic (charged) centers in the bulk material as a consequence of the applied electric field. Consequently, the escape probability of electrons by thermal excitation is increased by the application of the field. In Figure 4A curve a) shows the typical Coulombic potential as a function of the distance from the charged trap. As the electric field is applied, the electron needs smaller energy to escape the attraction of the trap which is reduced by an amount equal to β E . Once the electron gains the necessary ionisation energy by thermal excitation, it is promoted to the conduction band, or the quasi-conduction band in the case of amorphous materials (Hill, 1971). Thus, there exists a critical value of the applied field such that the potential wall is reduced to the ground state of the trap (see curve c) in Figure 4A). In such a case, no barrier is present and the charged donor traps are completely ionised, thus a saturation effect must be taken into account when high fields are involved.

In the classical formulation (Frenkel, 1938), and in its variations available in the literature, the conductivity σ of the bulk material is proportional to: σ ∝ e − q Φ − β E K T ∝ e α E (3) where E is the external applied field, K is the Boltzmann constant, T is the temperature, q Φ is the ionisation potential and β = q 3 / π ϵ ϵ 0 is the PF constant, where q is the charge of the electron and ϵ and ϵ 0 are the dielectric constants of the vacuum and of the amorphous material respectively.

Thus, the current density J is proportional to: J ∝ σ E ∝ E e α E (4)

As firstly recognized in (Ongaro and Pillonnet, 1989), the conductivity, σ , in equation 3 is an indefinitely increasing function of the applied field, which is in obvious contrast with the emptying of the donor traps underlying the PF effect described above. This inconsistency is due to the use of the Boltzmann function as a steadily adopted and convenient approximation of the Fermi–Dirac in the literature. Specifically, Ongaro and Pillonnet demonstrated that the use of the Fermi-Dirac function to describe the population statistics of the donor traps in their model leads to the description of the saturation of the PF effect at high fields. We report here their assumptions (Ongaro and Pillonnet, 1989; Harrell and Gopalakrishnan, 2002) which we consider valid to model our device in the following:

• The acceptor levels are positioned well-below the Fermi level (see Figure 4B) thus they are completely filled for all values of field and temperature.

The use of the Fermi-Dirac permits the description of the saturation and permits to avoid assumptions on the relative position of the energy levels of E d and E f , as in previous models.

Thanks to the assumption that all acceptors are filled, the equation of conservancy of the charge can be written in the form of n d = n + N a , which leads to the expression of the relative density of electrons in the conduction band n r = n / N c : n r = n N c = 1 2 e − η + z E 1 + s c e η − z E ∙ 1 + 4 s c − 1 e η − z E c 1 + s c e η − z E 2 − 1 (6) where η = q Φ / K T and z = β / K T . For the detailed derivation of the above equation, we refer to the original paper (Ongaro and Pillonnet, 1989) and to the more recent paper (Harrell and Gopalakrishnan, 2002). Figure 5 shows the saturation of the relative density of electrons, n_r, as a function of the applied field.

To fit the parameters of our model, we used the experimental data shown in Figure 7B. Simulation results for a complete cycle, −20 to 20 V, are illustrated in Figure 7A. The relevant features of the curves, such as the opening of the hysteresis for both non-zero voltages and currents and the saturation effect, are captured by our model. Figure 7 shows on top-left the current time evolution (red curve) for a triangular voltage signal (green curve) applied. The same dynamics is plotted in the I-V plane in the top-right of the figure. In the second row, the state variable w is plotted in time and against the applied voltage signal. In Figure 7B the experimental data are shown for comparison.

Figure 8 shows the simulated I-V curves for two voltage sweep velocities: 2.3 V/s (blue) and 1.0 V/s (orange). Simulation results are in good agreement with experimental measurements (see the experimental I-V characteristic shown in Figure 2) and (Profumo et al., 2023). Furthermore, a good agreement with parameters of (Chiu et al., 2005) is obtained. Fixing T at ambient temperature we obtained a Φp = 0.7 V (1.1 V in (Chiu et al., 2005)) and, given a characteristic length over which the voltage is applied of 200 nm (which is reasonable due to the confinement effect of the Au clusters and corresponding to a medium electric field of 1 MV/cm), a ϵ r = 6 , 5 ( ϵ r = 6 , 25 in (Chiu et al., 2005)). To obtain all curves shown in this section, an initial non-zero guess for the initial condition of the state variable was given as it is reasonable to assume that the initial DC bias of −20 V affects the initial state of the internal state variable.

Limitations of our modelling are believed to be due to the role of second order state variables, as the temperature. In fact, temperature affects both the mobility of ions and the electronic description, thus its role needs to be investigated for future modelling of this device. Furthermore, our modelling does not take into account spontaneous relaxation of the device.

An analysis from experimental data of the power exponent parameter G = d l n I / d l n V (Acha, 2017) at room temperature gives an insight into the different predominant conduction mechanisms in the two regimes of the hysteretic curve. In Figure 9 we report G curves as function of V 1 2 for the forward A) and backward B) branch of the positive voltage portion of the I-V curve, both for the experimental data (blue dots line) and the model proposed (red lines). The forward branch is well described by the model for high voltages while for low voltage a significant discrepancy arises. Furthermore, the peak in the forward gamma curve is often recognized in the case of PF conduction with a resistor in series and one in parallel (Acha, 2011; 2017; Acha et al., 2016). On the other hand, for the backward branch the model well describes the experimental curve besides for the raising of resistive switching activity (here clearly present in the raw data, note that some data points are out of the figure box). A zero-passing linear dependence (characteristic of pure Schottky interfaces) is also reported. These curves further support the hypothesis of the two mechanisms proposed in the memristive-like model. Nevertheless, they underline a poor accuracy of the model at low voltage which can be justified considering a higher noise in the experimental measurement at low voltages and due to capacitive effects (Zhou et al., 2020).

Remarkably, Poole-Frenkel (P-F) emission involves a process similar to Schottky emission, in which thermally excited electrons emit from traps into the dielectric’s conduction band. P-F emission is thus also known as internal Schottky emission (Zafar et al., 1998) and no marked distinction can be made between the bulk- and electrode-limited conduction mechanisms (Chiu, 2014) especially in such a nanogranular structure (Felix et al., 2011). The slower dynamics of ionic migration (when compared to the electrons’ almost-instantaneous dynamics) regulates the switching between the two conduction mechanisms involved, particularly promoted by the network responsiveness.

The low sensitivity of the in-situ apparatus at low current values prevented us from the evaluation of the gamma parameter at lower temperatures. Nevertheless, we reported in Figure 10A representative IV curves acquired at three decreasing temperatures, and on the right (Figure 10B) the Arrhenius plot of the current at 20V, for the backward and forward voltage ramps. These data confirm a hysteretic behavior of the current decreasing at lower temperatures. Furthermore, an exponential fit of the Arrhenius plot stress the accuracy of the proposed model, given the expected temperature dependence of both the Schottky and the Poole-Frenkel mechanisms. For the backward ramp where we propose a main contribute of the Schottky mechanism, a quadratic term in T is included (note that this term has not significant contribute to the fitted curve). The decrease of the hysteresis with temperature can be interpret also with the decreasing mobility of oxygen vacancies in the film (Aleskandrova et al., 2006).