Brain-inspired computing with fluidic iontronic nanochannels

Significance The brain’s computing principles (neurons connected by synapses) and information carriers (ions in water) both differ fundamentally from those of conventional computers. Building on this distinction, we present an aqueous memristor that emulates the brain’s short-term synaptic plasticity features through ion transport in water, mirroring the natural processes in the brain. This device, which is inspired by and understood through a theoretical model, is applied as a synaptic element for reservoir computing, a brain-inspired machine learning framework. Thus we implement a brain-inspired computing element in a brain-inspired fluidic medium, representing a considerable step toward computing devices that proverbially both walk and talk like the brain.

Neuromorphic computing aims to replicate the information processing of the human brain, which is orders of magnitude more energy efficient than conventional computing devices [1,2].This is paramount as the unsustainable trend of energy consumption by computers is growing at an exponential rate, driving investigations into brain-inspired computing paradigms [3].To pursue brain-like information processing, a device structure that goes beyond the conventional von Neumann architecture is necessary [4].To this end, memristors (memory resistors) have emerged as promising artificial analogues to biological synapses that enable brain-inspired circuit architectures [2,3,[5][6][7].
Despite the successful implementation of memristors in various conventional platforms, the vast majority of these devices consist (at least partially) of solid-state components, rely on only a single information carrier (usually electrons or holes), and only respond to electric driving forces [3,7].These limitations contrast with the brain's nimble synapses, which can utilize both electrical and chemical signals by relying on transport in an aqueous environment of various ionic and molecular species in parallel [8].In light of this disparity, an emerging and exciting approach seeks inspiration not only from the architecture of the brain, but also from its aqueous medium and ionic signal carriers [9].These so-called iontronic devices employ ions moving in an aqueous environment to carry information, offering the promise of multiple information carriers, chemical regulation, and bio-integrability [10].Consequently, various iontronic memristors have been presented [11][12][13][14] that can exhibit synaptic plasticity features [15,16], and utilize chemical regulation [17,18].Additionally, recent advancements have been made in employing iontronic devices for signaling and computing, with theoretical proposals [19][20][21], and demonstrations of traditional truth tables [22][23][24], respectively.Despite these prospects, the development of aqueous neuromorphic devices is still in its infancy and neuromorphic computing implementations remain a challenge [10,25,26].
Here we theoretically propose and experimentally realise the implementation of an aqueous volatile memristor as a synaptic element for iontronic neuromorphic computing.This device is stable over long periods of time, providing reliable and distinct responses to temporal inputs, enabling its use as computing element.Additionally, device fabrication is fast, cost-effective and easy via a soft-lithography process that is almost free-shaping.By constructing a channel of a certain chosen length, made easy by the flexible fabrication process, we can design our channel to feature a specific timescale chosen from a wide range, a desirable property of memristive devices [27].The volatile nature, i.e. decaying conductance memory when driving forces are removed, with adjustable memory retention times makes our memristor a promising candidate for reservoir computing, a brain-inspired machine learning framework which has drawn attention due to its capability of handling complex time series and sequential tasks [28][29][30][31][32][33].We implement benchmark protocols of classifying (handwritten) numbers that are encoded as temporal signals.Our aqueous channels process the time series, distinguishing them for subsequent in silico classification with a simple readout function, performing (at least) comparable with more conventional solid-state platforms employing similar protocols [32][33][34][35].
Our iontronic device is understood through a quantitative theoretical model that directly derives from continuum transport equations and identifies an inhomogeneous ionic space charge density, observed between the colloids [36], as (a general) main ingredient to induce salt concentration polarisation and consequent ion current rectification.Moreover, our theory elucidates how the voltage-driven net salt flux and accu- An action potential triggers neurotransmitter release (not depicted) from the presynaptic neuron (orange), binding to receptors on the postsynaptic neuron (yellow), potentially inducing ion transport and altering its membrane potential [8].The dynamic channel conductance is analogous to the synaptic strength.(e) Current measurements (blue) when four consecutive 5 V pulses and five read pulses (green) are applied.(f) Short-term plasticity features observed in the channel (blue) and predicted by the theory, where we show the full (numerical) solution for g(t)/g 0 (red, dashed) and the measurements this would correspond to in the experiment (red, dots).Four consecutive voltage pulses with ∆t smaller than the channel's memory retention time τ leads to facilitation (top) and depression (middle) for pulses of 5 V and −2.5 V, respectively.The short term characteristic is clearly visible when ∆t ≫ τ, in this instance no cumulative change in conductance is found (bottom).mulation, that underpin the (transient) concentration polarisation, surprisingly combine into a diffusionlike conductance memory timescale that quadratically depends on the channel length.Consequently, the theory correctly predicts the voltage-dependent (dynamic) conductance, thereby facilitating a great acceleration of the experiments by pointing out the relevant signal voltages, signal timescales, and suitable reservoir computing protocol.
The combination of i) a stable fluidic memristor that can be designed to feature a specific memory timescale, made with ii) an almost free-shaping and cost-effective soft-lithography fabrication process, iii) a theoretical model that quantitatively describes and predicts the device dynamics, and iv) the implementation of an aqueous iontronic device as an element for neuromorphic computing, forms a significant advancement towards developing iontronic devices that can facilitate the wealth of communication pathways harnessed by the brain.
1 Fluidic NCNM Memristor: Theory and Experiment Our experimental system, as shown in Fig. 1(a), consists of a tapered microfluidic channel of uniform height H = 5 µm and a width that linearly decreases over its length L = 150 µm from 2R b = 200 µm at the broad base to 2R t = 10 µm at the narrow tip.The channel, which connects two deep reservoirs containing an aqueous 10 mM KCl electrolyte, is filled with a rigid face-centered cubic (fcc) crystal structure at a near-closepacked volume fraction η ≃ 0.74 of charged silica spheres with radius a = 100 nm and approximate surface charge den-sity σ c = −0.01Cm −2 .To form the colloidal structure, the colloids, dispersed in a 70% ethanol solution, are injected into deep channel 1 (as in Fig. 1(a)) of 100 µm height, filling the tapered shallow channel up to the base through capillary action.The fluid halts at the interface between the shallow channel and deep channel 2 due to the Laplace-pressure and consequently evaporates, promoting nanoparticle self-assembly into a close-packed fcc.The device is then dried and prepared for the experiments by filling it with the aqueous electrolyte.The electrolyte in the space between the colloids forms a conducting nanochannel network membrane (NCNM) with pore spaces up to tens of nm in the tetrahedral and octahedral holes of the fcc lattice.This connected porous structure can support an ionic current I driven by a voltage drop V over the length of the channel, defined as tip minus base voltage (i.e.V = −V app ).More system parameters and characteristics are laid out in the SI.Some of us previously showed that these type of colloidfilled microchannels are ionic diodes with excellent ion current rectification (ICR) ratios of up to 55 for the channels on which the devices in this work are based [36], and even up to ∼ 1600, the highest reported value at the time, for a more recent version containing two colloidal structures of opposite surface charge [37].The physical phenomenon that underpins the ICR is the strong voltage-dependent salt concentration polarisation, hypothesised to be induced by an experimentally observed inhomogeneous ionic space charge density between the colloidal particles [36].Calculations on a standard colloidal Wigner-Seitz cell model, presented in detail in the Supplemental Information (SI), show that a small (essentially invisible) variation of order ∼ 1% in the colloidal packing fraction between base and tip of the channel can alter the colloidal surface charge density by several 10% for a fixed zeta potential.By assuming macroscopic electroneutrality [38,39], the ionic space charge in the (thin) electric double layers of the charged colloids varies similarly, thereby providing a natural explanation for the hitherto unexplained ionic charge density profile.In order to theoretically investigate the hypothesis that the inhomogeneous space charge density is responsible for the ICR, we employ standard Poisson-Nernst-Planck (PNP) equations for ionic transport to explain that the inhomogeneous ionic space charge density indeed leads to current rectification.Our theoretical framework, based on an efficient slabaveraging approach in tapered channel geometries [40,41], is described below, while the detailed calculations can be found in the SI.
The PNP equations form an effective theoretical framework to analyse ion transport in charged porous materials [42].However, the complex three-dimensional geometric structure of the NCNM, with features on length scales varying from the colloidal surface-surface distance all the way up to the channel length, introduces intricate numerical challenges for fully spatially resolved solutions of the PNP equations.To simplify, we consider slab-averages, i.e. the average along a cross section [38,40,41,43,44], of the electric potential and the ionic concentrations in the porous structure between the colloids.Although this sacrifices on nanoscale details, it does account for the pinched electric field lines towards the channel tip and for the spatial variation of the ionic charge density.Through this method we reduce the three-dimensional Nernst-Planck equation to a one-dimensional form, providing an expression for the total salt and charge flux through the channel.The divergence of the total salt flux qualitatively shows that the experimentally observed inhomogeneous ionic space charge density forms a source (sink) term of salt, resulting in salt accumulation (depletion) upon a positive (negative) applied voltage V .Quantitatively, a divergence-free steady-state condition on the total salt flux provides a differential equation for the voltage-dependent slab-averaged salt concentration profile, which we solve analytically.By viewing the channel as a series of conductive slabs, with the conductance of each slab proportional to the (now known) voltage-dependent salt concentration, we calculate the steady-state channel conductance g ∞ (V ) = I(V )/V .This describes how an increase (decrease) in salt in the channel at positive (negative) voltages makes the channel more (less) conductive.Our theory thus quantitatively confirms the experimental hypothesis that the ionic space charge distribution results in salt concentration polarisation and hence in current rectification [36].Moreover, leveraging the general analytical nature of our theory, we demonstrate that any inhomogeneous ionic space charge density in generic channels (provided they are well-described by slabaveraged PNP equations) is the key ingredient for a sourcesink term of salt and thus for current rectification, derived in detail in the SI.Therefore we not only provide a mechanistic insight as to how the space charge leads to current rectification in the channel of present interest, but this understanding could also explain current rectification in channels with other sources of space charge densities and with other geometries [23,37].Furthermore, this insight may provide inspiration for future design of devices that exhibit current rectification.
In Fig. 1(b) we plot the predicted steady-state current (red) and the experimentally observed current (blue), revealing a similar current rectification.The experimentally observed ICR ratio of 11 is lower than one of our earlier (slightly more complicated) microchannels [36], however it is sufficient for this work and we believe it will be straightforward to optimise our channel for higher ratios in the future as higher ratios were already achieved in similar channels [36,37].
Up to this point, we treated the system in steady-state.When extending our view to the device dynamics we need to consider the time it takes for ions to accumulate into or deplete out of the channel.Utilizing our aforementioned expression for the total salt flux through the channel, we calculate the net flux γV ′ into the channel upon a small applied voltage V ′ and find (see SI for details) that γ ∝ D/L, with L the channel length and D the ionic diffusion coefficient.The contributions to the net flux solely come from the conductive, i.e. voltage-driven, flux term in the Nernst-Planck equation.The proportionality to D/L is intuitive as the electric field strength in the channels is proportional to 1/L and all flux terms are proportional to the ionic mobilities and hence to D. With our expression for the slab-averaged salt concentration profile we also calculate the total change in salt αV ′ upon applying the small voltage V ′ , and we find α ∝ L. This proportionality to L again is intuitive, as the volume of the channel scales with L. The ratio α/γ between this total change in salt and net flux provides an estimate for the concentration polarisation timescale, given by where ξ ≈ 0.42 is an involved dimensionless number that depends on the ratio of the channel widths R t /R b and on the ratio of the internal space charge density at the tip and the base (full expression in SI).Eq. (1) shows that τ is a diffusionlike time, which is remarkable as no diffusion terms from the Nernst-Planck equation go directly into the derivation of Eq. ( 1).Nevertheless, from the aforementioned intuitive dependencies on D/L and L of the net salt flux and total change in salt, respectively, we surprisingly do retrieve a diffusionlike timescale from a voltage-driven process.Inserting our present system parameters, including D = 1.45 µm 2 ms −1 , yields a memory timescale of τ = 1.62 s.
With the the steady-state conductance g ∞ (V ) and the conductance memory timescale τ, we formulate our theory for the time-dependent channel conductance g(t) resulting from a time-dependent voltage V (t).We intuitively expect g(t) to relax towards the instantaneous static conductance g ∞ (V (t)) on a time scale τ, with a vanishing rate of change dg(t)/dt once g(t) = g ∞ (V (t)).We employ a simple equation of motion (which we more rigorously support in the SI) that characterises this relaxation process and write where Eq. ( 2) is a form that was also successfully applied in previous studies on various memristors [18,20,41,45] and Eq. ( 3) is Ohm's law.
To demonstrate the memristive properties of our device we impose a sinusoidal voltage V (t) of amplitude 5V and frequency f = 0.1 Hz as shown in Fig. 1(c, green).The resulting current-voltage (I − V ) diagram is shown in Fig. 1(c, blue), where a clear pinched hysteresis loop is observed, the hallmark of a memristor [46].This loop features a pronounced memory effect (i.e. more open hysteresis loop) compared to various fluidic memristive devices [12,17,18,47], allowing for a wide range of comparatively high conductances.Moreover, in Fig. 1(c, red) we see that our theory shows good agreement with experimental findings.
Memristors are recognised as artificial analogues to synapses, the connections between neurons [6].In neuronal communication, the change in membrane potential of a postsynaptic neuron, triggered by an influx of ions in response to a signal from a connected presynaptic neuron, is a measure for the synaptic connection strength [8].Similarly, our device's measured ion current, also a result of ion flux through a membrane, draws a parallel between the channel's conductance and synaptic strength in neurons, as schematically depicted in Fig. 1(d).A crucial aspect of neuronal functioning is shortterm plasticity (STP), which allows neurons to adjust their synaptic strength in response to recent input history, being of key importance in information processing [48].STP involves FIG. 2. Stability of our NCNM memristors.Averaged current (blue) measured over 50 subsequent voltage pulse train cycles, taking roughly 30 minutes.Each train consists of four write-pulses, of -2 V and 5 V respectively, interspersed with five read-pulses of 1 V.All cycles showed essentially the same behaviour, producing a spread with standard deviation for each current measurement around the average of maximally ∼ 7% (blue-gray), demonstrating the device stability.changes in the synaptic strength that decay over timescales ranging from milliseconds to minutes, where an increase of the synaptic strength is called (short-term) facilitation and a decrease (short-term) depression [48][49][50].To demonstrate that our fluidic memristor can mimic these aspects of neuronal STP, we apply four consecutive positive and negative "write-pulses" of 5 V and −2.5 V, respectively, with a 0.75 s duration, separated by intervals of ∆t = 0.75 s < τ smaller than the memory retention time τ.The asymmetric 5 V and -2.5 V voltages helped optimize the conductance response.We measure the channel conductance by applying small and short "read-pulses" of -1 V and duration 50 ms after each writepulse.In Fig. 1(e) we show the current measurements (blue) when four 5 V write-pulses and five -1 V read-pulses (green) are applied, with energy consumption of ∼ 1 − 10 µJ for the write-pulses and ∼ 10 − 100 nJ for the read-pulses.The readpulses are converted to the measured channel conductances shown in Fig. 1(f).As illustrated in Fig. 1(f, blue), our fluidic memristor exhibits both facilitation (top graph) and depression (middle graph), hence replicating the characteristic features of neuronal STP.These results are the average of three devices with two measurements per device, each showing quantitatively similar behaviour.In Fig. 1(f, bottom graph) the short-term character of the response is prominently visible when the interval between the pulses is much longer than the typical memory retention time τ and no cumulative change in conductance is observed.Our experimental findings are mostly in good agreement with Eq. (2) shown in Fig. 1(f, red), the only notable discrepancy being that the measured strength of depression is weaker than predicted.The overall agreement emphasizes the robustness and predictive power of the theoretical model in quantitatively characterizing the device properties.
To reliably perform neuromorphic reservoir computing our devices need to repeatedly produce the same response to signals.In Fig. 2 ity of our device by applying 50 subsequent cycles of four write-pulses of -2 V and 5 V respectively, taking roughly 30 minutes to complete.The resulting current, averaged over all 50 cycles (blue), shows a narrow spread with a standard deviation (blue-gray) of at most ∼ 7%.In the SI, we present additional measurements, including 26 cycles over 4 hours for all 16 different four-pulse voltage trains, reliably finding a similar conductance modulation each cycle with conductance standard deviations of typically a few % and no more than ∼ 10%.As fabrication is fast and easy, we normally constructed new devices for new sets of experiments.However, devices were in fact stable enough to be reusable, but did dry out if kept for long with our current way of storing, requiring cleaning the salt residue before using again.The device stability is an important feature that enables it to reliably distinguish a large number of different time series, underpinning the reservoir computing, as we discuss later.
The typical timescale τ over which the conductance memory is retained is an important property of memristive systems and the ability to incorporate a wide range of timescales is desirable [27].As per Eq. ( 1), we predict that the memory timescale of a device can freely be chosen from a wide range of options by constructing it with the appropriate length L or radii R t /R b , here we focus on the L-dependence.The prediction that τ is determined by a diffusionlike time ∝ L 2 /D, despite the channel being voltage driven, forms a specific, non-trivial and easily falsifiable prediction.To test this we fabricated channels of lengths 50 µm, 100 µm and 150 µm and determined for all three at which frequency f max of an applied sinusoidal voltage the area enclosed in the hysteresis loop in the I − V diagram is maximal.Using the relation 2π f max τ = 1, we can find an estimate for the timescale τ [51], and check whether τ ∝ L 2 .The natural relation 2π f max τ = 1 entails that maximal hysteresis is observed when the voltage changes over the typical memory time τ, i.e. providing enough time for the conductance to change, but not enough to reach its steady-state.We mathematically derived this relation for general memristors described by Eqs.[2]- [3] when a sinu-soidal voltage is applied [51].Additionally, since the equilibrium channel conductance is inversely proportional to the channel length g 0 ∝ 1/L we expect the overall current, and thus the overall loop area, to decrease for longer L at the same voltage.In our experiments we indeed find that channels of different lengths respond to different frequencies, and show overall decreasing conductances for increasing L, as can be seen by the various hysteresis loops in Fig. 3(a,b,c).By comparing the different values for f max we not only confirm that f −1 max ∝ L 2 , but that quantitatively we have good agreement with the theory.Thus by manufacturing channels of different lengths, facile via the flexible fabrication process, a wide range of memory timescales can be achieved.Therefore our device offers a versatility important for tasks that require processing of signals over various timescales [27].

Aqueous Reservoir Computing
The short-term memory properties of our fluidic memristor, with a memory retention time that can easily be chosen for each device from a wide range of options as shown in Fig. 3, make it a promising candidate for performing reservoir computing, a brain-inspired framework that leverages a fixed dynamic network, or "reservoir", to transform complex temporal input data into an output that can be easily classified [31].Unlike traditional neural network approaches, where the full network needs to be trained, in reservoir computing only a comparatively simple read-out function (here a single-layer neural network) that classifies the output of the reservoir requires training [28].These properties have driven interest in reservoir computing for analysing a variety of temporal signals [28,30].Generally, a device needs two properties for it to be applicable in reservoir computing, (i) a short-term memory and (ii) nonlinear dynamics [30], which our device satisfies as shown in Fig. 1

(f).
To demonstrate the reservoir computing capabilities of our fluidic memristor, we carry out an established benchmarking The confusion matrix on a test set of 2,000 samples, showing an overall accuracy of 81%, comparable with recent reported results using more conventional platforms [32,33].
protocol of classifying handwritten numbers using reservoir computing [32,33].To build up to this, we first employ a standard method of separating 4-bit strings and show that our memristor already performs remarkably well at this initial task compared to previous results using more conventional platforms [34,35].All 2 4 = 16 combinations are translated into a series of voltage pulses with a duration of 0.75 s, separated by intervals of 0.75 s, where a 0 and a 1 correspond to a voltage of -2.5 V and 5 V, respectively.In theory, applying the voltage pulse trains should yield 16 distinct conductance time traces g(t) shown in Fig. 4(a, top), effectively mapping the 16 possible input patterns onto the 16 different conductance values after the fourth pulse, purely by virtue of the device properties.Excitingly, when we apply this protocol in experiments we find the predicted 16 distinct conductance signatures as shown in Fig. 4(a, bottom), featuring the exact same ordering of conductance values and quantitatively similar changes in conductances, once more highlighting the predictive power of the theory.To obtain Fig. 4(a, bottom), all 16 different voltage pulse trains were repeated twice on three different devices, producing six measurements in total and yielding the average conductances shown in Fig. 4(a, bottom).All 6 runs displayed essentially the same behaviour, producing a spread of conductances after the fourth write-pulse for each pulse train, with standard deviations of around g/g 0 ∼ 0.06 − 0.26 (complete list in SI).Individual device variability is lower, with typically standard deviations of several g/g 0 ∼ 0.01 and no more than 0.16.As we perform an analog computing method, rather than discrete logic, the possible overlap between measured conductances is not an issue and the various time series can reliably encode (handwritten) numbers, as we show next.
To illustrate how the results shown in Fig. 4(a) can be leveraged to classify more complex data inputs with an explanatory example, let us consider the simple single-digit numbers 0-9, represented by black and white 4 × 5 pixel images.By converting a row of 4 pixels to a string of bits by letting a white pixel correspond to a "0" and a black pixel to a "1", we can encode the entire image with 5 strings of 4 bits, as shown in Fig. 4(b) for the number "2" (other digits are shown in the SI).These bit-strings then generate 5 distinct signature outputs, as we saw in Fig. 4(a).A single-layer fully connected 5×10 neural network is then trained in silico to classify the 5 measured conductances as numbers.This protocol is schematically illustrated in Fig. 4(c).Other types of simple readout functions could possibly also suffice.We trained our read-out network in silico using the results shown in Fig. 4(a, bottom).To incorporate the (device-to-device) variability, each individual pulse was subject to some noise newly drawn from a normal distribution with mean 0 and standard deviation given by the experimentally determined standard deviation for that specific voltage train.During training, we repeated this process 100 times for each of the numbers 0-9, achieving perfect classification of all 10 digits with noise-free inference measurements.If we also take the noise into account during inference, we still achieve an overall accuracy of 95%, highlighting the system's robustness against noise.Note that actual training is only performed on a simple and small neural network, that would otherwise not be capable of handling temporal inputs, while the "hard" work of separating the time-dependent signals is handled by the internal physics of our fluidic memristor.Ultimately, this successful classification of simple digit images serves as an explanatory proof-of-concept for the broader application of performing complex time-dependent data analysis tasks.
Building on our previous experiments with classifying simple digit images, we go a step further by classifying handwritten numbers from the well-known MNIST database.This database contains a large dataset of 28 × 28 pixel images of handwritten numbers and has become a standard dataset to test and demonstrate the classification capabilities of machine learning methods [52].We first converted each image into a 20 × 22 pixel black and white image by trimming the edges off and rounding the grayscales to either black or white pixels, as depicted in Fig. 4(d).Each image was then sectioned into pixel rows of four pixels, which can be encoded with a voltage pulse train as outlined before, leading to a total of 110 conductance states per image.The readout function consists of a single-layer fully connected 110 × 10 neural network, trained on a dataset of 20,000 samples.The training incorporated the conductance response and device-to-device noise using the results shown in Fig. 4(a), with the noise taken into account like we did for the classification 4 × 5 pixel digits.The noise hardly had any effect on overall accuracy, as can be seen by the nearly identical decrease of the loss function during training in Fig. 4(e).Via this rudimentary straightforward protocol, we achieved an accuracy of 81% on a test set of 2,000 samples, comparable with earlier reported accuracies of 83% and 85.6% resulting from the same protocol using solid-state memristors [32,33].In Fig. 4(f) the classification is schematically depicted in a confusion matrix, where we see how often each combination of true and predicted numbers occurred in the test set.
Our successful implementation of an aqueous iontronic device as a synaptic element for reservoir computing with a per-formance that is (at least) on par with more traditional platforms [32][33][34][35] is a promising demonstration of the potential that our fluidic platform offers for brain-inspired computing.The functionality extracted from our simple devices is remarkable, obviating the need for complicated circuits to distinguish the various time series of interest here, instead relying on the stable conductance modulation of individual devices.Consequently, as we advance towards circuits integrating multiple coupled NCNM devices, we anticipate that the robust individual device properties demonstrated here will enable the realization of expanded functionalities with relatively simple circuits.

Discussion and Conclusion
We implemented a fluidic iontronic volatile memristor as a synaptic element for neuromorphic reservoir computing, while the device relies on the same aqueous electrolyte medium and ionic signal carriers as biological neurons.Our memristor consists of a tapered microchannel that features a conducting network of nanochannels embedded in a rigid colloidal structure, forming a nanochannel network membrane (NCNM).Device fabrication is fast, cost-effective, and easy via an almost free-shaping soft-lithography process.The trait that underpins the conductance memory effect of the channel is its steady-state diode behaviour, for which NCNM devices have shown excellent performance [36,37], translating into a wide range of achievable conductances.Additionally, our device exhibited stable and reliable (dynamic) conductance modulation, enabling its use as a computing element.Moreover, the quadratic dependence of the memory timescale on the channel length offers a straightforward method to design a channel to feature a specific timescale, a sought after feature for advancing neuromorphic computing capabilities [27].
Our memristor is inspired and supported by a comprehensive theory directly derived from the underlying physical equations of diffusive and electric continuum ion transport.We experimentally quantitatively verified the predictions of our theory on multiple occasions, amongst which the specific and surprising prediction that the memory retention time of the channel depends on the channel diffusion time, despite the channel being constantly voltage-driven.The theory exclusively relies on physical parameters, such as channel dimensions and ion concentrations, and enabled streamlined experimentation by pinpointing the relevant signal timescales, signal voltages, and suitable reservoir computing protocol.Additionally, we identify an inhomogeneous charge density as the key ingredient for iontronic channels to exhibit current rectification (provided they are well-described by slab-averaged PNP equations).Consequently, our theory paves the way for targeted advancements in iontronic circuits and facilitates efficient exploration of their diverse applications.
For future prospects, a next step is the integration of multiple devices, where the flexible fabrication methods do offer a clear path towards circuits that couple multiple channels.Additionally, optimising the device to exhibit strong conductance modulation for lower voltages would be of interest to bring electric potentials found in nature into the scope of possible inputs and reduce the energy consumption for conductance modulation.From a theoretical perspective, the understanding of the (origin of the) inhomogeneous space charge and the surface conductance is still somewhat limited.These contain (physical) parameters that are now partially chosen from a reasonable physical regime to yield good agreement, but do not directly follow from underlying physical equations.We also assume that the inhomogeneous ionic space charge distribution is constant, while it might well be voltage-dependent.Lastly, our theoretical model treats the complex porous structure in terms of slab-averages, thereby possibly missing out on detailed features.These constraints of the theoretical model could account for some of the discrepancies between theory and experiment, which is notable in the steady-state current in Fig. 1(b) and the decrease in conductance in Fig. 1(f).For the purposes of this work our current approach is sufficient, however, a more in-depth study could offer a more profound understanding into the interesting features of the channel.
In conclusion, in order to narrow the gap between the promise of aqueous iontronic neuromorphic computation and its implementation, our work demonstrates the capabilities of a fluidic memristor by employing it as an artificial synapse for carrying out neuromorphic reservoir computing.Temporal signals, in the form of voltage pulse trains, that together represent (handwritten) numbers were distinguished by individual channels for subsequent in silico classification with a simple readout function, demonstrating (at least) comparable performance to more conventional solid-state platforms [32][33][34][35].Additionally, the device is fabricated with a cost-effective easy soft-lithography process.The achieved computing properties are inspired and supported by a quantitative predictive theoretical model of the device dynamics.Consequently, our work establishes a solid foundation, both theoretically and experimentally, for future investigations into fluidic memristive systems and their application in aqueous neuromorphic computing architectures, paving the way for computing systems that more closely resemble the brains fascinating aqueous processes.

Methods
The fabrication of microchannel and formation of the NCNM for the fluidic memristor is similar to previously reported methods [36,37] and is described in the SI in detail.A master for multi-layered channels (target heights are 5 µm for shallow channel and 100 µm for deep) was created using a multi-step UV exposure with negative photoresist (PR, SU-8 2005, 3050, Microchem Co., USA).After surface treatment of the master with (3,3,3-trifluoropropyl)silane (452807, Sigma-Aldrich, USA) for easy separation, Polydimethylsiloxane (PDMS, Sylgard, Dow Corning Korea Ltd., Korea) was poured and cured by heating.The detached PDMS device was bonded with a slide glass.The formation of NCNM was formed by a self-assembly of homogeneous nanoparticles with negative surface charge in the desired shallow channel using Laplace pressure to halt the solvent at the base and evaporation of solvent.A close-packed fcc was formed by the growth of the ordered lattice induced by the evaporation.
[47] W. Brown, M. Kvetny, R. Yang, and G. Wang We consider ionic transport through a triangular tapered channel of uniform height H, length L, and widths 2R t and 2R b at the (smaller) tip and the (larger) base, i.e.R t < R b , as illustrated in Fig. 1(a) of the main text.In our experiments we have L = 150 µm (unless stated otherwise), R t = 5 µm and R b = 100 µm.We introduce the axial coordinate x, with x = 0 at the base and x = L at the tip, and the cartesian width-coordinate y ∈ [0, R(x)], with y = 0 the symmetry axis and y = R(x) the half width of the channel given by R The channel is filled with a rigid close-packed fcc crystal of charged colloidal spheres (in the experiments with packing fraction η ≃ 0.74, radius a = 100 nm, and zeta potential ψ 0 ≈ −39 mV).The channel connects two large and deep aqueous reservoirs containing a 1:1 electrolyte at room temperature with salt bulk concentration ρ b (in the experiments ρ b = 10 mM KCl with a Debye length λ D = 3.1 nm).The ionic transport takes place through the electrolyte that fills the space between the colloids in the colloidal crystal and is driven by a time-dependent voltage V (t).We define the applied voltage V = V (L) − V (0) to be the voltage at the tip minus the voltage at the base, such that conductance enhancement always occurs when V > 0, regardless of which side is grounded.In the main text we have a grounded tip so we apply V app (t) = −V (t) at the base.The transport is described in terms of the Poisson-Nernst-Planck (PNP) equations, that relate the electrostatic potential Ψ(x, y, z,t) and the cationic and anionic concentrations ρ + (x, y, z,t) and ρ − (x, y, z,t), respectively, to the ionic fluxes j ± (x, y, z,t).In the space between the colloidal spheres we write where e denotes the elementary charge, ε = 80.23ε 0 the dielectric constant of water at room temperature T , k B the Boltzmann constant, and D the ionic diffusion coefficient that we take equal for the cations and the anions.The electrostatics is accounted for by the Poisson equation (S1), the conservation of ions by the continuity equation (S2), and the combination of Fickian diffusion and Ohmic conduction by the Nernst-Planck equation (S3).We neglect electro-osmotic fluid flow, which we expect to be relatively small due to the narrow constrictions of the geometry.The system of equations (S1)-( S3) is closed upon imposing blocking boundary conditions on all solid walls, n • j ± = 0, with n the (inward) normal on the walls of the channel and the colloids, together with Gauss' law n • ∇Ψ = −eσ /ε with σ the surface charge density (on the wall and on the colloidal surfaces).We also impose that ρ ± (x, y, z,t) equals the bulk concentration ρ b at the far end of either reservoir, that Ψ(0, y, z,t) = 0 and Ψ(L, y, z,t) = V (t) to account for the applied potentials.
The resulting closed set of PNP equations and boundary conditions can in principle be solved numerically by finite-element methods.The geometry of a colloidal crystal in a tapered channel, however, is computationally challenging as it requires a spatial resolution on the nm length scale of the electric double layer, on the 1-100 nm length scale of the pore structure in between the colloidal particles, and on the 10-100 µm length scale of the channel dimensions.By treating the complex porous structure homogeneous medium, these equations were modified and successfully solved numerically to describe the simplified physics inside the type of channel of interest here [1], however at a considerable computational cost and no analytic insights.Instead of a computationally costly numerical approach, we will derive analytical results straight from the Nernst-Planck equation (S3) to obtain an analytic approximation for the channel dynamics.This will yield a computationally significantly cheaper theoretical model which can be treated analytically to investigate the origin of the ion current rectifying properties of the channel, and to predict features such as its conductance memory properties.

1.2
Slab-averaged electric field, space charge, and salt concentration Our theoretical approach is based upon a methodology that we successfully developed and applied recently to quantitatively explain the steady-state and dynamic conductance properties of simpler channels filled with a homogeneous aqueous electrolyte [2][3][4].Here we show that this methodology can be extended to the nanoporous channel network that characterizes the channel we study in this work.The colloidal structure within the channel forms a (nearly) close-packed face centered cubic (fcc) crystal at a volume fraction η ≃ 0.74 as we saw before.With a colloid radius a = 100 nm, this means that the pores through which ions can be transported have diameters as large as several tens of nm in the octahedral and tetrahedral holes of the fcc-lattice [5,6], which is much larger than the Debye length of λ D ≈ 3.1 nm that characterises the thickness of the electric double layers.In other words, the channels are mostly in the regime of non-overlapping and hence fully developed thin electric double layers, bringing us into the scope of area-averaging techniques [2,4,7,8].Specifically, this justifies the same underlying assumption as in Refs.[2][3][4] that the local and voltage-dependent total salt concentration ρ s ≡ ρ + + ρ − , the total ionic space charge density ρ e = ρ + −ρ − , and the local electric potential Ψ can faithfully be represented by the y−z-slab-averaged functions ρ s (x,V ), ρ e (x), and Ψ(x), respectively, where we recall that the lateral coordinate x ∈ [0, L] runs from base to tip.Here we explicitly denote the dependence of ρ s on V , while refraining from denoting the explicit (linear) dependence of Ψ on V below for notational convenience.While we expect also a V -dependence of ρ e in a full calculation, we will restrict ourselves to a V -independent form below.
If the slab-averaged electric field lines cannot escape the tapered channel, a realistic assumption in our experiments as the dielectric constant of water is much higher than that of the wall-material, then the slab-averaged electric field component −∂ x Ψ(x) must be proportional to 1/R(x) on the basis of charge neutrality on the length scale beyond the Debye length.Since we define V = V (L) −V (0), i.e. the tip minus the base voltage, the applied voltage also satisfies L 0 ∂ x Ψ(x)dx = V .Combining the scaling with this property, we find with ∆R = R b − R t that This slab-averaged electric field in the channel is therefore proportional to the applied field V /L and gets progressively stronger closer to the tip.The total steady-state salt flux j s (x, y, z) = j + (x, y, z,t) + j − (x, y, z) can now be integrated over slabs in the y and z direction to obtain for the x-component of the total salt flux J Here we introduced the porosity ε fcc = 1 − η to take the volume into account that is excluded to the electrolyte by the colloidal fcc crystal, where we assume the colloidal particles to be impenetrable to the electrolyte.As we take slab averages one would expect an area term, instead of the porosity.However, microscopically, the available electrolyte area through which the ions can diffuse in the slab has a periodicity in x that is dictated by the lattice spacing, which is much smaller than the channel length L and can hence be ignored in our slab-averaged description.Therefore we consider each slab to have the same available surface area for ions, which in this simplified one-dimensional view is the porosity ε fcc [9].Eq. (S5) depends on the slab-averaged ionic space charge density ρ e (x), that picks up contributions from the electric double layers (EDLs) around the charged colloidal spheres.Given that EDLs have spatial extensions as small as the Debye length (here λ D = 3.1 nm) around the colloidal spheres (here of radius a = 100 nm), one might expect the slab-averaged space charge density to be a constant on the much larger length scale of the channel, with a magnitude ∝ ηZ with Z the charge of a colloid.Interestingly, however, earlier measurements presented in Ref. [1] on channels nearly identical to the ones we study here convincingly showed a heterogeneous rather than a homogeneous space charge density that could well be fitted by the monotonic functional form where ρ e,t and ρ e,b are the space charge density at the tip and the base, respectively, which may differ from each other even at zero applied voltage [1].In this work we have ρ e,t = 8.85 mM and ρ e,b = 8.00 mM for the steady-state results of Fig. 1(b) and ρ e,b = 6.3 mM for the other (dynamic) results, the choice for these values is further discussed in Sec. 2. Whereas we take the functional form of Eq.(S6) as experimental input from now on, the reason for this ionic charge heterogeneity on length scales of the channel length remains an interesting open question.Our hypothesis, that we underpin with standard Poisson-Boltzmann calculations of a single colloid in a spherical Wigner-Seitz cell in Sec. 1 .5,shows that a reduction of the colloidal packing fraction from 0.74 to 0.73 can already increase the colloidal charge density (and hence the average ionic charge in the EDLs) by as much as ∼ 20% at a fixed zeta potential.Thus, a small heterogeneity of the colloidal packing fraction from tip to base could provide a microscopic explanation for the spatial dependence of ρ e (x) by assuming macroscopic electroneutrality [7,9].Given the device fabrication (see Sec. 3) such a small asymmetry of η between tip and base is quite well possible although such a small deviation in packing fraction is difficult to observe experimentally.With Eqs.(S4) and (S6) every component of Eq. ( S5) is known except for ρ s (x,V ).We can use Eq.(S5) for J x (x) to find the steady-state salt concentration ρ s (x,V ) explicitly by imposing the steady-state condition ∂ x J x (x) = 0.This yields a differential equation for ρ s (x,V ), that can be solved analytically after inserting Eqs.(S4) and (S7), yielding We note that H, D, and ε fcc do not appear in Eq. (S7) since these drop out of the underlying differential equation ∂ x J x (x) = 0.More importantly, ρ s (x,V ) is voltage-dependent, resulting in the voltage-dependent channel conductance we derive next.

Static channel conductance
The tapered microchannels of our interest are well known to exhibit ionic rectification properties characterised by a static conductance g ∞ ≡ I(V )/V that is a nontrivial function of the applied static voltage V .By viewing the slabs of thickness dx at x ∈ [0, L] as a series of resistors with resistivities ∝ dx/ρ s (x), one can write the static conductance of the channels of present interest as where we made an approximation compared to the more accurate dependence on L/ L 0 (ρ s (x,V )) −1 dx, which reduces computational complexity and yields mostly the same results [2][3][4].The V -dependence stems from the salt concentration dependence on V as given by Eq.(S7).As in Ref. [4], we replace ρ s (x,V ) by max [0.2ρ b , ρ s (x,V )] in the actual (numerical) evaluations of Eq. (S8) in order to account for the possibility of unphysical negative concentrations that could follow from Eq.(S7) at strongly negative voltages in part of the density profiles.(In this regime the underlying assumption that the Debye length is much larger than than the channel dimensions breaks down).
The reference (zero-voltage) conductance g 0 of the channel can, in direct analogy with recent results from Refs.[2][3][4], be written as which includes the volumetric contribution ∝ ε fcc that depends on the channel-geometry parameters calculated with the total charge flux Additionally we also consider a surface contribution ∝ λ D /R pore , with R pore the effective radius of the pores embedded in the fcc crystal, which accounts for the excess conductivity due to the excess salt concentration in the colloidal EDLs [10] and which we determine below.This surface term is of direct relevance here due to the large internal surface in the channel as a result of the colloidal structure.However, the internal structure of the fcc crystal is complex with pores of varying sizes and shapes, and regions with fully developed EDLs in the pores with the size of several tens Treat electric field, charge density, and salt concentration as slab averaged Reduce Nernst-Planck Eq. to a one-dimensional form via integration Ψ(x,y,z) ρ s (x,y,z,V) ρ e (x,y,z)

Net flux
Total salt change Concentration polarisation timescale  Fig. S1.Schematic depiction of the steps taken in our theoretical analysis, involving 1) the slab-averaging assumption, 2) the integration of the salt flux over each slab to reduce the full three-dimensional Nernst-Planck equation to its one-dimensional slab-averaged form as in Eq. (S5), and 3) the solution of the steady-state condition of a divergence-free flux to obtain Eq. (S7).Lastly, 4) the combination of the salt concentration and the salt flux to extract a timescale for the salt concentration polarisation as detailed in Sec. 1 .4.
of λ D neighboured by regions where the EDLs are not fully developed since the shortest distance between colloids is smaller than a few λ D .Therefore, the surface conductance of the channel is a highly non-trivial property and investigating this in full detail falls outside the scope of this work.Instead we will employ an effective method where we treat the internal pores as a channel of radius R pore , with a total slab surface area 2ε fcc R(x)H.This approach then yields Eq. (S9) [10].We find good agreement with the experiments if we take R pore = 2λ D = 0.06a.Although this seems somewhat small compared to the radii of up to ∼ 10λ D for the larger pores, it is still in a reasonable regime given the fact that many regions in a close-packed fcc crystal feature much smaller distances between colloids.For the parameters used in the experiment (see Sec. 2) we have that g s /g 0 ≈ 0.38, so the conductance is still mostly dictated by bulk conductance, as expected with the relatively large pores compared to the Debye length.
The combination of Eq. (S8) and Eq.(S9) provides us with the steady-state conductance function g ∞ (V ) used in the main text.

Space charge density inhomogeneity is a salt source or sink term
The methodology we employ in Sec. 1 .2 to derive salt concentration polarisation can provide insights into designing new devices that also exhibit salt accumulation and depletion.As shown in Sec. 1 .3, the presence of salt concentration polarisation affects the channel conductance and will hence lead to current rectification.Our theoretical model explains current rectification in the colloid-filled channels we present here, but other types of channels also exhibit current rectification as a result of an inhomogeneous space charge [11,12].Our theoretical description might provide insights into such other systems as well, provided some assumptions are satisfied.(i) The underlying PNPS equations only describe continuum transport, so sub-nm or atomic scale systems [13,14] will fall outside the scope of this theoretical framework.(ii) The radial dependence of the potential and the concentrations must be relatively weak and/or short ranged, such that the slab-averaged salt-concentration ρ s (x,V ) and electric field −∂ x Ψ(x) are fair approximations.For instance, a single channel with strongly overlapping electric double layers throughout the entire channel will exhibit a significant salt concentration profile in the radial direction of the channel, possibly resulting in physical phenomena that a slab average ρ s (x,V ) does not take into account.
Eq. [S5] can be generalised to with Q e (x)dx the total ionic charge in a slab of infinitesimal thickness dx with volume A(x)dx at location x, and where we define ρ e (x) ≡ Q e (x)/A(x) the slab-averaged ionic charge density.Eq. (S10) describes a channel with cross-sectional area A(x) and the average area per slab available for ions a (which would be 1 for materials without any solid obstructions for the ions), where we remark that in the simplified one-dimensional view of treating each slab to feature the same available area one can use the porosity as we did in Eq. (S5) [9].In the present channel geometry we have A(x) = 2HR(x), but Eq. (S10) could also apply to hourglass shaped channels [12], T-shaped channels [11], or conical channels [2,4].If macroscopic charge neutrality is ensured within the channel, e.g. in the case of thin electric double layers, and if the electric field lines cannot leave the channel because of a weakly polarising outside medium, then the relation −∂ x Ψ(x) ∝ 1/A(x) must hold.A fully analytical solution as in Eq. (S4) might not always be possible, however the inverse proportionality with the channel area A(x) suffices for the qualitative mechanistic understanding discussed here.Consider, namely, a channel without any applied voltage such that the salt concentration ρ s (x,V ) is constant within the channel.Upon applying a voltage, the electric field −∂ x Ψ(x) will form quasi-instantaneously, so a short time after the voltage is applied we have a fully formed electric field, but still a constant salt concentration ρ s (x,V ).In this case the diffusion term in Eq. (S10) vanishes and due to the aforementioned proportionality −∂ x Ψ(x) ∝ 1/A(x), the only x-dependence that remains in Eq. ( S10) is ρ e (x).Therefore, if we take the divergence ∂ x J x (x) of the total salt flux in Eq. (S10) and apply the continuity Eq. [S2] we find This shows how any inhomogeneous ionic space charge density forms a source or sink for salt term upon applying a voltage, thereby inducing salt concentration polarisation and consequently current rectification.An important understanding is that the ionic charge density must be inhomogeneous, and not just the total ionic charge in the slab.The reason is that even though the total charge Q e (x) ≡ ρ e (x)A(x) in the slab for a constant space charge density ρ e (x) = ρ e could still be x-dependent due to its scaling with A(x), this dependence would cancel out with the 1/A(x) dependence of the electric field.The insight that an inhomogeneous charge density forms a source-sink term previously explained how a constant surface charge density σ in a conical geometry could induce current rectification [2] as the total (surface) charge is in this case given by Q e (x) = 2πR(x)σ dx.Since in this geometry A(x) = πR(x) 2 we see that Q e (x) ∝ A(x) and therefore ρ e (x) ≡ Q e (x)/A(x) ∝ 1/ A(x), exhibiting the required x-dependence.Additionally Eq. (S11) shows why merely a geometric inhomogeneity, such as a tapered geometry, is not enough to induce current rectification; it must go coupled with a spatially varying slab-averaged ionic charge density.The insight of Eq. (S11) could not only explain current rectification in channels with a space charge density step-function as in recent polyelectrolyte channels [12] and NCNM channels with colloids of opposing charge [11], it may also provide specific guidance to design current rectification properties in future iontronics.

Typical conductance memory retention time
As detailed in Refs.[3,4], the process of ion accumulation and depletion is not instantaneous.To investigate this timescale for the channel of interest here we can apply the same approach as in Refs.[3,4].We consider two quantities, the total number of ions N = L 0 2R(x)Hε fcc ρ s (x,V )dx in the channel and the net salt flux J x (0) − J x (L) into the channel.The change of N given by Eq. (S7) upon a small voltage perturbation where α > 0 for our parameters, in agreement with the enhanced (reduced) conductance of a positive (negative) potential V ′ found in the experiments and as can be seen in Fig. 1(b).
At V = 0 the concentration profile is at equilibrium, so for a small voltage perturbation V ′ we can assume ρs (x) = 2ρ b .With this assumption the first term in Eq. (S5) vanishes.The net salt flux into the channel, J x (0) − J x (L), is then determined by the remaining conductive terms where γ > 0 for our parameter choices, again in agreement with the enhanced (reduced) conductance of a positive (negative) potential V ′ found in the experiments and as can be seen in Fig. 1(b).We see that the net flux is determined by the difference in space charge density between the tip and the base (ρ e,t − ρ e,b ).This is consistent with Eq. (S5), where we saw that the total salt flux is proportional to the charge density in the absence of a diffusion term.Therefore, this reconfirms the underlying mechanistic insight that a inhomogeneous charge density is what drives a net salt flux and consequent salt concentration polarisation, as also argued in Sec. 1 .3.1.
The typical time it takes for ion depletion or accumulation, and thus the typical memory retention timescale, is then approximated by τ = α/γ.This yields From the above calculation it has now become clear that τ ∝ L 2 /(4D), which corresponds to the time it takes to diffuse over distances of the order of the channel length L. The remaining involved terms form a dimensionless number For our standard parameters we find ξ ≈ 0.42.This simplifies our notation considerably such that we arrive at the final concise expression Eq. (S14) for τ which can be found in the main text as Eq.[1], For our standard parameter set as laid out in Sec. 2, we have τ ≈ 1.62 s.
To then arrive at Eq. [2] of the main text, we make the natural assumption that the time-derivative of the dynamic conductance ∂ t g(t) depends on the difference with the corresponding steady-state conductance, i.e. ∂ t g(t) = f g ∞ (V (t)) − g(t) for some function f , where f (0) must vanish based on stability arguments.By expanding f up to first order we find ∂ t g(t) ∝ g ∞ (V (t)) − g(t) with the proportionality constant naturally given by the typical timescale τ of the underlying salt concentration polarisation that drives the conductance change, given by Eq. (S14) (Eq.[1] in the main text).

Cell model
In Sec. 1 .2we showed how the inhomogeneous charge distribution leads to salt concentration polarisation and thereby to ion current rectification.However, Eq. (S6) presented in Ref. [1], which describes the inhomogeneous space charge, follows from empirical measurements without a microscopic explanation of the actual origin of this important feature of the channel.Here we will leverage the well-established Poisson-Boltzmann cell-model [15,16] to propose a tentative explanation of the emergence of the observed inhomogeneous space charge.
We consider a dispersion of charged colloidal spheres of radius a at packing fraction η in a 1:1 electrolyte of Debye length λ D , such that the volume per particle equals (4π/3)b 3 with b ≡ aη −1/3 > a.The environment of each particle fluctuates due to colloidal Brownian motion, and hence the calculation of the profile of the electric potential ψ(r; {R}) is a complicated many-body problem that depends on the instantaneous configuration {R} of the colloidal particles.This problem can be reduced tremendously, however, if we assume each sphere to be at the center of an electrically neutral and spherically symmetric Wigner-Seitz cell of radius b.The dimensionless electrostatic potential φ (r) ≡ eψ(r)/k B T is then the same in each cell, and can for r ∈ [a, b] be described by the Poisson-Boltzmann equation with boundary conditions where a prime denotes a radial derivative and where φ 0 is the fixed dimensionless zeta potential in units of k B T /e ≃ 25 mV.This is a closed system of equations that can easily be solved numerically for fixed a, b, λ D , and ψ 0 or rather for fixed dimensionless η, a/λ D and φ 0 .The resulting surface charge density eσ of the spheres is then from Gauss's law given by σ = −φ ′ (a)/4πλ B with λ B = e 2 /4πεk B T = 0.72 nm the Bjerrum length of water.
For dispersions that are dilute enough that b − a ≫ λ D , the EDLs described by Eqs.(S15) are fully developed.In this regime it is well known that the colloidal particles obtain their maximum (absolute) surface charge [17].At high packing fraction where b−a ≃ λ D , the colloidal surfaces discharge and loose their charge completely in the (unphysical) limit b = a.The situation for the colloidal dispersion of our system is intricate, since for a = 100 nm and η = 0.74 we have b/a ≃ 1.105 and hence b − a ≃ 3λ D .In this regime the assumption of spherical symmetry of the EDLs is highly questionable, because the surface-surface distances between a central particle in a close-packed fcc crystal at η = 0.74 vary between smaller than λ D in the 12 directions of its nearest neighbours to about 10λ D in the directions of the 8 tetrahedral and the 6 octahedral holes of the fcc crystal.We can therefore expect that the particles only significantly discharge in the vicinity of nearest-neigbour contact.Within the sphericalcell approximation, we mimic this anisotropy of the colloidal surface charge by averaging over these 12 + 8 + 6 = 26 directions, with equal weight for simplicity, to define the effective surface charge density σ * = (12σ 2 + 14σ ∞ )/26, with the surface charges σ 2 and σ ∞ of the spherical cell with b half the nearest-neighbour distance and b ≫ a the dilute limit, respectively.In Fig. S3 we show the dependence of σ * /σ ∞ on the packing fraction η for the present system parameter with φ 0 = −1.53,which is such that σ C = σ ∞ = 0.01 C/m 2 .It shows that decreasing the packing fraction from η = 0.74 to 0.73, which corresponds to reducing b by as little as 0.5 nm (i.e. by ≃ λ D /6), the colloidal charge can increase by as much as ∼ 20%.A reduction to η ≃ 0.67 and η ≃ 0.63 (which corresponds to b − a = λ D and 2λ D ) increases the surface charge by about 80% and 100%, respectively, compared to that at η = 0.74.On the basis of these simple estimates, we argue that the slab-averaged space charge in the channel (that compensates the colloidal surface charge) can therefore also exhibit the same relative change between tip and base.The variation we use of ρ e,b /ρ e,t = 0.90 for the steady-state results of Fig. 1(b) and ρ e,b /ρ e,t = 0.71 for the other (dynamic) results falls well within this reasonable range of change predicted by the explanation we offer here.We leave an additional estimate on an alleged spatial variation of φ 0 in the channel to explain the heterogeneous space charge as future work.

System parameters
The channel has a base radius of R b = 100µm, a tip radius of R t = 5µm and a height of H = 5µm.The channel connects two reservoirs with a bulk concentration ρ b = 10 mM of aqueous KCl electrolyte, yielding a Debye length of λ D ≈ 3.1 nm.The colloids have a radius of a = 100 nm and carry a uniform charge density of σ C = −0.01C/m 2 , squeezed together during the device fabrication to form a face centered cubic crystal with a close-packed packing fraction η = 1 − ε fcc ≈ 0.74, where ε fcc is the porosity.The maximum space charge at the tip is determined in the same way as in Ref. [1], i.e. ρ e (L) = ρ e,t = 4πa 2 σ C n/(eε fcc Ω) = 8.85 mM with n = 3 4 ηΩ/(πa 3 ) the number of colloids in the channel and Ω the overall channel volume.The space charge concentration at the base is assumed to be ρ e (0) = ρ e,b = 8 mM in the steady-state calculations for Fig. 1(b) and ρ e,b = 6.3 mM in the rest of the manuscript for the time-dependent calculations.These values were chosen such that (i) the change in space charge density does not exceed the reasonable regime predicted by our tentative explanation for the the inhomogeneous space charge in Sec. 1 .5,i.e. the difference between ρ e,b and ρ e,b is not more than ∼ 40%, and (ii) a good agreement is found with the experiments.Although the space charge density ρ e (x) of Eq. (S6) has a clear physical meaning, its precise form and parameters values are not entirely clear and may contain a voltage-dependence that we ignore in the present study and leave for future investigations.We stress, however, that our present overall results do not depend strongly on the detailed value of ρ e,b and ρ e,t , as it only marginally changes the conductance properties in the parameter regime we consider here and the overall behaviour relevant to reservoir computing remains.However, we do note that our theory predicts that some difference between the space charge density at the tip and at the base is necessary for any current rectification, so even though the precise values of ρ e,b and ρ e,t are not of major importance for our overall results, our theory still predicts it is crucial that there is some difference between the two.Lastly, the effective diffusion coefficient is also taken to be similar to Ref. [1] with D eff = ε fcc D = 0.38 µm 2 ms −1 .Fig. S4 shows schematic diagrams describing the fabrication procedures for microfluidic devices using soft lithography: A negative photoresist (PR) was applied to the silicon wafer (SU-8 2005; Microchem Co., Westborough, Massachusetts, USA) using a spin coater and then soft-baked.The PR was then exposed to UV to create a shallow channel and the wafer was hard baked.The unexposed PR was removed to form the shallow 5 µm channel.A deep 100 µm channel was then formed using a negative photoresist, SU8-3050 (Microchem Co., Westborough, Massachusetts, USA), in the same process as above.After completion of the master mold, the surface was treated with (3,3,3-trifluoropropyl)silane (452807; Sigma-Aldrich, St. Louis, Missouri, USA).Polydimethylsiloxane (PDMS; Sylgard; Dow Corning Korea Ltd., Gwangju-si, Gyeonggi-do, Republic of Korea) was then poured over the master mold and heated on a hot plate at 95 • C for 1 hour.The reservoir of the PDMS device was punched out with a 1.5 mm medical punch.The surface of the PDMS and the slide glass were treated using a plasma device (Cute-MP; Femto Science, Hwaseong-si, Gyeonggi-do, Republic of Korea) and bonded together.As shown in Fig. S5, the diluted nanoparticles (with a carboxylic (-COOH) end groups on the surface and with radius a = 100 nm) dispersed in a 2 µL 70% ethanol solution were injected into the deep channel, such that the solution including the nanoparticles fills the shallow channel through the capillary.Due to the neck pressure at the interface between the shallow and deep channels, the flow of the particle dispersion stopped, after which evaporation of the solvent (70% of the ethanol) was induced through the deep channel.This way additional nanoparticles were transported toward the neck by the convective flow that compensates for the loss of solvent by evaporation.This influx of particles promotes the growth (and the compression) of the ordered fcc lattice in the shallow channel.When the self-assembled nanoparticles completely filled the shallow channel with a wet close-packed fcc lattice, the remaining diluted nanoparticles were removed from the deep channel by suction.The finished device was dried for one day and then used for the experiment.

Voltage pulse measurements
To obtain Fig. 4(a) we applied all 2 4 = 16 different voltage trains of four pulses in sequence with 30 s between each pulse train.This was performed on three different devices with two cycles per device, resulting in 6 independent measurements per voltage train, all individually shown in Fig. S6.The average of these measurement is shown in Fig. 4(a), the standard deviations calculated with the six conductance values after the fourth (last) pulse for each pulse train are shown in Table I.These values are used in the main text to take the measured (device-to-device) variability into account.Therefore we stress that all (device-todevice) variability visible in Fig. S6 is incorporated in the main text.
To ensure that our devices also remain reliable and stable over longer periods of time, we performed various additional cycles of pulse trains over a single device.In Fig. S7 we show a 50 cycle repeat, which lasted about 30 minutes, of the pulse train corresponding to the bitstring 0101.We see a remarkable stability, with essentially the same current response each cycle (which we also show as Fig. 2 in the main text), yielding a narrow spread for each pulse and reliably measured altered conductances with conductance standard deviations of g/g 0 ∼ 0.03 − 0.05.Moreover, we cycled through all 16 different bitstrings for a total of 26 cycles, lasting around 4 hours, and found essentially the same current response for each cycle, even after the device had been cycling for 4 hours.The resulting average normalized conductances per pulse, per cycle, are shown in Fig. S8, with the corresponding standard deviations in the range g/g 0 ∼ 0.02 − 0.15.
We repeated a similar protocol, but now with the ground at the base and the applied voltage at the channel tip.In this instance, a "0" corresponds to a pulse of 2 V, while a "1" corresponds to a pulse of -5 V. Pulse duration and interval are still 0.75 s, the read pulses are 1 V of 50 ms duration.The result is shown in Fig. S9, showing a similar clear separation of the different bit-strings.To again ensure the device stability, the bit-strings 1111, 0101 and 0011 were repeated 5 times, which we present in Fig. S10.Here we show the 5 individual measurements (light grey), the average of the measured current (black) and the calculated normalized conductances in the bar plots, averaged over the 5 measurements.The error bars depict the measured standard deviations, where we find good reproducibility for each voltage pulse train.
In Fig. S11 we schematically show how we use the read pulses to calculate the channel conductance.We calculate the difference in (average) current during a read pulse and the measured current just before the read pulse to calculate the conductance.
In the main text we classified simple single-digit images consisting of 4 × 5 black and white pixels.In Fig. 3(b) we only showed the "2" as an example, the other digits used are depicted in Fig. S12.

Timescale measurements
In Fig. 2   Each column of two figures corresponds to a device.The average for each different pulse train is used to create Fig. 4(a).Standard deviations for the conductance measurement after the fourth (last) pulse were also calculated with these six measurements and shown in Table I

FIG. 1 .
FIG. 1. Features and properties of our iontronic memristor through theory and experiment.(a) Schematic (left) and pictures (right) of the device.The channel connects two reservoirs of aqueous KCl electrolyte and incorporates a rigid colloidal structure, forming a network of nanochannels between the colloids.(b) Steady-state I − V curve observed in experiments (blue) and predicted by our theory (red), showing a similar current rectification property.(c) Dynamic I − V curve in response to a sinusoidal voltage over the channel (top, green).The theory (bottom, red) and the experiments (bottom, blue) both exhibit a similar pinched hysteresis loop.(d) Simplified schematic of synaptic signal transmission.An action potential triggers neurotransmitter release (not depicted) from the presynaptic neuron (orange), binding to receptors on the postsynaptic neuron (yellow), potentially inducing ion transport and altering its membrane potential[8].The dynamic channel conductance is analogous to the synaptic strength.(e) Current measurements (blue) when four consecutive 5 V pulses and five read pulses (green) are applied.(f) Short-term plasticity features observed in the channel (blue) and predicted by the theory, where we show the full (numerical) solution for g(t)/g 0 (red, dashed) and the measurements this would correspond to in the experiment (red, dots).Four consecutive voltage pulses with ∆t smaller than the channel's memory retention time τ leads to facilitation (top) and depression (middle) for pulses of 5 V and −2.5 V, respectively.The short term characteristic is clearly visible when ∆t ≫ τ, in this instance no cumulative change in conductance is found (bottom).

FIG. 4 .
FIG. 4. Our iontronic memristor as a reservoir computing element.(a) Theoretical prediction (top) and experimental observation (bottom) of the relative change in conductance as a response to the 2 4 = 16 different possible bit-strings, where a "0" and "1" correspond to pulses of −2.5 V and 5 V, respectively.Three separate devices were used to find the average conductance and typical variation in response to each unique voltage train.(b) Depiction of how a "2" can be transformed into 5 distinct bit-strings (other digits are depicted in the SI).(c) Schematic of how the number "2" from (b) is translated to 5 voltage trains, yielding 5 conductance values after the fourth (last) pulse.The conductances and variabilities from (a) were then used to train a single-layer fully connected 5 × 10 neural network in silico, that converts the conductances to a classification of the number "2".(d) Nine examples of handwritten numbers from the MNIST database [52], where the (originally 28 × 28 pixel) images are trimmed to 20 × 22 pixel images, the grayscales are rounded to either white or black pixels and then the image is segmented into 110 bit-strings.(e) The loss function (mean squared loss) during training per training round when the experimentally found noise, which is experimentally quantified using the (device-to-device) variabilities found in our result in (a), of the devices is not taken into account (orange) and when it is taken into account (blue).(f) The confusion matrix on a test set of 2,000 samples, showing an overall accuracy of 81%, comparable with recent reported results using more conventional platforms[32,33].

Fig
Fig. S2.The memory time proportionality function ξ as a function the internal space charge ratio ρ e,b /ρ e,t in (a) and the tip and base radii ratio R t /R b in (b), for three different values of R t /R b and ρ e,b /ρ e,t , respectively.In both (a) and (b) the blue graph represents R t /R b = 0.05 and ρ e,b /ρ e,t = 0.71, respectively, the system parameters used in the main text, except for Fig. 1(b) which uses ρ e,b /ρ e,t = 0.90.

Fig
Fig. S3.Dependence on the colloidal packing fraction η of the effective surface charge density σ * of a colloid in an fcc crystal as a fraction of the surface charge density σ ∞ = σ C = 0.01 C/m 2 of a free isolated colloid.An increase of up to ∼ 100% in σ * is already visible when the cell radius b in the direction of the nearest neighbour just extends one or two Debye lengths λ D beyond the colloid of radius a, i.e. when the nearest neighbour distance b − a ∈ [0, 2λ D ].

Fig. S5 .
Fig. S5.Schematic depiction of in situ nanoparticle assembly in the microchannel.
we show measured hysteresis loops for 3 different frequencies for channels of lengths 50 µm, 100 µm and 150 µm.Additional measurements were conducted for intermediate frequencies in order to find upper and lower bound for the frequency which exhibits the most open hysteresis loop.All loops are shown in Fig. S13, Fig. S14 and Fig. S15 for channel lengths of 50 µm, 100 µm and 150 µm respectively.

0
Fig.S6.Six measurements of the normalized channel conductance for all 2 4 = 16 different voltage pulse trains, obtained from three devices.Each column of two figures corresponds to a device.The average for each different pulse train is used to create Fig.4(a).Standard deviations for the conductance measurement after the fourth (last) pulse were also calculated with these six measurements and shown in TableI.

Fig. S7 .
Fig. S7.A 50-cycle repeat of the 0101 voltage pulse train.With (left) the current measurement during all cycles, (top right) all cycles overlaid with the light grey spread the raw overlaid data, the dark grey line the average of the measurements, and the variability characterised by the standard deviation in light blue.(bottom right) The resulting normalized conductances determined by the five read-pulses with their respective standard deviations.

Fig. S9 .
Fig. S9.Measurements of the normalized channel conductance for all 2 4 = 16 different voltage pulse trains, but here the channel base is grounded and the voltage is applied at the tip.

Fig. S10 .
Fig.S10.Voltage pulse train measurements with the ground at the tip and the applied voltage at the channel base.A "0" corresponds to a pulse of 2 V, while a "1" corresponds to a pulse of -5 V. Pulse duration and interval are 0.75 s, with read pulses of 1 V and 50 ms duration.Voltage pulse trains corresponding to the bit-strings 1111, 0101 and 0011 were repeated 5 times, where we show the 5 individual measurements (light grey), the average of the measured current (black) and the calculated normalized conductances in the bar plots, averaged over the 5 measurements.The error bars depict the measured standard deviations.

Fig. S11 .
Fig.S11.Schematic depiction of how the read-pulses are used to calculate the channel conductance.The current measurements during the read pulse are averaged.The difference with the average current just before the read pulse then yields the channel conductance after dividing by the applied voltage.Before each voltage pulse train, a read pulse is applied to obtain the base conductance g 0 , which is used to normalise the measurements.

Fig. S12 .
Fig. S12.Simple single digit numbers used for classification in the main text as shown in Fig. 3(b) and Fig. 3(c).

Fig. S13 .
Fig. S13.Current-voltage hysteresis loops for the 50 µm channel as a result of a sinusoidal voltage over the channel of amplitude 5 V for the various frequencies shown.Measurements were gathered during five periods, depicted as the light grey graphs, where the average of the measurements is shown as a black graph.Enclosed areas were calculated with f max = 0.9 Hz (highlighted yellow) exhibiting the most open hysteresis loop.

Fig. S14 .
Fig. S14.Current-voltage hysteresis loops for the 100 µm channel as a result of a sinusoidal voltage over the channel of amplitude 5 V for the various frequencies shown.Measurements were gathered during five periods, depicted as the light grey graphs, where the average of the measurements is shown as a black graph.Enclosed areas were calculated with f max = 0.215 Hz (highlighted yellow) exhibiting the most open hysteresis loop.

Fig. S15 .
Fig. S15.Current-voltage hysteresis loops for the 150 µm channel as a result of a sinusoidal voltage over the channel of amplitude 5 V for the various frequencies shown.Measurements were gathered during five periods, depicted as the light grey graphs, where the average of the measurements is shown as a black graph.Enclosed areas were calculated with f max = 0.1 Hz (highlighted yellow) exhibiting the most open hysteresis loop.In the first run of measurements, 0.05 Hz yielded a hysteresis loop with an enclosed area very close to 0.1 Hz.Both frequencies have been tested various more times, where the 0.1 Hz loops were consistently more open, one of these experiments was used for Fig. 1(c).The other 0.05 Hz measurements are shown here.
, we demonstrate the stability and reproducibil-Memory timescale τ for all channels determined by the frequency f max for which the enclosed area inside the loop is maximal via the relation 2π f max τ = 1 [51].The measured frequencies around each f max yield upper and lower bounds, shown here as error bars.

TABLE I .
. The standard deviations of the normalized channel conductance g/g 0 for each different bit-string, determined through the six measurements per bit-string shown in Fig.S6.