Physical echo state network based on the nonlinearity and dynamic response of ambipolar heterostructure transistors

Domain: physics
Published: 2026-02-28
Editorial outcome: major_revision

Abstract
In the field of neuromorphic computing, time-series prediction poses a significant challenge to recurrent neural network architectures, often requiring task-specific customization that limit the development of general-purpose computing platforms. In this work, we implement a physical echo-state network (ESN) using ambipolar organic–inorganic heterostructure transistors to form its reservoir layer. Leveraging the ambipolar nature of the transistor, its variable-resistance region enables sparse matrix operations, while the saturation region provides tanh-like nonlinearity, making it well-suited for implementing both synaptic weighting and neuronal activation in an ESN. Additionally, its dynamic response naturally introduces temporal attributes. Thus, it can serve as a neuromorphic computing model for time-series tasks. Without the involvement of dynamic mechanisms, it is capable of performing image recognition, time-series prediction, and multimodal recognition tasks. When dynamic mechanisms are incorporated, the model achieves an accuracy of 96.98% on the MNIST handwritten digit dataset and 86.67% on the Fashion-MNIST dataset. This work offers a neuromorphic computing architecture, providing insights for tasks such as nonlinear mapping and time-series prediction.

Introduction
With the rapid development of artificial intelligence, computer chips are increasingly presenting challenges to humanity1–5. With the advancement of manufacturing processes, the size of transistors is approaching physical limits, meaning it is difficult to increase the density of transistors per unit area6. This has led to rising chip manufacturing costs, significant heat generation, and increased energy consumption7,8. In the von Neumann architecture, the separation of memory and processor causes memory bottleneck issues. To solve these problems, neuromorphic architectures have been proposed9–11. These novel computing architectures are inspired by the human brain and train neural network in a manner that mimics brain function12–15. Recently, researchers have utilized physical devices to simulate biological synapses and neurons, connecting these devices to construct neuromorphic neural networks. These networks emulate the brain’s working mechanisms, significantly enhancing efficiency and reducing energy consumption.

In a traditional sense, physical devices have been used as synapses to implement physical weights and as neurons to perform data mapping functions16–22. Neural networks composed of these devices have significantly advanced the field of image recognition. However, there is relatively limited research on using recurrent neural networks to address time series problems. This is partly because the training process is more complex, requiring a greater range of conductive states from physical devices. To overcome this challenge, researchers have begun exploring physical reservoir computing23–28. This approach involves constructing a reservoir with memory to replace the memory units in recurrent neural networks. Previous studies have utilized the nonlinear conductive states of physical devices to capture information from time series29. While this method is effective, it also presents limitations, such as difficulties in addressing binarization and multi-bit issues. An echo state network (ESN) is a specialized type of recurrent neural network that includes a reservoir layer30–32. The reservoir layer is constructed with a sparse matrix with fixed weights and neurons, which do not require training. The data is processed by the reservoir, resulting in reservoir states with high-dimensional features. The neural network excels at time series prediction tasks and can also classify image data that has been transformed into a time series. Consequently, reservoir layers constructed with physical devices hold significant potential.

In this work, we design ambipolar transistors based on poly(3-hexylthiophene) (P3HT)–zinc oxide (ZnO) heterostructures for implementing a physical echo-state network (ESN). Serving simultaneously as synaptic and neuronal elements within the ESN reservoir, these transistors generate high-dimensional reservoir states from electrical inputs, which are subsequently classified by a simple readout layer for recognition or prediction tasks. Across diverse machine learning benchmarks, the reservoir states produced by the ambipolar ESN consistently surpass those of unipolar counterparts, highlighting the significance of ambipolar conduction in reservoir construction. We evaluate the physical ESN in both static and dynamic operational modes. In static mode, the system achieves 95.27% accuracy on the MNIST handwritten digits dataset and 86.10% on Fashion-MNIST. For time-series prediction, it delivers normalized root-mean-square errors (NRMSE) of 0.15%, 0.22%, and 0.19% along the x-, y-, and z-axes of the Lorenz attractor, respectively. In typhoon trajectory forecasting, the ESN exhibits optimal performance with a 500×500 reservoir, and attains 89.7% accuracy in multimodal classification tasks. In dynamic mode, data are preprocessed through the intrinsic temporal response of the device, embedding temporal correlations before ESN training. This approach enhances system performance, yielding 96.98% accuracy on MNIST digits and 86.67% on Fashion-MNIST. Moreover, in Lorenz attractor prediction, the dynamic mode achieves NRMSEs of 0.165, 0.28, and 0.206% across the three axes—performance comparable to that attained with a forgetting factor in conventional ESNs. Together, these results demonstrate that ambipolar transistors constitute a promising hardware platform for high-performance neuromorphic computing.

Results
The neural network comprises three distinct layers: the input layer, the reservoir layer and the readout layer (Fig. 1a). The input layer is used for inputting data and transforming it to other dimensions. It allows the data to be fed into the reservoir layer. The reservoir is a layer constructed with a sparse matrix and neurons. It has the function of extracting features to map data to high dimensions and can also record historical information. It is responsible for outputting the reservoir states of the ESN. The readout layer is the layer that interprets the reservoir states. The reservoir states can be quickly solved through linear regression33–35. At the same time, other machine models can also be utilized for processing, such as Support Vector Machines (SVM) or deep neural networks. The cross-sectional morphology of the transistors with various thicknesses of ZnO are shown (Fig. 1b). The thickness of P3HT deposited on ZnO of varying thicknesses is almost consistent, close to 70 nm. The thickness of ZnO has a minor effect on the deposition of P3HT, and the semiconductor performance of P3HT is almost uniform. The high-resolution image shows that the thicknesses of ZnO are 17.8, 12.2, and 8.0 nm, respectively, corresponding to the n-type, ambipolar and p-type transistors. The surface morphology of the ZnO and P3HT layers of the transistors are captured (Fig. 1c). The AFM images reveal that both the ZnO and P3HT layers are relatively smooth. The root mean square (RMS) surface roughness of the ZnO layers for the n-type transistor, ambipolar transistor and p-type transistor are 935, 718 and 704 pm, respectively. Correspondingly, the surface roughness of the P3HT layers are 4.67, 3.36 and 4.12 nm, respectively.

The output performance of the transistors is illustrated (Fig. 1d). The p-type silicon wafer acts as the gate, and the 100 nm SiO2 layer serves as the gate dielectric. The ZnO layer is deposited on top, serving as the n-type channel material for the transistor. Subsequently, the P3HT layer is deposited on the ZnO layer, acting as the p-type channel material. An Au electrode is then deposited on the P3HT layer. The ambipolar transistor can transport electrons or holes depending on the applied gate voltage. For the n-type transistor, when the gate voltage (Vg) is less than 0, electrons in the ZnO layer are captured by O2− ions, and holes in the P3HT layer are released by the conjugated chain36–38. Given that the ZnO layer is 17.8 nm thick, the electron concentration is higher than the hole concentration, resulting in the transistor being in an off state. When Vg is greater than 0, electron pairs in the ZnO layer are released by O2− ions, and holes in the P3HT layer are captured by the conjugated chain. As the drain-source voltage (Vds) increases, the drain-source current (Ids) rises linearly until it reaches the pinch-off point, after which Ids tends to saturate, indicating the transistor is in an on state. Ids reaches approximately 0.93 µA at Vg = 30 V and Vds = 30 V. For the ambipolar transistor, when Vg is less than 0, electrons in the ZnO layer are captured by O2− ions, and holes in the P3HT layer are released by the conjugated chain. As Vds increases, Ids rises linearly until the pinch-off point, after which it saturates, indicating an on state. Ids is approximately −0.33 µA when Vg = −30 V and Vds = −30 V. When Vg is greater than 0, electrons in the ZnO layer are released by O2− ions, and holes in the P3HT layer are captured by the conjugated chain. Ids again rises linearly with Vds until saturation, indicating an on state. Ids reaches about 0.33 µA at Vg = 30 V and Vds = 30 V. The balanced performance of the ZnO and P3HT layers allow the transistor to exhibit ambipolar properties. For the p-type transistor, when Vg is less than 0, electrons in the ZnO layer are captured by O2− ions, and holes in the P3HT layer are released by the conjugated chain. As Vds increases, Ids rises linearly until the pinch-off point, after which it saturates, indicating an on state. Ids is approximately −0.34 µA when Vg = −30 V and Vds = −30 V. When Vg is greater than 0, electron pairs in the ZnO layer are released by O2− ions, and holes in the P3HT layer are captured by the conjugated chain. Due to the ZnO layer being 8.0 nm thick, its electron concentration is lower than the hole concentration in the P3HT layer, resulting in the transistor being in an off state. Due to the presence of traps within P3HT and ZnO, its transfer characteristic exhibits hysteresis (Supplementary Fig. S1).

In fabricating the array device, we first employed a 100-nm-thick Al2O3 layer as the gate dielectric. The array structure (Fig. 1e) and the corresponding hardware system (Fig. 1f) are shown, along with the memory characteristics of the array (Supplementary Figs. S2–S6). To further optimize device performance, we reduced the dielectric thickness to 30 nm and 10 nm (Supplementary Fig. S7). Thinning the dielectric layer effectively lowers the operating voltage while preserving ambipolar behavior. Transistors fabricated with the 100 nm Al2O3 dielectric also exhibit well-defined output characteristics consistent with their device type (Supplementary Fig. S8). To advance integration density, we further scaled the channel length down to 2 µm (Supplementary Fig. S9). With the dielectric thickness reduced to 10 nm, the device achieves an ultralow switching energy of 202 fJ (Supplementary Fig. S10 and Supplementary Table S1), providing a strong foundation for future low‑power electronics. Moreover, the device shows clear dynamic response (Supplementary Figs. S11, S12), supporting its potential in neuromorphic computing applications. Notably, the device’s transfer characteristic resembles the tanh function (Fig. 1g), enabling simultaneous data mapping and signal bounding. Specifically, in the linear region, it can perform sparse matrix multiplication, whereas in saturation, it offers nonlinearity and amplitude constraint. Its dynamic properties further introduce temporal correlation processing. Collectively, these features highlight the advantages of using ambipolar transistors to implement physical echo‑state networks for neuromorphic computing.

To demonstrate the capabilities of the physical ESN, we first conduct a letter recognition task involving six letters (O, Z, W, T, X, U). Each letter is repeated 2000 times and assigned a label from 0 to 5. The dataset is shuffled and split into 10,000 training and 2000 test samples. Letter images are flattened into time-series signals (Supplementary Fig. S13) and fed into an input layer with randomized signed weights, converting grayscale values into bipolar features. These features are then mapped to voltage inputs for a reservoir layer implemented as a 16×16 transistor matrix (Fig. 2a, b). Sparse connectivity is emulated by activating only 10% of the synaptic connections. Two key characteristics distinguish transistor-based synaptic multiplication. First, ambipolar transistors exhibit two linear regions under positive and negative gate biases, enabling a wider and more accurate multiplicative range compared to unipolar (n- or p-type) transistors, which operate in only one linear region. Second, ambipolar devices achieve balanced channel conduction through complementary ZnO and P3HT layers, mitigating the inherent conductivity imbalance of unipolar transistors. In a comparative matrix multiplication experiment (Supplementary Fig. S14), synaptic arrays encoded matrix elements as Vds (first matrix) and Vg (second matrix). Ambipolar synaptic arrays produced results closest to digital computation, outperforming unipolar counterparts (Supplementary Fig. S15).

When transistors act as neurons, their conductance characteristics and saturation behavior critically influence the reservoir. Ambipolar transistors show symmetric current saturation under both gate polarities beyond the pinch-off voltage, providing an inherent nonlinear activation function analogous to tanh. Under ±30 V time-series inputs (Supplementary Fig. S16), n-type neurons output only positive currents, p-type only negative currents, while ambipolar neurons generate balanced ambipolar responses. Consequently, reservoir states differ markedly: n-type ESNs produce outputs from 0 to 0.75 µA, p-type from 0 to –0.1 µA, and ambipolar ESNs from –0.1 to 0.1 µA (Fig. 2c). Despite these differences, all ESN variants achieve 100% accuracy in the letter recognition task (Fig. 2d and Supplementary Fig. S17). The complete system operation is further illustrated (Supplementary Fig. S18).

The MNIST datasets are used for classification39–41. These images are all 28×28 size and are fed into the physical ESN in the same way as above. 60,000 images are used as the training set, while the remaining 10,000 images serve as the test set. Thank to the inherently excellent nonlinearity of ambipolar transistors, even when performing direct image straightening, the data can be mapped to a high-dimensional space. This allows the readout layer to easily interpret the reservoir states, enabling a single-step solution with high accuracy. The ESN achieves an accuracy of 95.27% on the grayscale MNIST handwritten digit dataset and 86.10% on the MNIST fashion dataset (Supplementary Fig. S19–S23). This is not the ultimate limit of the architecture. To achieve higher accuracy, scaling up the sparse matrix size is required. Given its already high accuracy, further improvements would necessitate additional devices, introducing inherent limitations. Another method involves prepossessing the datasets with the dynamic response of the device and then using the prepossessing datasets to complete the recognition task. Since the dynamic response can provide features for the datasets, effectively adding dynamic nonlinearity, it can also yield favorable results. We further binarize the MNIST dataset and use 4-bit encoding to map the images, followed by processing with a physical ESN for recognition. Different pulse intervals in the 4-bit encoding are employed to impart more features (Supplementary Figs. S24 and S25). Information encoding is first established through a four-bit scheme (Fig. 3a). Different types of nonlinearity are achieved by adjusting the interval pulses. For example, when the interval is increased, the responses produced by the “0,0,0,0” and “1,0,1,0” encodings are different (Fig. 3b). This gives us further opportunities for exploration. When images are processed using encodings with different pulse intervals, different image mappings are obtained (Fig. 3c). When these mappings are processed by the physical ESN, they produce different reservoir space states (Fig. 3d). Both the digital dataset and the fashion dataset exhibit discrete reservoir states (Supplementary Fig. S26). However, their time attributes are incorporated due to preprocessing with different pulse intervals. This results in the potential for creating better performance opportunities. When the mappings converted from six pulse intervals are read out by the physical ESN, the accuracy reaches 96.98% (Fig. 3e). Similarly, on the fashion dataset, after dynamic processing, the accuracy also increases relatively to 86.67% (Supplementary Fig. S27), higher than 86.10%. Therefore, introducing nonlinearity and dynamic mechanisms through physical methods in MNIST dataset processing may serve as an effective approach to enhance performance while reducing computational costs.

For neuromorphic computing systems, the capability to forecast temporal and spatial sequences is essential for ESNs42–44. To evaluate this, we apply our ESN to predict the Lorenz attractor, using data at time step t as input and at t+2 as target. This task represents a second key challenge in time‑series prediction, an area where hardware‑based architectures are actively studied. We explore two training strategies. The first adopts the forgetting‑factor method commonly used in ESNs (Supplementary Fig. S28), where the factor α substitutes for explicit dynamic mechanisms. By tuning α, the physical ESN achieves normalized root‑mean‑square errors (NRMSE) of 0.165, 0.28, and 0.206% along the x‑, y‑, and z‑axes, respectively. The forgetting factor ensures continuity in reservoir states across time (Supplementary Figs. S29–S31). The second approach directly embeds temporal attributes by preprocessing the data through the dynamic response of the physical device before training (Fig. 4a). Here, training proceeds with a forgetting factor of 1, equivalent to a static mapping setup. The device’s dynamic response performs a nonlinear transformation on the Lorenz attractor, and the resulting time‑varying Ids serves as input to the ESN (Fig. 4b–d). By varying the sampling intervals during preprocessing, different dynamic features can be extracted for prediction. Though the reservoir states under this method appear discrete, they retain encoded temporal information (Supplementary Fig. S32). This preprocessing strategy yields NRMSE values of 0.165%, 0.28%, and 0.206% on the x‑, y‑, and z‑axes, respectively (Fig. 4e)—performance that closely matches, and in some dimensions slightly exceeds, that of the forgetting‑factor method. Consequently, the system achieves low prediction error and generates smooth Lorenz attractor trajectories (Fig. 4f).

The system is further applied to predicting typhoon-center imagery and trajectories. Using grayscale image sequences and longitude–latitude data from Typhoon “Kanu” (2023), the ambipolar ESN achieved the lowest image-prediction NRMSE of 10.88% and maintained balanced dynamic activity between positive and negative states even when the reservoir was scaled to 2000×2000 (Supplementary Figs. S33–S35). For trajectory prediction, the ambipolar ESN with a 500×500 reservoir best fitted the test data, notably attaining a stable latitude‑prediction error in the range of 8.34–8.46% (Supplementary Fig. S36). These results collectively demonstrate that the ambipolar ESN excels in multi‑dimensional spatiotemporal forecasting, owing to its rich state diversity and the synergistic interplay between positive and negative activations, which together enhance dynamic‑system modeling.

To further validate the multimodal classification capability of our system, we employ electrocardiogram (ECG) signals. The signals are converted into time–frequency images via short-time Fourier transform (STFT)45,46 to form the dataset (Fig. 5a). This process yields two complementary data modalities: the original 1D waveforms and the derived 2D spectrograms. Both modalities are processed independently by the reservoir, and the resulting states are classified through a shared readout layer. Recognition can be performed on either modality alone (signal mode or image mode) or on their combined representations (dual mode) (Fig. 5b). When processed by the physical ESN, the signal and image modalities produce distinct reservoir state patterns (Fig. 5c). Combining these states enhances overall discriminability. We systematically evaluate ESN performance across reservoir scales. The ambipolar ESN consistently achieves the highest accuracy in all modes: signal mode (82.80–88.38%), image mode (84.06–87.46%), and dual mode (87.62–89.77%), with accuracy generally improving as the reservoir size increases until saturation beyond 500×500 (Fig. 5d). The p‑type ESN delivers stable, scale‑insensitive performance, particularly in dual mode where accuracy remains nearly constant (87.39–88.35%). In contrast, the n‑type ESN peaks in dual mode at 87.81% with a 200×200 reservoir, beyond which performance degrades, revealing limited scalability. These results underscore the ambipolar ESN as the most robust and scalable architecture for multimodal processing, achieving superior accuracy while maintaining consistent performance scaling. This advantage stems from its more balanced reservoir state distribution across positive and negative activations (Supplementary Figs. S37–39).

Finally, a task utilizing the Mosei dataset for multimodal recognition was conducted. The audio, text, and visual features from Mosei were processed using a pooling mechanism and then classified with SVM. For tri-modal classification (Supplementary Fig. S40), the accuracy without the ambipolar ESN is 71.33%, whereas with the ambipolar ESN, it improves to 76.07%. This demonstrates that the ambipolar ESN significantly enhances the accuracy of multimodal recognition. For dual-modal classification, the accuracies for visual/audio, visual/text, and text/audio combinations are 70.86, 76.46, and 76.01%, respectively. For monomodal classification, the accuracies for text, audio, and visual modalities are 76.48, 71.01, and 71.03%, respectively. Using a 10 nm alumina dielectric layer allows the transistor array to operate at lower voltages, yet introduces process-related challenges in uniformity and endurance (Supplementary Fig. S41). When the pulse interval is increased sufficiently, the dynamic response of the array gradually transitions toward static behavior (Supplementary Fig. S42). In image recognition tasks, relying solely on the final current value after dynamic mapping proves inadequate. For instance, in 5‑bit English letter recognition using only the steady‑state response as the training feature, we observed misclassification patterns (Supplementary Figs. S43, S44), attributable to the loss of transient information. Compared to binarized inputs, conventional 4‑bit digital encoding yields lower accuracy due to feature compression (Supplementary Fig. S45). However, by fully utilizing the temporal 4‑bit response data, recognition accuracy improves from 84.92% to 87.94% (Supplementary Fig. S46), confirming that dynamic nonlinearity enhances binarized image classification. To further boost accuracy, we introduced sparse projection and activation‑based feature extraction, combining these derived features with raw data for training. This strategy led to additional gains in recognition performance (Supplementary Figs. S47–S55). The array architecture can physically implement such operations through its tunable linear regions and inherent nonlinearity, supporting efficient and accurate image recognition.

Nevertheless, while the system exhibits strong recognition capability, further accuracy improvements demand increased computational resources. Moreover, time‑series prediction currently relies on a forgetting factor, implying that the computational model must support online iteration. These limitations can be mitigated by leveraging the device’s intrinsic dynamic response. The additional features generated through dynamics enhance physical ESN performance without digital overhead. In time‑series prediction, the device’s dynamic equation mapping enables single‑pass inference, eliminating the need for online iterative updates.

In summary, ambipolar transistors, with their balanced linear characteristics and bounded output response, are particularly suitable for emulating both the linear weighting of synapses and the nonlinear activation of neurons in a reservoir layer. Since the weights in the reservoir layer are fixed and require no training, such devices represent a promising candidate for constructing recurrent neural networks with inherent dynamic constraints, highlighting considerable research potential. In contrast, p-type and n-type transistors exhibit limitations in building high-performance reservoirs due to their weaker capability to bound signal amplitudes. Under static operation, the system has demonstrated strong potential in tasks such as image recognition, time-series prediction, and multimodal classification. In dynamic mode, by leveraging the intrinsic transient response of the device, temporal correlations can be embedded into the data, leading to further improvements in image recognition performance. Particularly for time-series forecasting, the dynamic mechanism can replace the forgetting factor used in conventional ESNs, transforming iterative online computation into single-step forward processing, thereby enhancing overall computational efficiency.

Methods
Fabrication of the device and array
The bottom gate electrodes were fabricated via physical vapor deposition. Aluminum oxide was subsequently deposited as the gate dielectric using atomic layer deposition. For the semiconductor layers, ZnO was prepared by dissolving the precursor in ammonia to form a uniform solution, while P3HT was dissolved in toluene and thoroughly mixed by stirring and heating to ensure homogeneity. The substrate was placed on a spin coater, and the ZnO precursor solution was dispensed onto it. Spin-coating was performed sequentially at 1500 rpm, 2000 rpm, and 2500 rpm, each for 30 s. After spin-coating, the substrate was annealed at 150 °C for 40 min to remove the solvent and crystallize the ZnO thin film. The P3HT solution was then applied onto the ZnO layer, with the volume carefully controlled to ensure film uniformity, followed by spin-coating at 1500 rpm for 30 s. The resulting film was annealed at 100 °C for 20 min. Device and array structures were defined by covering the thin film with a mask and depositing electrodes via physical vapor deposition. All fabrication steps were conducted in a nitrogen atmosphere. The completed array was encapsulated by covering it with a quartz glass sheet and sealing the edges with UV-curable adhesive.

Control of the array
Specific rows are selected by applying a feature signal that connects the input voltage to the drain terminals of all transistors in the chosen row via a digital-to-analog converter (DAC). Simultaneously, specific columns are addressed by applying a column selection signal to the gate terminals of all transistors in the selected column. This addressing scheme enables individual control of each unit within the array, facilitating the implementation of matrix multiplication and activation functions.

Measurement of the array
A custom program was used to sequentially access and test each node in the array, thereby evaluating the operational performance and stability of the array. To measure the output, a resistor with a known value (1 MΩ) was connected in series at the output terminal of the transistor array, allowing the output current to be converted into a voltage signal for further analysis.

Device characterization
The electrical properties of the device were tested using a semiconductor analyzer (Keysight B2902B). The device was first cutting using Focused Ion Beam Scanning Electron Microscopy (Helios 5UC 03040707), and then the cross-sectional morphology was taken by Field Emission Transmission Electron Microscopy (JEOL-F200). The surface morphology of the device was characterized using an Atomic Force Microscope (Cypher S/Oxford Instruments Asylum Research).

Equation
The formula for the status of the reservoir is given as follows:
U(t) = U(t−1)(1−α) + αtanh(Xt Win + U(t−1) Wres) (1)

Where U(t), U(t−1), α (alpha), Xt, Win, Wres are represent for reservoir state at time step t, reservoir state at the previous time step at t−1, forgetting rate, input at time step, input weight matrix and reservoir weight matrix.

The solution method for the reading layer uses the linear regression formula:
Wout = (UTU)−1(UTylabel) (2)

Wout is the readout layer matrix, U is the reservoir state, and ylabel is the label value.

The formula for the Lorentz attractor is given as follows:
dx/dt = σ(y−x)
dy/dt = x(ρ−z)−y (3)
dz/dt = xy−βz

σ, ρ and β are respectively the Prandtl number, the Rayleigh number and the geometric factor.

In addition, the gradient descent algorithm is used in the prediction of the typhoon center, with Adam as the optimizer. The loss function used is the Root Mean Square Error (RMSE), and its formula is as follows:
RMSE = 1/n Σ (ylabel − ypredict)2 (4)

ypredict is the predicted result.

The formula for NRMSE is given as follows:
NRMSE = 1/(ymax − ymin) * sqrt(1/n Σ (ylabel − ypredict)2) (5)

ymax, ymin are max value of ylabel, min value of ylabel, respectively.

Dataset
The MNIST handwritten digits and fashion datasets are provided by the National Institute of Standards and Technology. Each dataset consists of 60,000 images for training and 10,000 images for testing. The Lorenz attractor curve comprises 10,000 points, with the first 5000 points used as the training set and the remaining 5000 points as the test set. In this task, the input consists of points 1 to 4998, while the target includes points 3 to 5000. The images of the typhoon center and its movement trajectory are derived from Typhoon Khanun in 2023. The ECG dataset is sourced from the Massachusetts Institute of Technology-Beth Israel Hospital, and includes 17 categories based on their annotations. To ensured dataset validity, ECG signals with more than 100 samples are selected, resulting in a total of 11,811 samples. The first 10,000 samples are used for training, and the remainder are used for testing. The CMU-MOSEI dataset, developed by the Multimodal Communication and Language Analysis Lab at Carnegie Mellon University (CMU), uses binary classification labels: 0 for values less than 0, and 1 for values greater than or equal to 0. In addition, the support vector machine employed in the final section of the study uses a linear kernel function, with the C parameter set to 0.0001.

Networks
The English alphabet recognition task and Lorenz attractor prediction are achieved through a 16×16 ambipolar transistor array, with read/write operations on the array controlled by DACs and an ADC via STM32. MNIST handwritten digits, the Fashion dataset, Lorenz attractor prediction, ECG data classification, and multimodal recognition with the MOSEI dataset are implemented using Python+PyTorch simulations. The simulation is performed by an operational model of the device. The specific physical ESN network structure is provided in the Supplementary Table S2.

Data availability
The physical data generated in this study have been deposited in the Zenodo database [https://doi.org/10.5281/zenodo.18475506].

Code availability
The code is available at Zenodo [https://doi.org/10.5281/zenodo.18459291].

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "This paper describes the implementation of physical reservoir computing by utilizing ZnO/P3HT ambipolar transistor array as a physical reservoir. The proposed ZnO/P3HT device can be configured in either n-type, ambipolar, or p-type. The authors found that for the purpose of vector-matrix multiplication and accurate physical reservoir computing, the ambipolar characteristic is highly desirable. The single-polar n-type or p-type devices resulted in single-signed outputs, leading to limited reservoir states. As demonstrations of the approach, tasks such as letter recognition, MNIST dataset classification were done as an example of processing static data. Another set of tasks such as prediction of typhoon route or classifying ECG signals were presented as demonstrations of learning dynamic systems.\n\nThis paper reveals the importance of ambipolar properties in building neural networks, and presents a novel method of utilizing transistors for physical reservoir computing. However, most importantly, the descriptions of how the system works are unclear. I believe that they should be addressed before additional review and publication.\n- It seems that a single device is used repeatedly to implement the physical array operation. Given the device structure, construction of the array is not straightforward. This should be clearly stated. I wonder if the computation is done by simulation (based on the model of the device operation), or by measuring the response from the single device repeatedly to construct the response from the array.\n- What is the exact structure of the network including the reservoir? It seems that the input data scale is first reduced to match the size of the reservoir input using a single-layer neural network. Then, what does the 'step' mean in figures for reservoir state output? (Fig. 2c, 3b, 3c) Probably, the step means the sample number.\n- Once again, detailed network structure for temporal and spatial sequence forecasting is not clear from the text. How is temporal data being processed?", "reviewer_2": "The manuscript presents a physical Echo State Network utilizing ambipolar organic-inorganic heterostructure transistors. Key achievements include high accuracy on multiple tasks, such as the typhoon trajectory prediction. The ambipolar transistors demonstrate superior performance compared to n-type and p-type variants, attributed to their balanced conductivity and dual linear regions, enabling effective feature mapping and reservoir state diversity. This work bridges neuromorphic hardware design with machine learning applications, offering a pathway toward energy-efficient computing platforms.\n\nThe manuscript can be improved by addressing the following concerns:\n1. As is introduced in the manuscript, the PRC involves constructing a reservoir with memory to replace the memory units in recurrent neural networks. It would be reasonable to discuss the memory characteristic of the devices presented. My understanding is that different from the \"ideal\" transistors for CMOS ICs with no or little hysteresis, the transistors developed in this work should present hysteresis in the transfer curves due to ion immigration or trap states, which can be the source of coherent memory.\n2. The size of the physical reservoir, for example \"500*500\", is not clearly explained. Does the size of the physical reservoir correspond to the array dimensions of the devices? How are these device arrays fabricated, controlled, and measured? I would suggest that the authors provide a description of the operation of the system.\n3. The writing can be further improved by condensing the manuscript.\nOther suggestions for the authors:\n4. The color bar is missing in several figures, such as Figure 2.\n5. The text in Figure 2a is horizontal, which seems unnecessary.", "reviewer_3": "In this manuscript, the authors present an ambipolar organic-inorganic heterostructure transistor with potential applications in reservoir computing systems. While various applications have been demonstrated, the work does not introduce any fundamental innovations or address key challenges in the field. Additionally, the manuscript suffers from significant issues in grammar and logical flow, which further undermine its clarity and overall presentation. The current version does not meet the standards expected of journals in the Nature Communications. Below are my detailed concerns regarding the manuscript:\n1. The manuscript lacks comprehensive device testing. The authors need to provide data on endurance, retention, uniformity, switching speed, and environmental stability to evaluate whether the device is sufficiently robust and suitable for computing applications.\n2. The authors mention in the introduction that transistor sizes are approaching physical limits. What is the actual size of the proposed device? Does this new design contribute to reducing transistor dimensions or increasing the density of transistors per unit area?\n3. Reservoir computing typically requires specific temporal dynamics. Do the proposed devices exhibit any intrinsic temporal properties essential for reservoir computing?\n4. Can the operating voltage of this device be reduced in future implementations? The authors should provide an analysis of potential strategies for lowering power consumption.\n5. It is unclear what hardware was used for the later application demonstrations. If the applications were performed through simulations, the authors should explicitly state this to help readers better assess the significance of the work.\n6. How were the input pixels translated into time-series sequences in Figure 2? A more detailed explanation is necessary.\n7. The authors appear to read out the current at each time step, suggesting that all temporal information originates from the input signal. However, the transistor array itself does not seem to perform any nonlinear transformation in the temporal domain. Can the authors clarify this?\n8. If the prediction is accurate and the dataset is uniform, why do the confusion matrices exhibit different color distributions? The authors should explain these discrepancies.\n9. The MNIST and Fashion MNIST datasets are inherently static. How did the authors convert these static datasets into temporal voltage signals for processing within the reservoir computing framework?\n10. The authors need to demonstrate that the reservoir computing component is indeed contributing to the classification performance. Additionally, the classification network used should be clearly described. There is a concern that the high accuracy is primarily attributed to the readout network rather than the reservoir itself."}, "author_response": "We sincerely thank the reviewer for your positive evaluation of our work and for highlighting the novelty and significance of our approach. We also appreciate your constructive feedback regarding the clarity of the system’s operational description. In response to your concerns, we have substantially revised the relevant sections of the manuscript to provide a clearer and more detailed explanation of how the physical reservoir computing system operates. We hope that these revisions address your concerns and improve the overall clarity and readability of the manuscript. Thank you again for your valuable comments, which have helped us to strengthen our work.\n\nComment 1:\nIt seems that a single device is used repeatedly to implement the physical array operation. Given the device structure, construction of the array is not straightforward. This should be clearly stated. I wonder if the computation is done by simulation (based on the model of the device operation), or by measuring the response from the single device repeatedly to construct the response from the array.\n\nResponse to Comment 1:\n1.1 It seems that a single device is used repeatedly to implement the physical array operation. Given the device structure, construction of the array is not straightforward.\nThank you very much for your insightful comment regarding the implementation of the physical array operation. In our work, to ensure the practical applicability of our design, a physical ESN system operating on a 16×16 array was developed. To further investigate the potential of larger-scale arrays, we have also performed extensive simulations. These simulations allow us to model and analyze the behavior of the system in larger configurations, providing valuable insights into its scalability and performance beyond the current experimental limitations.\nWe have clarified these points in the revised manuscript to provide a more comprehensive understanding of both the experimental and simulation aspects of our work. Thank you again for your valuable feedback.\n\nFig. 1 Structure and characterization of the ambipolar transistor. a Schematic diagram of ESN. b TEM image of the cross-sectional morphologies under 100,000x magnification: n-type transistor, amipolar transistor and p-type transistor cross-section. High-resolution cross-sectional morphologies of ZnO layers in n-type transistor, ambipolar transistor and p-type transistor. c Surface morphologies observed under AFM: ZnO layers in n-type, ambipolar, and p-type transistors, as well as P3HT layers in n-type, ambipolar, and p-type transistors. d Output curves of n-type, ambipolar, and p-type transistors using 100 nm SiO as the dielectric layer. e Optical image of the ambipolar transistor array. f Optical image of the physical ESN system. g Output characteristic of ambipolar transistor in the array and tanh function.\n\n1.2 I wonder if the computation is done by simulation (based on the model of the device operation), or by measuring the response from the single device repeatedly to construct the response from the array.\nThank you for your thoughtful question regarding whether the computations were performed through simulation or by direct measurement from the physical device. For the English letter recognition and Lorenz attractor (16×16) tasks, we conducted rigorous validation using the physical 16×16 array. The device responses were directly measured from the physical system, allowing us to demonstrate the capability of our array to output reservoir states. Due to current technological limitations in fabricating large-scale arrays, it is not yet feasible to experimentally realize arrays of the size required for certain tasks, such as MNIST dataset recognition. For these larger-scale tasks, we employed predictive modeling and simulations based on experimentally validated device models. These simulations enabled us to explore the potential performance and scalability of our approach in more extensive configurations.\nWe have clarified these distinctions between experimental and simulation-based results in the revised manuscript. Thank you again for your valuable feedback.\n\nComment 2:\nWhat is the exact structure of the network including the reservoir? It seems that the input data scale is first reduced to match the size of the reservoir input using a single-layer neural network. Then, what does the 'step' mean in figures for reservoir state output? (Fig. 2c, 3b, 3c) Probably, the step means the sample number.\n\nResponse to Comment 2:\n2.1 What is the exact structure of the network including the reservoir?\nThank you for your question regarding the exact structure of the network. The network consists of an input layer, a hidden layer, and an output layer. The input layer contains input nodes and an input weight matrix, where the input nodes receive external input data, and the input weight matrix transmits the weighted input data to the hidden layer. The hidden layer is the core of the ESN, comprising a large number of hidden nodes and an internal weight matrix. The hidden nodes are interconnected to form a complex dynamic system capable of capturing long-term dependencies in time-series data, while the internal weight matrix is randomly initialized and remains unchanged during the training process. The output layer maps the states of the hidden layer to the final output results through an output weight matrix, which is the only parameter that needs to be learned during the ESN training process.R1,R2\nReference:\nR1. Lukoševičius, M. & Jaeger, H. Reservoir computing approaches to recurrent neural network training. Comput. Sci. Rev. 3, 127-149 (2009).\nR2. Boccato, L., Lopes, A., Attux, R. & Von Zuben, F. J. An extended echo state network using Volterra filtering and principal component analysis. Neural Netw. 32, 292-302 (2012).\n\n2.2 It seems that the input data scale is first reduced to match the size of the reservoir input using a single-layer neural network.\nThank you for your observation regarding the role of the input layer in our network. The input layer is an important component of neural networks, responsible for receiving and transmitting input data. From the perspectives of functionality and computational capability, the input layer cannot be regarded as an independent single-layer neural network. It serves more as a data interface, ensuring that input data can correctly enter the network for processing. Thus, the input layer itself does not possess the core computational characteristics of a neural network. This design is consistent with the ESN framework, as described in \"The 'echo state' approach to analysing and training recurrent neural networks-with an Erratum note\".R3\nReference:\nR3. Jaeger H. The “echo state” approach to analysing and training recurrent neural networks-with an erratum note. Bonn, Germany: German national research center for information technology gmd technical report, 13, 148(34) (2001).\n\n2.3 Then, what does the 'step' mean in figures for reservoir state output?\nThank you for your question regarding the meaning of 'step' in the figures for reservoir state output. For a time series, it corresponds to the time step n. Its input at each time step generates the corresponding reservoir space state, as well as the output. As you mentioned, it is indeed a sample at a certain moment.R3\n\nComment 3:\nOnce again, detailed network structure for temporal and spatial sequence forecasting is not clear from the text. How is temporal data being processed?\n\nResponse to Comment 3:\nThank you very much for your valuable comment. We appreciate your interest in the detailed network structure and the processing of temporal data. To clarify, when forecasting temporal and spatial sequences, such as in the Lorenz attractor prediction experiment, our system operates as follows:\nFor algorithmic logic, when dealing with time series problems, one data point is fed into the network at each time step. This data will be dimensionally augmented by a fixed random one-dimensional matrix.\nFor time series data, the following is provided:\nu = [a ⋯ a ]\n0 t\nwhere u represents the dataset, and a to a represent the elements of the entire dataset\n0 t\nat a certain moment.\nb\nin\nb\nn\nWhen t=0, the reservoir space state (Res) equals a vT.\nRes = tanh(a vT +Res W )\nt t t−1 res\nW is the fixed weight matrix.R3\nres\nAt the hardware level, input data undergoes linear transformation before being encoded into pulse sequences. The magnitude of the data is represented by the intensity of the pulses. This encoding method enables the array to process data in a sparse and efficient manner.\nBy mapping the elements of the input matrix to the transistors in the reservoir layer, matrix multiplication can be transformed into spike transmission and integration operations between nodes. For example, the multiplication of input reservoir matrix Res and reservoir layer matrix W can be achieved by mapping the elements of\nt−1 res\nRes to spike intensities, inputting them from the drain-source terminal array, and\nt−1\nadjusting the weights from the gate-source terminal. R4,R5\nReference:\nR4. Lawrence, S., Yandapalli, A. & Rao, S. Matrix multiplication by neuromorphic computing. Neurocomputing 431, 179-187 (2021).\nR5. Kim, S., Lee, Y., Kim, H.-D. & Choi, S.-J. Precision-extension technique for accurate vector–matrix multiplication with a CNT transistor crossbar array. Nanoscale 11, 21449-21457 (2019).\nTo facilitate a better understanding of these aspects, we have added detailed explanations and illustrative figures in the supplementary materials. Additionally, we have included experimental results of Lorenz attractor predictions implemented with our hardware system.\nThis operation employs a 16×16 physical array to accomplish the prediction of the Lorenz attractor. As a time-series signal, we utilize the data at time t, denoted as a.\nt\nFirst, it is expanded from a fixed 1×16-dimensional matrix into 16-dimensional data. Then, this 16-dimensional data is processed by our array operation. The array represents the matrix of W . The weights of W are fixed and exhibit sparsity. Since the array is ambipolar, it\nres res\ncan perform both positive and negative computational operations, with its weights derived from the gate voltage application. Leveraging the physical characteristics of the linear region, it accomplishes multiplication computations.\nThe input is applied as a voltage from the drain, facilitated by a multi-precision DAC. Simultaneously, another multi-precision DAC is used to apply voltage for weight configuration. Thus, by utilizing the ambipolar property of the array, the multiplication results for each column can be obtained by reading the I (V /R ) of each column\nds read ref\noutput via an ADC, thereby completing the computation. The matrix computation results are then converted back into voltage by a DAC. At this stage, the first row of the array is selected to perform the activation function operation. The I subsequently read by the\nds\nADC represent the reservoir spatial states. Furthermore, since the network requires memory for time-series prediction, dynamic response is employed to establish short-term memory and link historical correlations after converting I into voltage.\nds\nThese reservoir spatial states are processed via the Echo State Network (ESN) approach and then fed into the readout layer. Figure S18 illustrates the schematic diagram of the operational logic.\nFig. S18 operational logic schematic diagram.\nStep 1, Terminal dimension elevation operation.\nStep 2, Array multiplication operation.\nStep 3, Array activation function operation.\nStep 4, Array short-term memory operation.\nStep 5, Terminal readout operation.\nFig. 4 Temporal and spatial sequence forecasting. a The NRMSE prediction of the Lorenz attractor. b The reconstructed Lorenz curve. c NRMSE of the Lorenz attractor as predicted by n-type ESN, ambipolar ESN and p-type ESN with varying alpha values along the x-axis, y-axis, and z-axis. d Lorenz attractor curves as predicted by n-type ESN, ambipolar ESN and p-type ESN.\nWe predicted its behavior using the array, achieving NRMSE values of 0.24%, 0.49%, and 0.43% along the x, y, and z axes, respectively as shown in Fig. 4a and the reconstructed prediction curve is presented in Fig. 4b.\nIn addition, we provide all the network structures in this work as shown in Table S2.\nTable S2 Network structure of physical ESN for identification tasks.\nReadout layer\nInput layer Reservoir layer\n(train methods\nTask (Matrix (Matrix Operation\nand Matrix\ndimensions) dimensions)\ndimensions)\nLinear\nEnglish Letter (25,16) (16,16) regression Hardware\n(42, 6)\nLinear\nLorenz attractor (1,16) (16,16) regression Hardware\n(18,1)\nLinear\nLorenz attractor (1,100) (100,100) regression Simulation\n(102,1)\nAdam\nTyphoon image (625,200) (200,200) Simulation\n(826, 50265)\nLinear\nTyphoon track (10,20) (20,20) regression Simulation\n(31,1)\nLinear\nMNIST\n(784, 800) (800, 800) regression Simulation\ndigits/fashion\n(1785,10)\nLinear\nECG (100,150) (150,150) regression Simulation\n(251,1)\n\nResponse to Reviewer 2:\nWe sincerely thank you for your positive evaluation of our work and for highlighting its significance in bridging neuromorphic hardware design with machine learning applications. We greatly appreciate your constructive feedback and have carefully considered each of your suggestions. In response, we have thoroughly revised the manuscript and addressed your concerns point by point. We hope that these revisions have improved the clarity and quality of our work.\n\nComment 1:\nAs is introduced in the manuscript, the PRC involves constructing a reservoir with memory to replace the memory units in recurrent neural networks. It would be reasonable to discuss the memory characteristic of the devices presented. My understanding is that different from the \"ideal\" transistors for CMOS ICs with no or little hysteresis, the transistors developed in this work should present hysteresis in the transfer curves due to ion immigration or trap states, which can be the source of coherent memory.\n\nResponse to Comment 1:\nThank you very much for your insightful comment regarding the memory characteristics of our devices. In response to your suggestion, we have conducted additional hysteresis measurements on the ambipolar transistors. As anticipated, we observed clear hysteresis in the transfer curves, which we attribute to ion migration and trap states within the device. This hysteresis behavior indeed serves as a source of intrinsic memory, distinguishing our devices from conventional \"ideal\" CMOS transistors. We have included the results and a detailed discussion of these hysteresis measurements in the supplementary materials to further clarify the memory characteristics of our devices.\nThank you again for your valuable feedback, which has helped us to strengthen our manuscript.\n\nFig. S1 Transfer characteristics of the ambipolar transistor. (a) Transfer curve measured at negative drain-source voltage. (b) Transfer curve measured at positive drain-source voltage.\nThere may be a large number of trap states at the interface between P3HT and ZnO. These trap states can capture and release charge carriers, leading to hysteresis in the transfer curves. Specifically, when the direction of voltage scanning changes, the previously captured carriers will gradually be released, altering the current response of the device. P3HT is a p-type semiconductor, while ZnO is an n-type semiconductor, and the heterojunction formed between them affects charge transport. In the P3HT/ZnO system, the separation and transport efficiency of electrons and holes have a significant impact on the hysteresis of the transfer curves. If charge transport at the interface is incomplete or if there is charge accumulation, it can lead to hysteresis. The transfer characteristic curves of the device have been tested. The voltage scanning procedure starts from +30 V to -30 V, and then scans back from -30 V to +30 V. The device is tested under different V conditions. It is observed that hysteresis curves appear across different V\nds ds\nvalues.1-3\nReference:\n1. Ba, M. et al. Controlling of Hysteresis by Varying ZnO-Nanoparticles Amount in P3HT:ZnO Hybrid Thin-Film Transistor: Modeling. J. Electron. Mater. 52, 1203-1215 (2023).\n2. Zhang, L., Lin, H., Wu, Y. & Zhuo, S. Interfacial electron transfer of P3HT/PDI/ZnO nanocomposite and its application in visible-light detection. Chem. Phys. Lett. 661, 224-227 (2016).\n3. Liu, H., Li, M., Yang J. et al. In situ construction of conjugated polymer P3HT coupled hierarchical ZnO composite with Z-scheme enhanced visible-light photocatalytic activity. Mater. Res. Bull. 106, 19-27 (2018).\n\nComment 2:\nThe size of the physical reservoir, for example \"500×500\", is not clearly explained. Does the size of the physical reservoir correspond to the array dimensions of the devices? How are these device arrays fabricated, controlled, and measured? I would suggest that the authors provide a description of the operation of the system.\n\nResponse to Comment 2:\n2.1 The size of the physical reservoir, for example \"500×500\", is not clearly explained. Does the size of the physical reservoir correspond to the array dimensions of the devices?\nThank you for your question regarding the definition of the \"size of the physical reservoir\" in our work. In physical echo state networks (ESNs), the \"size of the physical reservoir\" (e.g., \"500×500\") refers to the number of neurons in the physical reservoir, which is directly related to the dimensions of the sparse matrix. The dimensions of the sparse matrix indicate that there are 500 neurons in the physical reservoir, each of which can connect to up to 500 other neurons. The size of the physical reservoir has a significant impact on the feature representation capability, computational complexity, and dynamic behavior of the ESN. As you mentioned, it indeed represents the dimensionality of the reservoir.\n\n2.2 Does the size of the physical reservoir correspond to the array dimensions of the devices? How are these device arrays fabricated, controlled, and measured? I would suggest that the authors provide a description of the operation of the system.\nThank you for your question regarding the relationship between the physical reservoir size and the device array dimensions, as well as the fabrication, control, and measurement of the device arrays. The size of the physical reservoir directly corresponds to the dimensions of the device array for the English letter recognition and Lorenz attractor tasks. We have added a comprehensive description of the array fabrication, control, and measurement procedures, as well as the overall system operation, in the revised manuscript and supplementary materials to provide greater clarity.\nFabrication of the Device and Array\nThe bottom gate electrodes were fabricated via physical vapor deposition. Aluminum oxide was subsequently deposited as the gate dielectric using atomic layer deposition. For the semiconductor layers, ZnO was prepared by dissolving the precursor in ammonia to form a uniform solution, while P3HT was dissolved in toluene and thoroughly mixed by stirring and heating to ensure homogeneity. The substrate was placed on a spin coater, and the ZnO precursor solution was dispensed onto it. Spin-coating was performed sequentially at 1500 rpm, 2000 rpm, and 2500 rpm, each for 30 seconds. After spin-coating, the substrate was annealed at 150 °C for 40 minutes to remove the solvent and crystallize the ZnO thin film. The P3HT solution was then applied onto the ZnO layer, with the volume carefully controlled to ensure film uniformity, followed by spin-coating at 1500 rpm for 30 seconds. The resulting film was annealed at 100 °C for 20 minutes.\nDevice and array structures were defined by covering the thin film with a mask and depositing electrodes via physical vapor deposition. All fabrication steps were conducted in a nitrogen atmosphere. The completed array was encapsulated by covering it with a quartz glass sheet and sealing the edges with UV-curable adhesive.\nControl of the Array\nSpecific rows are selected by applying a feature signal that connects the input voltage to the drain terminals of all transistors in the chosen row via a digital-to-analog converter (DAC). Simultaneously, specific columns are addressed by applying a column selection signal to the gate terminals of all transistors in the selected column. This addressing scheme enables individual control of each unit within the array, facilitating the implementation of matrix multiplication and activation functions.\nMeasurement of the Array\nA custom program was used to sequentially access and test each node in the array, thereby evaluating the operational performance and stability of the array. To measure the output, a resistor with a known value (1 MΩ) was connected in series at the output terminal of the transistor array, allowing the output current to be converted into a voltage signal for further analysis.\n\nComment 3:\nThe writing can be further improved by condensing the manuscript.\n\nResponse to Comment 3:\nThank you very much for your suggestion. We have carefully considered your comments and revised the manuscript accordingly. Several figures and sections of content have been removed from the main text and included in the supplementary materials to improve the conciseness of the article.\n\nFig. S20 The reservoir states of ESNs with reservoir size=800 when recognizing the MNIST dataset. For the digit dataset, the reservoir states of (a) n-type ESN, (b) ambipolar ESN, and (c) p-type ESN. For the fashion dataset, the reservoir states of (d) n-type ESN, (e) ambipolar ESN, and (f) p-type ESN.\nWe illustrate the reservoir states for different devices on both the handwritten and fashion datasets. The reservoir states of n-type ESN on handwritten digit dataset are between 0.93 uA to 0.94 uA and on fashion dataset is 0.94 uA. The reservoir states of ambipolar ESN on handwritten digit dataset are between -0.13 uA and 0.13 uA and on the fashion dataset range from 0.15 uA to -0.17 uA. The reservoir states of p-type ESN on handwritten digit dataset are between 6.33 uA to 12.46 uA and on fashion dataset is between 6.91 µA and 12.15 µA.\nThe reservoir states of the n-type ESN are nearly constant, indicating a lack of variability and ineffective feature conversion. Consequently, the reservoir layer constructed with n-type transistors fails to produce valid features. In contrast, the reservoir states of the ambipolar ESN fluctuate between positive and negative values, whether applied to the digits dataset or the fashion dataset. This fluctuation maintains a relative balance within the range of values, making the ambipolar ESN's reservoir states effective.\nThe reservoir states of the p-type ESN generally reside in the positive region. As the data dimensionality recognized by the ESN increases, these states shift from a bias towards negative values to a bias towards positive values, yet they remain effective for recognition tasks.\n\nFig. S27 (a) Prediction image for typhoon centers by n-type ESN, ambipolar ESN and p-type ESN. (b) NRMSE of typhoon center predictions by n-type ESN, ambipolar ESN and p-type ESN. (c) Predictions of typhoon trajectories by n-type ESN, ambipolar ESN and p-type ESN.\nA dataset was created using videos of Typhoon Kanu from 2023. Grayscale images with a size of 225 × 225 at time t were compressed to 25 × 25 and used as input, with the grayscale image at time t+1 serving as the target. The dataset consists of 251 images, initially used for training and subsequently for evaluation. Fig. S27(a) displays the predicted image of the typhoon center. The NRMSE for the n-type, ambipolar, and p-type ESNs are 12.17%, 10.88%, and 11.73%, respectively. The ambipolar ESN achieves the lowest NRMSE, while the n-type ESN has the highest NRMSE and struggles with high-dimensional data. To further assess performance, the NRMSE of ESNs with varying reservoir sizes was collected. Fig. S27(b) presents the NRMSE for different reservoir sizes, with all three subplots sharing the same color bar. The NRMSE of the n-type ESN is predominantly in the red area, indicating higher error, while the ambipolar ESN's NRMSE is in the blue area, indicating lower error. The p-type ESN's NRMSE falls in the medium range. Overall, the ambipolar ESN demonstrates the lowest NRMSE. We illustrate the reservoir states of ESNs during typhoon image prediction (Fig. S26). For a reservoir size of 200 × 200, the n-type ESN exhibits states ranging from 0.11 µA to 0.93 µA, the ambipolar ESN ranges from 0.15 µA to -0.16 µA, and the p-type ESN ranges from -0.02 µA to 0 µA. When the reservoir size increases to 2000 × 2000, the n-type ESN stabilizes at 0.94 µA, the ambipolar ESN ranges from -0.16 µA to 0.14 µA, and the p-type ESN ranges from 0.01 µA to 0 µA. In multi-dimensional tasks, the n-type ESN initially shows positive values, but at a size of 2000 × 2000, it becomes constant and loses activity. The ambipolar ESN maintains a balance between positive and negative values, ensuring diversity and activity in its states. The p-type ESN shows minimal variation in its states.\nThe trajectory of Typhoon Kanu is predicted using a dataset comprising longitude and latitude data. The location at time step t serves as the input, while the location at t+2 is the target. Data from steps 1 to 50 are used for training, and steps 51 to 90 are used for testing. Fig. S27(c) displays the prediction results, indicating that with a reservoir size of 500 × 500, the ambipolar ESN most accurately matches the real trajectory on the test set.\nThe NRMSE for different reservoir sizes is also analyzed (Fig. S28). For longitude prediction, the n-type ESN's NRMSE ranges from 4.93% to 17.74%, the ambipolar ESN ranges from 6.99% to 16.93%, and the p-type ESN is approximately 9.31%. For latitude prediction, the n-type ESN's NRMSE ranges from 13.90% to 16.56%, the ambipolar ESN ranges from 8.34% to 8.46%, and the p-type ESN ranges from 13.49% to 56.24%. The n-type ESN achieves the lowest NRMSE for longitude but performs poorly in latitude prediction. The ambipolar ESN, while slightly less effective than the n-type ESN for longitude, excels in latitude prediction, especially at a reservoir size of 500 × 500, where it achieves a reasonable NRMSE. Thus, predictions at this reservoir size are the most accurate (Fig. S29). We further illustrate that at a reservoir size of 500 × 500, the ambipolar ESN demonstrates superior performance.\n\nComment 4:\nThe color bar is missing in several figures, such as Figure 2.\n\nResponse to Comment 4:\nThank you very much for your comment. We have added color bars to the relevant figures, including Figure 2, 3 and 5, as suggested.\n\nFig. 2 An ESN constructed with transistors for letter recognition. c Reservoir state output of an ESN. d Confusion matrix of ESN recognizing English letters.\nFig. 3 Physical ESN used for MNIST dataset recognition. a Confusion matrices of n-type ESN, ambipolar ESN and p-type ESN on the digits dataset. b Confusion matrix of n-type ESN, ambipolar ESN and p-type ESN on the fashion dataset.\nFig. 5 Multimodal identification task by our ESN. c Reservoir state of the ECG image and the ECG signal.\n\nComment 5:\nThe text in Figure 2a is horizontal, which seems unnecessary.\n\nResponse to Comment5:\nThank you very much for your suggestions. We have made revisions in accordance with your suggestions.\n\nFig. 2 An ESN constructed with transistors for letter recognition. a The transformation process of an English letter at the output layer.\n\nResponse to Reviewer 3:\nThank you very much for your valuable comments and constructive feedback. We have carefully reviewed your concerns and have made substantial revisions to the manuscript to address your suggestions. In particular, we have worked to improve the clarity, logical flow, and language throughout the text. We have also revised the presentation of our results and further clarified the novelty and significance of our work. We hope that these changes have addressed your concerns and improved the overall quality of the manuscript. Thank you again for your thoughtful review.\n\nComment 1:\nThe manuscript lacks comprehensive device testing. The authors need to provide data on endurance, retention, uniformity, switching speed, and environmental stability to evaluate whether the device is sufficiently robust and suitable for computing applications.\n\nResponse to Comment1:\nThank you very much for your insightful comment. In response, we have conducted comprehensive performance tests on the device, including evaluations of endurance, retention, uniformity, switching speed, and environmental stability. The results of these tests have been added to the revised manuscript and supplementary materials to demonstrate the robustness and suitability of our device for computing applications.\n\nFig. S2 Endurance characteristics of the device.\nWe performed cyclic endurance testing on the device, subjecting it to 10⁵ continuous read-write cycles. After 10⁵ cycles, key performance metrics, such as the current on/off ratio, exhibited no significant degradation. The specific data are presented in Fig. S2, which demonstrates that the device maintains stable switching characteristics even after extensive cycling. These results indicate that the device possesses excellent endurance, can reliably withstand repeated read-write operations, and is well-suited for computing applications.\n\nFig. S3 Retention characteristics of the device.\nWe conducted long-term retention tests on the devices at room temperature. The results show that the stored information exhibited no significant leakage or degradation over a period of 10⁴ seconds. These findings demonstrate that the device possesses excellent retention characteristics, effectively preserving stored information during extended power-off periods. This makes the device highly suitable for non-volatile memory applications.\n\nFig. S4 The uniformity of the device array. (a) Low conductance state of the array. (b) High conductance state of the array.\nWe performed tests on multiple devices within a 16×16 physical array, measuring parameters such as the low and high conductance states. The results show that the standard deviations for both the low and high conductance states are very low, demonstrating excellent device uniformity across the array.\n\nFig. S5 The switch speed of the device. (a) From high conductance state to low conductance state. (b) From low conductance state to high conductance state.\nWe evaluated the switching speed of the device by performing switching operations and measuring the corresponding switching time. The device exhibited a switching time of less than 66 ns, demonstrating a rapid response to pulse signals. This fast switching speed meets the requirements for high-speed computational applications, further confirming the device's suitability for such uses.\n\nFig. S6 Environmental stability of the device.\nTo address the potential oxidation and moisture absorption of P3HT in air, the device was encapsulated and tested under ambient conditions. As shown in Fig. S6, the transfer characteristic hysteresis loops of the device in air remain nearly unchanged over a period of seven days. This demonstrates the device's excellent environmental stability.\n\nComment 2:\nThe authors mention in the introduction that transistor sizes are approaching physical limits. What is the actual size of the proposed device? Does this new design contribute to reducing transistor dimensions or increasing the density of transistors per unit area?\n\nResponse to Comment2:\n2.1 The authors mention in the introduction that transistor sizes are approaching physical limits. What is the actual size of the proposed device?\nThank you very much for your interest in our research. The device array features a channel length of 50 microns, and we have also fabricated devices with a channel length of 2 microns. While these dimensions are larger than those achieved by current silicon-based transistors, our device offers unique advantages for neuromorphic computing applications due to its distinct operational mechanisms. Fig. S9(a) presents the surface morphology of the device with a 2-micron channel length. Additionally, the corresponding transfer and output characteristic curves for this transistor are also provided.\n\nFig. S9 2-micron channel device performance. (a) Transfer curve under negative V . (b) Transfer curve under positive V . (c) Output curve of the p-type region, and (d) Output curve of the n-type region.\nFig. S9 (a) and (b) show the transfer curves of the device. As V increases from its initial value, I begins to change. In the n-type conduction region, when V exceeds the threshold voltage (V ), the device operates predominantly through electron transport, resulting in a significant increase in Ids with increasing V . As V continues to increase and surpasses another threshold voltage (V ), the device enters the p-type conduction region, where hole transport becomes dominant and I increases again. Fig. S9(c) and S9(d) present the output characteristics of the device. Due to the ambipolar nature of the transistor, which allows for the simultaneous transport of both electrons and holes, the output curves exhibit bidirectional current transport characteristics. Specifically, the output current initially increases, then decreases, and subsequently increases again. This behavior arises from the interplay between electron and hole transport under varying applied voltages.4,5\nReference:\n4. Stein, E., Nahor, O., Stolov, M. et al. Ambipolar blend-based organic electrochemical transistors and inverters. Nat. Commun. 13, 5548 (2022).\n5. On, S., Kim, Y.-J., Lee, H.-K. & Yoo, H. Ambipolar and anti-ambipolar thin-film transistors from edge-on small-molecule heterostructures. Appl. Surf. Sci. 542, 148616 (2021).\n\n2.2 Does this new design contribute to reducing transistor dimensions or increasing the density of transistors per unit area?\nIn most neuromorphic computing systems, matrix calculations typically require a pair of devices for each node. One device handles long-term potentiation and the other handles long-term depression, allowing for the computation of both positive and negative values by taking the difference between them. This approach effectively doubles the number of devices needed for matrix computation tasks.R6–R8 In contrast, our device exhibits both n-type and p-type characteristics, enabling simultaneous computation of positive and negative values within a single device. As a result, our design requires only half the number of devices compared to conventional approaches. Furthermore, due to the sparsity and specificity of ESN, often only 5%–10% of the nodes are actively involved in computation. Therefore, to achieve matrix calculations of the same dimensionality, our design theoretically requires only 2.5%–5% of the total device count, significantly increasing the integration density and reducing the overall footprint.\nReference:\nR6. Kudithipudi, D., Schuman, C., Vineyard, C.M. et al. Neuromorphic computing at scale. Nature 637, 801–812 (2025).\nR7. van de Burgt, Y., Melianas, A., Keene, S.T. et al. Organic electronics for neuromorphic computing. Nat. Electron. 1, 386–397 (2018).\nR8. Tang, J. et al. Bridging Biological and Artificial Neural Networks with Emerging Neuromorphic Devices: Fundamentals, Progress, and Challenges. Adv. Mater. 31, 1902761 (2019).\n\nComment 3:\nReservoir computing typically requires specific temporal dynamics. Do the proposed devices exhibit any intrinsic temporal properties essential for reservoir computing?\n\nResponse to Comment3:\nThank you very much for your insightful comment. Indeed, our device exhibits intrinsic temporal dynamics that are essential for reservoir computing. We conducted detailed measurements of the device’s response under various input signals. The experimental results demonstrate that the device displays both time delay and memory effects, which are key temporal properties required for reservoir computing. Furthermore, by adjusting the frequency and amplitude of the input signals, we can modulate the device’s dynamic response, thereby optimizing its performance in reservoir computing tasks. These findings confirm that our device is well-suited for applications requiring specific temporal dynamics.\n\nFig. S11 (a) The paired-pulse facilitation performance of the device, (b) PPF index of the device under different voltages, (c) 6-bit pulse encoding of the device (V =+5V pulse width 10 ms), (d) potential and depression performance of the device.\nPaired-pulse facilitation (PPF) is a dynamic response characteristic of devices and an important manifestation of short-term plasticity. The calculation formula for the PPF index is as follows:30,31\nPFF = A /A ×100 %\n2 1\nWhere A and A represent the amplitude of excitatory postsynaptic conductance under two adjacent input pulses, respectively. As shown in Fig. S11(a), when pair +5V pulses with a width of 10 ms are applied to the device, a progressive relationship appears in I . Fig. S11(b) illustrates the variation of the PPF index with the interpulse interval (Δt) between adjacent input pulses. As Δt increases, the PPF index exhibits an exponential decline. This phenomenon resembles the characteristics of biological synapses and can be well fitted by a double exponential function32:\nPPF = C ×exp(−∆t/τ )+C ×exp(−∆t/τ )\n1 1 2 2\nHere, C and C represent the initial facilitation intensities of the fast and slow phases respectively, while τ and τ correspond to the relaxation times of the two phases.\n1 2\nFig. S11(c) illustrates the 6-bit pulse encoding of the device. By precisely controlling the pulse parameters, the input signal is encoded into 64 distinct states, each corresponding to a unique dynamic response pattern. The dynamic mechanism of the physical reservoir enables the system to efficiently accomplish complex temporal signal processing tasks.\nFig. S11(d) shows that the synaptic weights are modulated by continuous pulses of different quantities of +5V gate voltage, which is similar to the phenomenon of neurotransmitter increase under constant external stimulation in biological synapses. The increase in the number of pulses leads to a rise in the injected electron concentration, thereby causing conductance enhancement. This type of modulation of synaptic weights induced by V pulses belongs to potential performance. When no gate voltage is applied afterward, the electron pair concentration gradually decreases, and the device exhibits a depression state.33,34\nReference:\n30. Yang, J.-T. et al. Artificial Synapses Emulated by an Electrolyte-Gated Tungsten-Oxide Transistor. Adv. Mater. 30, 1801548 (2018).\n31. Zucker, R. S. & Regehr, W. G. Short-Term Synaptic Plasticity. Annu. Rev. Physiol. 64, 355-405 (2002).\n32. He, H.-K. et al. Photonic Potentiation and Electric Habituation in Ultrathin Memristive Synapses Based on Monolayer MoS . Small 14, 1800079 (2018).\n2\n33. Mou, X. et al. Analog memristive synapse based on topotactic phase transition for high-performance neuromorphic computing and neural network pruning. Sci. Adv. 7, eabh0648.\n34. Liu, Z., Zhang, Q., Xie, D. et al. Interface-type tunable oxygen ion dynamics for physical reservoir computing. Nat. Commun. 14, 7176 (2023).\n\nComment 4:\nCan the operating voltage of this device be reduced in future implementations? The authors should provide an analysis of potential strategies for lowering power consumption.\n\nResponse to Comment 4:\n4.1 Can the operating voltage of this device be reduced in future implementations?\nYes, the operating voltage of the device can be further reduced in future implementations. By employing thinner aluminum oxide dielectric layers, such as 30 nm or even 10 nm, the operating voltage can be significantly lowered. As demonstrated in Fig. S7, the output and transfer curves of the device with a 10 nm aluminum oxide dielectric layer clearly show reduced working voltages. This approach offers a promising pathway for achieving low-power operation in future device designs.\n\nFig. S7 Transfer characteristics of devices with different dielectric layers. (a) Transfer characteristics with 30 nm aluminum oxide, (b) Transfer characteristics with 10 nm aluminum oxide.\nFig. S7 illustrates the transfer characteristics of the device. It can be observed that the operating voltage is strongly influenced by the thickness of the dielectric layer. As the thickness of the Al₂O₃ dielectric layer decreases, the electric field intensity increases significantly. This enhanced field enables more effective modulation of the charge distribution at the semiconductor surface, allowing the same level of charge accumulation or depletion to be achieved at lower voltages.\n\nFig. S8 Output characteristics of n-type, ambipolar and p-type transistors with 100 nm Al O . (a) output characteristics of n-type transistor in negative region, (b) output characteristics of ambipolar transistor in negative region, (c) output characteristics of p-type transistor in negative region, (d) output characteristics of n-type transistor in positive region, (e) output characteristics of ambipolar transistor in positive region, (f) output characteristics of p-type transistor in positive region.\n\nFig. S10 The transistor performance with 10 nm Al O . (a) Transfer curve under negative V , (b) Transfer curve under positive V , (c) Power consumption under negative V , (d) Output curve under positive V , (e) Output curve under negative V , (f) Power consumption under positive V .\nAmbipolar charge transport phenomena are observed at an operating voltage of 3 V, as evidenced by the presence of both p-type and n-type operation modes in the measured transfer I–V characteristic curves (Fig. S10(a) and (b)). Specifically, under negative V bias, hole transport dominates and enhances the current, while under positive V bias, electron transport becomes predominant, resulting in increased current. This behavior produces characteristic V-shaped curves for V under both positive and negative biases.\nFig. S10(d) and (e) show the output curves of the ambipolar transistor in p-type and n-type operation modes (i.e., V < 0 and V > 0, respectively). In these output curves, a linear increase in I is observed at lower Vgs biases, indicating that the majority carrier concentration (holes in p-type mode or electrons in n-type mode) increases with rising |V |. This phenomenon is consistent with previous reports.\nFig. S10(c) and (f) illustrate the device's power consumption. The power consumption is calculated as P = |V | × |I | × t , where V and I are the absolute values of the gate-source voltage and drain-source current, respectively, and t is the pulse duration of a single pulse. The power consumption for positive and negative spikes is 812 fJ and 735 fJ, respectively.\n\n4.2 The authors should provide an analysis of potential strategies for lowering power consumption.\nTo highlight the potential of our device for low-power applications, we have compiled the power consumption data of various neuromorphic devices reported in recent years. Table S1 summarizes the power consumption of these neuromorphic computing devices. It can be seen that, by reducing the operating voltage, our device achieves and maintains a relatively low power consumption level, comparable to or even lower than that of other state-of-the-art devices.\n\nTable S1. Power consumption of neuromorphic devices.\nMaterials Type Power Dataset acc Ref.\nMoS h-BN transistor ≈450 pJ MNIST 91.1 % 6\n2/\nMoS /Nb O transistor ≈6 pJ MNIST 94.2% 7\n2 2 5\nMoSe transistor ≈140 nJ MNIST ~70 % 8\n2\nα-In Se transistor ≈4 nJ MNIST 91.9 % 9\n2 3\n180nm CMOS transistor ≈77.24pJ 10\nIGZO transistor 0.15 nJ MNIST 73% 11\nCsPbBr QDs transistor 1nJ MNIST 73% 12\nIGZO/HfZrOx transistor 5.69nJ MNIST 91.1% 14\nAlGaN/GaN transistor 18.5 pJ 15\nSi doped-GaN transistor 33.4 pJ 16\nGa O /ZnO heterojunction 236 pJ MNIST 92.31% 17\n2 3\nPTCDA/MoS transistor 10 pJ 18\nGraphene/SWNTs phototransistor 2.5nJ 20\npentacene/PMMA/CsPb\nflash memory 1.4 nJ 21\nBr QD\nITZO phototransistor 2.5 nJ 90.3% 23\nPorphyrin/Sinanowires phototransistor 1 nJ 24\nP3HT/CsPbBr QD phototransistor 250 nJ 25\nx\nTFT-CN/DTT-8 phototransistor 5.3nJ 27\nFe O /WS memristor 72 fJ 94.3% 28\n2 3 2\nTiN transistor 30fJ 29\nThis\nP3HT/ZnO transistor 812 fJ MNIST 95.27%\nwork\nReference:\n6. Oh, S., Lee, J.J., Seo, S., Yoo, G. & Park, J.H. Photoelectroactive Artificial Synapse and Its Application to Biosignal Pattern Recognition. npj 2D Mater. Appl. 5 (1), 95 (2021).\n7. Nam, J. H. et al. Low Power MoS /Nb O Memtransistor Device with Highly Reliable\n2 2 5\nHeterosynaptic Plasticity. Adv. Funct. Mater. 31, 2104174 (2021).\n8. Wang, Y. et al. Near-Infrared-Irradiation-Mediated Synaptic Behavior from Tunable Charge-Trapping Dynamics.Adv. Electron. Mater. 6, 1900765 (2020).\n9. Wang, L. et al. Exploring Ferroelectric Switching in α-In2Se3 for Neuromorphic Computing. Adv. Funct. Mater. 30, 2004609 (2020).\n10. Kim, J., Kim, J. N., Kim, Y., Hwang, S. & Koo, M. Design of a 180 nm CMOS Neuron Circuit with Soft-Reset and Underflow Allowing for Loss-Less Hardware Spiking Neural Networks. Adv. Intell. Syst. 6, 2300460 (2024).\n11. Sun, J. et al. Optoelectronic Synapse Based on IGZO-Alkylated Graphene Oxide Hybrid Structure. Adv. Funct. Mater. 28, 1804397 (2018).\n12. Lv, Z. et al. Mimicking Neuroplasticity in a Hybrid Biopolymer Transistor by Dual Modes Modulation. Adv. Funct. Mater. 29, 1902374 (2019).\n13. Yu, R. et al. Electret-Based Organic Synaptic Transistor for Neuromorphic Computing. ACS Appl. Mater. Interfaces 12, 15446-15455 (2020).\n14. Kim, M.-K. & Lee, J.-S. Ferroelectric Analog Synaptic Transistors. Nano Lett. 19, 2044-2050 (2019).\n15. Kai, C. et al. AlGaN/GaN-Based Optoelectronic Synaptic Devices for Neuromorphic Computing. Adv. Opt. Mater. 11, 2202105 (2023).\n16. Mo, Y., Luo, B., Dong, H. & Hou, B. Light-stimulated artificial synapses based on Si-doped GaN thin films. J. Mater. Chem. C 10, 13099-13106 (2022).\n17. Wang, J. et al. Piezo-phototronic effect modulated optoelectronic artificial synapse based on a-Ga O /ZnO heterojunction. Nano Energy 120, 109128 (2024).\n2 3\n18. Wang, S. et al. A MoS /PTCDA Hybrid Heterojunction Synapse with Efficient\n19. Wang, Y. et al. Silicon-Nanomembrane-Based Broadband Synaptic Phototransistors for Neuromorphic Vision. Nano Lett. 23, 8460-8467 (2023).\n20. Qin. S., Wang F, Liu Y. et al. A light-stimulated synaptic device based on graphene hybrid phototransistor. 2D Mater, 4, 035022 (2017).\n21. Wang, Y. et al. Photonic Synapses Based on Inorganic Perovskite Quantum Dots for Neuromorphic Computing. Adv. Mater. 30, 1802883 (2018).\n22. Yin, L. et al. Synaptic silicon-nanocrystal phototransistors for neuromorphic computing. Nano Energy 63, 103859 (2019).\n23. Kim, J. et al. Light-stimulated artificial photonic synapses based on solution-processed In-Sn-Zn-O transistors for neuromorphic applications. J. Alloy. Compd. 903, 163873 (2022).\n24. Li, X. et al. Complementary Photo-Synapses Based on Light-Stimulated Porphyrin-Coated Silicon Nanowires Field-Effect Transistors (LPSNFET). Small 17, 2101434 (2021).\n25. Gupta, G. K., Kim, I.-J., Park, Y., Kim, M.-K. & Lee, J.-S. Inorganic Perovskite Quantum Dot-Mediated Photonic Multimodal Synapse. ACS Appl. Mater. Interfaces 15, 18055-18064 (2023).\n26. Mi, Y.-C., Yang, C.-H., Shih, L.-C. & Chen, J.-S. All-Optical-Controlled Excitatory and Inhibitory Synaptic Signaling through Bipolar Photoresponse of an Oxide-Based Phototransistor. Adv. Opt. Mater. 11, 2300089 (2023).\n27. Zhu, X. et al. High-Contrast Bidirectional Optoelectronic Synapses based on 2D Molecular Crystal Heterojunctions for Motion Detection. Adv. Mater. 35, 2301468 (2023).\n28. Ghafoor, F. et al. Novel Solution-Processed Fe O /WS Hybrid Nanocomposite\n2 3 2\nDynamic Memristor for Advanced Power Efficiency in Neuromorphic Computing. Adv. Sci. 12, 2408133 (2025).\n29. Xu, H., Shang, D., Luo, Q. et al. A low-power vertical dual-gate neurotransistor with short-term memory for high energy-efficient neuromorphic computing. Nat. Commun. 14, 6385 (2023).\n\nComment 5:\nIt is unclear what hardware is used for the later application demonstrations. If the applications are performed through simulations, the authors should explicitly state this to help readers better assess the significance of the work.\n\nResponse to Comment 5:\nThank you very much for your valuable comment. For the application demonstrations, we utilized a 16×16 hardware array to perform both English letter recognition and Lorenz attractor prediction tasks. The hardware setup is shown in Fig. 1f. All results presented in these sections were obtained from experiments conducted on this physical hardware array, rather than through simulations. We believe this hardware-based demonstration further highlights the practical significance and potential of our device for real-world neuromorphic computing applications.\n\nFig. 1 Structure and characterization of the ambipolar transistor. a Schematic diagram of ESN. b TEM image of the cross-sectional morphologies under 100,000x magnification: n-type transistor, amipolar transistor and p-type transistor cross-section. High-resolution cross-sectional morphologies of ZnO layers in n-type transistor, ambipolar transistor and p-type transistor. c Surface morphologies observed under AFM: ZnO layers in n-type, ambipolar, and p-type transistors, as well as P3HT layers in n-type, ambipolar, and p-type transistors. d Output curves of n-type, ambipolar, and p-type transistors using 100 nm SiO as the dielectric layer. e Optical image of the ambipolar transistor array. f Optical image of the physical ESN system. g Output characteristic of ambipolar transistor in the array and tanh function.\nOther demonstrations were carried out using computer simulations. We have added the following explanation to the experimental section:\nThe English alphabet recognition task and Lorenz attractor prediction (16×16) were performed using a 16×16 ambipolar transistor hardware array, with read/write operations managed by DACs and an ADC controlled via STM32. In contrast, tasks such as MNIST handwritten digit recognition, Fashion dataset classification, Lorenz attractor prediction, ECG data classification, and multimodal recognition with the MOSEI dataset were implemented through simulations using Python and PyTorch. These simulations were based on an operational model of the device. The specific physical ESN network structure used in these experiments is detailed in Table S2. To facilitate evaluation, we have also uploaded the code for the MNIST handwritten digit recognition and Lorenz curve prediction code simulation portion. You can find it at the following URL:\nhttps://github.com/WenMinZhong/physical_ESN\n\nComment 6:\nHow are the input pixels translated into time-series sequences in Figure 2? A more detailed explanation is necessary.\n\nResponse to Comment6:\nThank you for your comment. We have provided additional clarification below.\nTo convert the input pixels into time-series sequences, we first flatten each 5×5 two-dimensional image into a one-dimensional vector containing 25 elements. For the entire training set, these vectors are then arranged sequentially to form a time-series dataset, progressing from the first character to the last. At each time step t, the system receives a 25-dimensional input vector corresponding to one image. This process is repeated for each image in the dataset, so that the time-series sequence is constructed by feeding the images one after another as input vectors over time. The specific operational steps are illustrated in Fig. S13.\n\nFig. S13 Schematic diagram of converting English letters into time series.\n\nComment 7:\nThe authors appear to read out the current at each time step, suggesting that all temporal information originates from the input signal. However, the transistor array itself does not seem to perform any nonlinear transformation in the temporal domain. Can the authors clarify this?\n\nResponse to Comment 7:\nThank you for your review comments. When employing hardware, the dynamic mechanism of the device is utilized. As a computational system for time series, we have also taken into account its dynamic response and nonlinear transformation behavior.\nFor time-varying data with continuous temporal attributes, the ESN requires leaky integration units. This leakage unit can accomplish its task by utilizing the dynamic response behavior of previous devices.\nx(n+1) = (1−δCp)∗x(n)+δCf(W u(n+1)+W ∗x(n))\nin res\nwhere C is a time constant and p is the leaking decay rate. δ is time step. We implemented this process using a 1ms pulse width and time step, thus allowing C and δ to be set to 1. The parameter a originates from the device's forgetting characteristics. We demonstrate the forgetting characteristics of the device under a 5V pulse. Its relaxation factor τ is 2.56 ms, as shown in Fig. S12(a). Fig. S12(b) displays its a-value when functioning as a dynamic neuron. Through this method, dynamic responses are achieved for nonlinear mapping in ESN.\n\nFig. S12 (a) Dynamic response of the device. (b) Leakage rate of the device at different times.\n\nComment 8:\nIf the prediction is accurate and the dataset is uniform, why do the confusion matrices exhibit different color distributions? The authors should explain these discrepancies.\n\nResponse to Comment 8:\nThank you very much for your comment, and we apologize for any confusion caused by this issue. Ideally, a dataset should maintain uniformity. However, in our case, after customizing the dataset, we randomly selected two thousand samples as the test set. This random selection introduced a certain degree of non-uniformity in the sample distribution across different classes, which is reflected in the varying color intensities in the confusion matrices. Despite this imbalance, the prediction results remain accurate. As shown in the original confusion matrix results, all predictions fall along the diagonal, indicating 100% accuracy. The differences in color distribution are solely due to the unequal number of samples per class in the test set, not due to misclassification. This is clearly illustrated in the figure.\n\nFig. 2d The confusion matrix for English letter recognition, as previously mentioned in the article content.\nFurthermore, we reconstructed a uniform dataset and found that the accuracy remained at 100%. Therefore, we conclude that, although our original dataset exhibited some variability, this had minimal impact on the final results.\n\nFig. S17 The confusion matrix for English letter recognition under a uniform test set.\n\nComment 9:\nThe MNIST and Fashion MNIST datasets are inherently static. How did the authors convert these static datasets into temporal voltage signals for processing within the reservoir computing framework?\n\nResponse to Comment 9:\nThank you very much for your comments. In our paper, we adopted a one-time processing method. This involves flattening the 60,000 training samples. Thus, a 28*28 image is transformed into a one-dimensional matrix with 784 elements. The 60,000 samples are then converted into 60,000 time steps and fed into the system sequentially along the temporal dimension. In this way, each sample undergoes dimension transformation by W and is then converted into temporal voltage input to the system.\nin\nTo better illustrate this process, we provide a schematic diagram of the specific operations as shown in Fig. S19.\n\nFig. S19 Time series schematic diagram of the MNIST handwritten digit dataset.\nFor example, consider selecting 10 samples from the MNIST dataset. Each 28×28 image is first flattened into a 784-dimensional vector, resulting in a 10×784 matrix, where 10 corresponds to the number of time steps. If the sample size is increased to 60,000, the time step length becomes 60,000, and the ESN can then output 60,000 reservoir states accordingly.\n\nComment 10:\nThe authors need to demonstrate that the reservoir computing component is indeed contributing to the classification performance. Additionally, the classification network used should be clearly described. There is a concern that the high accuracy is primarily attributed to the readout network rather than the reservoir itself.\n\nResponse to Comment 10:\nThank you very much for your comments. This device outputs the reservoir space state, which indeed can enhance classification performance. Our system features a hidden layer constructed with transistors that has recurrent connections, enabling it to retain dynamic information from input data. For the MNIST dataset, although the images themselves are static, they can be transformed into sequential inputs for the reservoir by unfolding the image pixels row-wise or column-wise. The dynamic memory capability of the reservoir helps capture spatial relationships between pixels, such as stroke continuity and local shape features. This memory function allows the reservoir to better comprehend the overall structure of the images, thereby improving recognition accuracy.\nThe readout layer for classification is a matrix of 1785×10. Its training method employs the traditional linear regression calculation approach rather than the gradient descent method (the manuscript's experimental section provides the training methodology and relevant formulas).\nWe are indeed concerned about the issues you have raised. As you mentioned, perhaps our classification accuracy is solely attributed to the readout layer. Therefore, we have carefully and cautiously examined this matter. However, upon thorough investigation, we found that if only the contribution of the readout layer is considered, the accuracy on the MNIST handwritten digit dataset would be merely 85.46%, and the accuracy on the Fashion dataset would only reach 81.24%. In contrast, the accuracy reported in our manuscript is higher than these values. Thus, we believe our system is capable of enhancing classification performance.\n\nFig. S22 The accuracy rate of recognizing handwritten digits and fashion datasets only using the readout layer.\nTo verify that the accuracy does not solely originate from the readout layer, we provide classification tasks using only the readout layer. The accuracy rates on the digits and fashion datasets are 85.46% and 81.24%, respectively, which are lower than our results (Fig. S22)."}, {"round": 2, "reviewer_comments": {"reviewer_1": "I appreciate the authors' thorough revision of the manuscript. However, after reading the revised manuscript as well as the supplementary information, I still cannot obtain a clear understanding of how the suggested physical reservoir computing approach works. For example, for the English letter classification (the first example) task, Table S2 states the input layer dimension is (25 x 16), which, in my opinion, should be 25 x 1. In Fig. S19, I do not understand how different numbers are combined to form a single time series (t1-t10) input, where the output cannot be defined.\nWhile I agree that the topic of physical reservoir computing is highly important and the use of ambipolar transistors for this purpose can benefit researchers in the field, the current manuscript does not provide clear illustrations and details of the experimental procedure.", "reviewer_2": "The authors have addressed all my comments, and I recommend the manuscript to be accepted for publication.", "reviewer_3": "The authors have addressed most of my concerns. However, several aspects of the experimental methodology and descriptions remain unclear. Please see the detailed comments below:\n1. The label in Fig. S5 appears to be incorrect.\n2. The authors have successfully fabricated a thinner Al₂O₃ layer that enables lower operating voltage. Given this advancement, it is unclear why this improved device was not used in the main experiments. Are there any potential limitations? Clarifying this point would be highly valuable for guiding future research in the field.\n3. While the Fig. S5 states that the device operates at nanosecond-level switching speeds, the experiments were performed using microsecond-long pulses. This discrepancy should be addressed.\n4. The calculation method for power consumption remains unclear. The reported values cannot be reproduced from the data provided. The authors should include a step-by-step description or an example calculation.\n5. The working flow of the network described in Fig. S13 is still not clearly explained. The input image is flattened and directly passed to an array, raising the question of why a reservoir is necessary. A simple linear layer could classify the inputs in a single step. The authors should further justify the use of the reservoir architecture.\n6. It is unclear why the authors treat different samples in MNIST as different time steps. This approach is confusing. Do the authors intend for the reservoir to extract temporal features across samples? Further clarification and justification are needed.\n7. Fig. S22 does not clearly demonstrate the effectiveness of the reservoir. It is also unclear whether the dimension of the linear layer was modified. If not, how were the dimensional mismatches handled? The authors should provide a more detailed explanation of the setup and clarify how the comparison was made."}, "author_response": "Reviewer #1 (Remarks to the Author):\nI appreciate the authors' thorough revision of the manuscript. However, after reading the revised manuscript as well as the supplementary information, I still cannot obtain a clear understanding of how the suggested physical reservoir computing approach works. For example, for the English letter classification (the first example) task, Table S2 states the input layer dimension is (25 x 16), which, in my opinion, should be 25 x 1. In Fig. S19, I do not understand how different numbers are combined to form a single time series (t1-t10) input, where the output cannot be defined.\nQ1 For example, for the English letter classification (the first example) task, Table S2 states the input layer dimension is (25 x 16), which, in my opinion, should be 25 x 1.\nResponse:\nThank you very much for your careful review and for your valuable feedback.\nWe sincerely appreciate your recognition of our efforts in revising the manuscript. We apologize for any remaining confusion regarding the physical reservoir computing approach, especially in relation to the English letter classification task.\nWe have revised the manuscript by adding additional content and dividing the recognition of English letters into four distinct regions to enhance clarity and facilitate understanding. Below, we provide detailed explanations for each of these four regions.\nFig. S13 Schematic diagram of converting English letters into time series.\nThe first part of our approach involves flattening a 5×5 image into a 25×1 vector. Following the sequence of the dataset, we assign each letter an order along the time scale: the first letter is denoted as t1, the second as t2, and so forth. We understand that reviewers may have concerns regarding the treatment of image data as a time series. To address this, we have consulted relevant literature that supports this methodology. Specifically, some studies utilize the number of samples in the entire dataset as the time scale, while others divide a single sample into multiple segments and concatenate them to form a time sequence. For example, related research by Zhou et al.,[R1] Sun et al.,[R2] and Han et al.[R3] has demonstrated the application of ESN to image processing tasks using similar strategies. In our study, we adopt the most straightforward approach by directly flattening the image as the input to the model.\n. [R1] Zhou, M., Lun, S. & Li, M. Rapid training echo convolution network for image recognition. Information Sciences 696, 121750, (2025).\n. [R2] Sun, J., Li, L. & Peng, H. in 2021 International Conference on Neuromorphic Computing (ICNC). 165-170.\n. [R3] Han, M., Zhao, S., Yin, H., Hu, G. & Ghadimi, N. Timely detection of skin cancer: An AI-based approach on the basis of the integration of Echo State Network and adapted Seasons Optimization Algorithm. Biomedical Signal Processing and Control 94, 106324, (2024).\nWe appreciate your insightful question regarding the input layer, and would like to provide further explanation as to why it is configured as 25 × 16 instead of 25 × 1. This design choice is dictated by the architecture of ESN, as illustrated in the provided diagram, which depicts a linear transformation layer. Mathematically, this layer functions as a fixed matrix that transforms the dimensionality of the input data, either expanding or reducing it, to ensure compatibility with our physical ESN. This approach is discussed in detail in the study by Sun et al., where the corresponding architecture is also presented (please see the following information). Therefore, in order to match the dimensions of the hardware array (16 × 16), the input layer must be set to 25 × 16.\n“The architecture of ESN is not complicated, it consists of three parts, called input layer, intermediate reservoir, and output layer. Fig. 1 shows that there are I neurons in the input part and E neurons set up the output layer in Echo State Network, that is to say a(i) ∈ RIa and y(i) ∈ RKg respectively. Here i = 1, 2, 3 · · · , i is the discrete time, is the amount of data points in the training datase. In the intermediate reservoir, there are E × E neurons and b(i) ∈ REb means the state vector of the E ×\nb b b\nE matrix. The procedures of feeding input data a(i) into the ESN and updating the internal state b(i) in the reservoir can be seen as inputting data a(i) from low-dimensional space to highdimensional space b(i).”\nFig. 1. Echo State Network and its structure, input layer, reservoir, and output layer.\nThe numbers of their neurons are I, E, K. Its output weight W is the only need to be\nout\ntrained (Sun et al.).\nIn order to match the dimension of the reservoir, it will not be 25×1 unless the reservoir size is 1×1.\nThe third part pertains to the operation of our hardware, which leverages the nonlinearity of the ambipolar transistor array to process data after dimensionality reduction, thereby achieving high-dimensional nonlinear transformation. The operational logic of this section is illustrated in Figure S18. Following processing by the hardware array, a nonlinear output representing the reservoir state is generated. Given that the hardware array is configured as 16 × 16, the dimensionality of its output is 16.\nThe fourth part concerns the operation logic of the readout layer. In this stage, the original data, the reservoir state, and a bias term of 1 are concatenated and subsequently used for training via linear regression. The underlying principle is to utilize the high-dimensional state features produced by the hardware array to enhance the nonlinearity of the data, thereby improving recognizability. To accommodate these inputs, the dimensionality of the readout layer is set to the sum of the original data dimension, the reservoir state dimension, and the bias term. As a result, the readout layer has dimensions of (25 + 16 + 1) × 6. The detailed architecture is provided in Table S2.\n\nQ2 In Fig. S19, I do not understand how different numbers are combined to form a single time series (t1-t10) input, where the output cannot be defined.\nResponse:\nWe apologize for any confusion our previous explanation may have caused. The application of ESN for image recognition on the MNIST dataset is supported by existing literature. In particular, we refer to the work reported by An et al.,[R4] which provides a relevant foundation for our approach.\nFig. 4. Architecture of the proposed deep-DFR (An et al.).\nThey described it in this way “As demonstrated in Fig. 4, the deep-DFR structure for image classification first deploys a fully connected layer as the W layer\nin\nwith dimensions [784, 800] or [3072, 800] for MNIST and CIFAR10 dataset respectively. ” “Then the W layer is built as a fully connected layer with dimensions\nout\n[800, 10].”\nAs illustrated, the first dimension of the input layer (W ) corresponds to the\nin\nimage dimension, while the second dimension represents the reservoir dimension. We would like to note that, as described in our manuscript, the input dimensions for W\nin\nare consistent with those reported in the referenced work when the reservoir size is set to 800.\n[R4] An, Q., Bai, K., Liu, L., Shen, F. & Yi, Y. C. A unified information perceptron using deep reservoir computing. Comput. Electr. Eng. 85, 106705 (2020).\nWith regard to the temporal attributes of the MNIST dataset, we refer to the article by Sun et al. for clarification. In their work, they describe it as follows:\n\"Different from processing sequence problems, when ESN is used for static classification, it is the input of discrete image samples.\" Therefore, images can be sequentially input into the ESN in a certain order (from t1 to t60000), but there is no sequential dependency between the images themselves.\nThe output is described as follows when training data: \"the transfer and update process is as follows. REb stand for the update status.\" Here, REb represents the reservoir state space. In simple terms, one image input corresponds to one reservoir output.\nTo further facilitate the understanding, we have added additional content to Fig. S19. Specifically, we have provided a more detailed explanation of the image flattening process, as illustrated in Fig. S19(b), where the images are directly flattened and then input into the reservoir. Fig. S19(c) presents the specific architecture of the system. After flattening, the images are sequentially fed into the physical ESN according to the original order of the dataset. At this stage, the ambipolar transistor array generates the corresponding reservoir state current data. These data, together with the original input data and a bias term (1), are concatenated for training the readout layer. Consequently, when the reservoir size is set to 1000, the resulting shape of the readout layer is 1785 × 10.\nWe have supplemented the following content in the SI:\nFig. S19 (a) Schematic diagram illustrating the time series representation of the MNIST handwritten digit dataset. (b) Illustration of the image flattening process applied to the MNIST handwritten digits. (c) Schematic of the physical ESN architecture used for recognizing the MNIST dataset.\nIn addition, we have incorporated further explanations and clarifications in the revised manuscript to address your concerns. We hope these additions will help clarify our methodology and enhance the reader’s understanding of our approach.\nThe process of flattening the images, as shown in Fig. S19 (b). These images are directly flattened and placed into the reservoir. Fig. S19 (c) shows the specific architecture. After the images are flattened, they are sequentially input into the physical ESN according to the original order of the dataset. At this point, the ambipolar transistor array will output the corresponding reservoir state current data. These data, along with the original data and bias data (1), are stacked together for training the readout layer. Therefore, when the size of the reservoir is 1000, the shape of the readout layer is 1785 × 10.\n\nReviewer #2 (Remarks to the Author):\nThe authors have addressed all my comments, and I recommend the manuscript to be accepted for publication.\nResponse:\nWe sincerely appreciate your time and thoughtful review of our manuscript. Your constructive comments and suggestions have significantly improved the quality of our work.\n\nReviewer #3 (Remarks to the Author):\nThe authors have addressed most of my concerns. However, several aspects of the experimental methodology and descriptions remain unclear. Please see the detailed comments below:\nQ1. The label in Fig. S5 appears to be incorrect.\nResponse:\nThank you very much for your careful review and for pointing out the issue with the label in Fig. S5. We apologize for this oversight. We have corrected the label in the revised version of the supplementary materials (Fig. S5a) to ensure accuracy and clarity. Thank you again for bringing this to our attention.\nWe have modified the following contents in SI:\nFig. S5 (a) The switch speed of ambipolar transistor with pulse width of 2 μs.\nQ2. The authors have successfully fabricated a thinner Al O layer that enables lower\n2 3\noperating voltage. Given this advancement, it is unclear why this improved device was not used in the main experiments. Are there any potential limitations? Clarifying this point would be highly valuable for guiding future research in the field.\nResponse:\nThank you for your insightful comment. We appreciate your recognition of our successful fabrication of a thinner Al O layer, which indeed enables lower operating\n2 3\nvoltage. The reason we did not use this improved device in the main experiments is as follows:\nWhile the thinner Al O layer offers the advantage of reduced operating voltage,\n2 3\nit also introduces certain limitations. Specifically, we observed that devices with thinner Al O layers exhibited increased leakage current and reduced long-term\n2 3\nstability during repeated operation cycles. These issues affect the reliability of the experimental results. When we currently use thinner devices to make arrays, due to insufficient repeatability, a large amount of noise occurs, resulting in poor recognition accuracy. As shown in the Fig. S34, we use a 10nm Al O dielectric layer array to\n2 3\nrecognize a dataset of English letters, with an accuracy rate of 79.50%.\nWe have supplemented the following content in the SI:\nFig. S34. Accuracy of physical ESN achieved by ambipolar transistor array with 10 nm oxide layer on the English letter data set.\nWhen we currently use thinner devices to make arrays, due to insufficient repeatability, a large amount of noise occurs, resulting in poor recognition accuracy. As shown in Fig. S34, we use a 10 nm Al O dielectric layer array to recognize a\n2 3\ndataset of English letters, with an accuracy rate of 79.50%.\nWe have supplemented the following content in the manuscript:\nDue to the insufficient stability of the 10 nm Al O dielectric layer, there are\n2 3\nchallenges regarding the long-term reliability of the ambipolar transistor array. As a result, the performance of the physical ESN constructed with thin Al O dielectric\n2 3\nlayer needs to be improved in the future (Supplementary Fig. S34).\nTherefore, for the main experiments, we chose to use devices with a thicker\nAl O layer to ensure optimal balance between performance and stability. We agree\n2 3\nthat further optimization of the Al O layer is an important direction for future\n2 3\nresearch, and we are currently working on strategies to mitigate the limitations associated with thinner layers.\nQ3. While the Fig. S5 states that the device operates at nanosecond-level switching speeds, the experiments were performed using microsecond-long pulses. This discrepancy should be addressed.\nResponse:\nThank you very much for your careful observation and for highlighting this important point. We have conducted additional experiments using nanosecond-scale pulses (as short as 200 ns) and confirmed that the device maintains its switching capability at these shorter timescales as shown in Fig. S5 (b).\nFig. S5 (b) Switch speed of the device with 200 ns pulse width.\nQ4. The calculation method for power consumption remains unclear. The reported values cannot be reproduced from the data provided. The authors should include a step-by-step description or an example calculation.\nResponse:\nThank you very much for your valuable feedback regarding the calculation of power consumption. We apologize for the lack of clarity in our previous description and for any difficulties in reproducing the reported values.\nOur calculation method is completed as follows. We took five data points when we applied a 100-microsecond pulse, and then we added them up. For positive pulses:\nThe collected I are 4.742, 7.340, 9.010, 9.680 and 9.841 nA. For negative pulses:\nds\nThe collected I are-2.735, -6.037, -8.216, -9.645 and -10.12 nA. According to the\nds\nformula:\nI =∑5 I (1)\nds i=1 ds i\nP=I ×t ×V (2)\nds pulse ds\nThe absolute values of the pulse width (t ) are 100 μs and the V is 0.02 V.\npulse ds\nThe sum of positive and negative I pulses is 40.615 nA and -36.75 nA respectively.\nds\nTherefore, the power consumption is approximately 812 fJ and 735 pJ.\nHowever, we once again refer to the literature (“Polyvinyl alcohol electrolyte-gated oxide transistors with tetanization activities for neuromorphic computing”, page 6, “The energy consumption (P) could be estimate dusing the following equation: P = V × I × t” and “An extremely low-power-consumption\nds peak\nreconfigurable two-dimensional tellurene artificial synapse for bio-inspired wearable edge computing”, page 5, “Furthermore, to assess its potential as an artificial synaptic device with edge-computing capabilities, we calculated the power consumption using the following equation: E = V ×I ×t ”), which universally considers\nconsumption read peak pulse\nthe peak current value as the entire pulse. So, we selected the I with the maximum\nds\npulse duration for calculation and adjusted the results accordingly. For positive pulses, the power consumption is 9.841 nA × 0.02 V × 100 μs ≈ 197 fJ. For negative pulses, the power consumption is 10.12 nA × 0.02 V × 100 μs ≈ 202 fJ.\nTherefore, we have corrected this result in Table S1 and Fig. S10 (c) and (f).\nFig. S10 (c) Power consumption under positive V , (f) Power consumption under\nds\nnegative V .\nds\n[R5] Huang, Y. J. et al. Polyvinyl alcohol electrolyte-gated oxide transistors with tetanization activities for neuromorphic computing. Journal of Materials Chemistry C 12, 5166-5174 (2024).\n[R6] You, B. et al. An extremely low-power-consumption reconfigurable two-dimensional tellurene artificial synapse for bio-inspired wearable edge computing. Journal of Materials Chemistry C 12, 6596-6605 (2024).\nTable S1. Power consumption of neuromorphic devices.\nMaterials Type Power DATASE accurac Ref.\nT y\nMoS h-BN transistor ≈450 pJ MNIST 91.1 % 6\n2/\nMoS /Nb O transistor ≈6 pJ MNIST 94.2% 7\n2 2 5\nMoSe transistor ≈140 nJ MNIST ~70 % 8\n2\nα-In Se transistor ≈4 nJ MNIST 91.9 % 9\n2 3\n≈77.24p\n180nm CMOS transistor 10\nJ\nIGZO transistor 0.15 nJ MNIST 73% 11\nCsPbBr QDs transistor 1nJ MNIST 73% 12\nIGZO/HfZrOx transistor 5.69nJ MNIST 91.1% 14\nAlGaN/GaN transistor 18.5 pJ 15\nSi doped-GaN transistor 33.4 pJ 16\nGa O /ZnO heterojunction 236 pJ MNIST 92.31% 17\n2 3\nPTCDA/MoS transistor 10 pJ 18\n2\nphototransisto\nSi/P3HT/PbS 15.8 pJ 19\nrs\nphototransistor\nGraphene/SWNTs 2.5nJ 20\npentacene/PMMA/CsPbBr\n3\nflash memory 1.4 nJ 21\nQD\nSi NCs phototransistor 9.14 nJ 94% 22\nITZO phototransistor 2.5 nJ 90.3% 23\nPorphyrin/Sinanowires phototransistor 1 nJ 24\nP3HT/CsPbBr QD phototransistor 250 nJ 25\nx\nTFT-CN/DTT-8 phototransistor 5.3nJ 27\nFe O /WS memristor 72 fJ 94.3% 28\n2 3 2\nTiN transistor 30fJ 29\nThis\nP3HT/ZnO transistor 202 fJ MNIST 95.27% wor\nk\nQ5. The working flow of the network described in Fig. S13 is still not clearly explained. The input image is flattened and directly passed to an array, raising the question of why a reservoir is necessary. A simple linear layer could classify the inputs in a single step. The authors should further justify the use of the reservoir architecture.\nResponse:\nWe appreciate the reviewer’s insightful question regarding the necessity of the reservoir in our architecture. We agree that, for simple classification tasks, a linear layer alone could achieve satisfactory results. However, our intention in employing such straightforward tasks is to provide a clear demonstration of the operational mechanisms and underlying principles of the reservoir architecture.\nAs shown in Figure 2c in manuscript, the spatial states generated by the reservoir output illustrate the additional features produced by this framework. These features are nonlinear and represent dynamic, high-dimensional data that a simple linear layer cannot capture. While this simple dataset may not fully showcase the advantages of the reservoir architecture, it serves as an effective example to help readers understand how the reservoir processes and transforms input data. Furthermore, we have conducted subsequent experiments on the MNIST dataset, which further validate the capabilities and effectiveness of our architecture in more complex scenarios. To improve clarity, we have also added supplementary explanations to Fig. S13. Thank you again for your valuable feedback.\nFig. S13 Schematic diagram of converting English letters into time series.\nThe 5 × 5 letter image is flattened into a 25-dimensional vector in sample order as a temporal input, which is then passed through a fixed linear layer of 25 × 16 to match the dimensions to a 16 × 16 physical reservoir array. The nonlinearity of the transistor array is used to obtain a 16-dimensional reservoir state. The original 25-dimensional data, reservoir state, and bias are concatenated into a 42-dimensional feature, which is then mapped to 6 classes for recognition using linear regression.\nQ6. It is unclear why the authors treat different samples in MNIST as different time steps. This approach is confusing. Do the authors intend for the reservoir to extract temporal features across samples? Further clarification and justification are needed.\nResponse:\nWe sincerely appreciate the reviewer’s thoughtful comment regarding our approach of treating different MNIST samples as distinct time steps. Our motivation for transforming images into sequences was to assess the reservoir’s capability for sequential modeling. As highlighted in the existing literature, the use of ESNs for processing MNIST data is well established and provides valuable insights. ESNs are often incorporated into various neural network architectures to improve image recognition accuracy. Consequently, it is a common practice to convert MNIST data into a time series format when applying ESNs, and the resulting accuracy is regarded as a measure of the physical ESN’s performance.\nIn the work by Sun et al., a similar methodology was adopted. They explicitly stated, “Different from processing sequence problems, when ESN is used for static classification, it is the input of discrete image samples.” This approach is frequently employed for image recognition tasks using ESNs. In conventional ESN implementations, it is standard practice to convert the rows or columns of an MNIST image into a time series, thereby enabling the network to capture temporal features within each sample.\nFig. 2. Inputs and outputs of a ESN for image classification, with N = 28 (Schaetti. et\ni\nal).\nIn this approach, the time series characteristics are derived from the internal structure of each image, and the dataset inherently possesses temporal attributes. We applied the same processing method to our physical ESN, and the results are presented below: the model achieved accuracies of 91.92% and 84.59% on the digit and fashion datasets, respectively. While these results surpass those obtained with a simple linear layer, they are still lower than the 94.38% and 85.33% accuracies reported in the main article.\nWe believe this decrease in accuracy may be attributed to the introduction of artificial temporal relationships between samples, which do not naturally exist and may not be suitable for this type of data. Based on these findings, we opted for a cross-sample processing approach in our study. We hope this clarification addresses your concerns, and we have revised the manuscript to provide a more detailed explanation and justification for our methodology. Thank you again for your valuable feedback.\nWe have supplemented the following content in the SI:\nFig. S22 (a) The accuracy of identifying the digits and fashion datasets using a column-by-column installation method for the physical ESN.\nAs show in Fig. S22 (a), the 28 pixels of each column vector in the MNIST dataset as a time-series signal and achieved recognition accuracies of 91.92 % on the digits dataset and 84.59 % on the fashion dataset.\nWe have also supplemented the following content in the manuscript:\nIn addition, we encoded the 28 pixels of each column vector in the MNIST dataset as a time-series signal and achieved recognition accuracies of 91.92 % on the digits dataset and 84.59 % on the fashion dataset (Supplementary Fig. S22(a))\n[R7] Schaetti, N., Salomon, M. & Couturier, R. in 2016 19th IEEE Intl Conference on Computational Science and Engineering (CSE), IEEE 14th Intl Conference on Embedded and Ubiquitous Computing (EUC), and 15th Intl Symposium on Distributed Computing and Applications for Business Engineering (DCABES).\n484-491.\nQ7. Fig. S22 does not clearly demonstrate the effectiveness of the reservoir. It is also unclear whether the dimension of the linear layer was modified. If not, how were the dimensional mismatches handled? The authors should provide a more detailed explanation of the setup and clarify how the comparison was made.\nResponse:\nThank you for your thoughtful comments and for highlighting these important points. We apologize for any confusion caused by our previous description.\nTo address your concerns, we confirm that the dimensions of the linear layer were indeed modified to accommodate the unique architecture of the ESN.\nSpecifically, the linear layer’s dimensions were expanded to include not only the original feature dimensions of each dataset, but also the reservoir state dimensions and a bias term. The rationale behind this design is to leverage the reservoir’s ability to introduce high-dimensional, nonlinear features, thereby enhancing classification accuracy.\nTo further clarify the experimental setup, we have added schematic illustrations for both English-letter recognition and the MNIST dataset (see Fig. S13 and Fig. S19 in the revised manuscript). For additional verification, we performed linear classification using both the original feature dimension (784) and the expanded dimension (1785, corresponding to the reservoir output) on two datasets, without involving the reservoir. When processing MNIST data with the ESN, the readout layer dimension is 1785. Therefore, in Fig. S22 (b), we present classification results using linear layers of both 784 and 1785 dimensions. The results indicate that simply increasing the dimensionality of the linear layer, without the reservoir, does not lead to improved performance.\nWe have supplemented the following content in the SI:\nFig. S22 (b) The classification results for the digits and fashion datasets with 784- and 1785-neuron linear layers.\nFig. S22 (b) shows the classification results of 784 and 1785 linear layers on the digits and fashion datasets. For the digits dataset, their accuracy rates were both 85.34%. For the fashion dataset, their accuracy rates were 81.12%.\nWe hope this explanation clarifies the setup and the basis for our comparisons. Thank you again for your valuable feedback."}, {"round": 3, "reviewer_comments": {"reviewer_1": "I do appreciate the revision work from the author. Here are my comments.\n1. The method of assigning different input images to different time steps is still confusing. Let's say, 'U' is T1, and 'X' is t2. Then, driving the reservoir with 'U' will leave some memory to it, and this will affect classifying the input 'X'. This does not make sense, as classifying 'U' should not affect classifying 'X'.\n2. Otherwise, probably the reservoir does not have the memory capacity, and only performs non-linear transformation of the input data. (In such case, the word 'step' will be confusing and should be replaced with 'sample #', etc.) However, in the supplementary information, there are discussions about the temporal response of the device (e.g. Fig. S11). So this might not be the case.\n3. In summary, my only question is how different numbers/characters can form a 'time series', and led to the successful classification. Usually, in other literatures, people disassembles a single image into multiple time-series segments (such as binary 4 digit series), to perform the physical reservoir computing.", "reviewer_3": "While I appreciate the authors' efforts in their response, the authors have not convincingly demonstrated the necessity or effectiveness of the proposed reservoir component. And I believe some of the results remain questionable or potentially incorrect.\nFundamentally, the MNIST task is a static image recognition problem as the authors mentioned in their reply. And in their setup, each time step corresponds to a single input image, and no temporal sequence is involved. This effectively eliminates any temporal dynamics, which are the core motivation for using a reservoir. Without sequential input, the reservoir cannot perform dynamic processing, and its role is reduced to that of a static transformation layer. Therefore, applying a reservoir in this context appears conceptually unjustified.\nIn addition, the system architecture is different from the conventional reservoir design. Between the input and output, the authors introduce an additional fully connected layer, creating a multi-layer feedforward structure. This design significantly complicates the interpretation of results, as it becomes unclear whether the reported performance is attributable to the reservoir or to the fully connected layer. Given that a standard fully connected neural network is already known to perform very well on MNIST, it is likely that the performance gain originates primarily from this component rather than from the reservoir itself.\nTo support the role of the reservoir, the authors present results comparing single-layer and two-layer fully connected networks, claiming that their system achieves better accuracy. However, these results are questionable. For one, the reported accuracy of the baseline networks is surprisingly low, far below what is commonly achieved on MNIST. A single-layer fully connected network typically exceeds 90% accuracy, and a two-layer network with sufficient neurons generally surpasses 95%. Furthermore, the authors report identical performance and confusion matrices for both the two-layer and three-layer networks, which is highly implausible given normal variance in training. These inconsistencies suggest that the experiments do not reliably isolate or validate the reservoir’s contribution and may even indicate issues in the experimental procedure itself."}, "author_response": "Reviewer #1\n1. The method of assigning different input images to different time steps is still confusing. Let's say, 'U' is T1, and 'X' is t2. Then, driving the reservoir with 'U' will leave some memory to it, and this will affect classifying the input 'X'. This does not make sense, as classifying 'U' should not affect classifying 'X'.\nResponse:\nThank you for raising this important point. In our implementation, the reservoir state is always reset to its initial resting value (or to a randomly selected state from the same distribution) before each new input sequence is presented. This ensures that the reservoir does not retain any memory from previous inputs. As a result, when an image such as 'U' is presented at time step t1, the reservoir processes it starting from a clean, unbiased state. Similarly, when a different image such as 'X' is presented at t2, the reservoir is again initialized to the same baseline condition. This procedure prevents any cross-contamination between trials, ensuring that the classification of one sample does not influence the outcome of another. The specific reset protocol is illustrated by the interval of conditional pulses, as shown in the figure below.\nFig. S35. Device I-t performance at different pulse intervals.\n2. Otherwise, probably the reservoir does not have the memory capacity, and only performs non-linear transformation of the input data. (In such case, the word 'step' will be confusing and should be replaced with 'sample #', etc.) However, in the supplementary information, there are discussions about the temporal response of the device (e.g. Fig. S11). So this might not be the case.\nResponse:\nThank you for your insightful comment. In response to Reviewer #3’s request, we conducted additional experiments, the results of which are presented in the Supplementary Fig. S11. These data confirm that the reservoir does possess short-term memory capacity, which can be utilized for tasks such as English letter recognition. However, we found that the memory effect is not sufficiently robust for optimal performance in our primary application. As a result, we chose not to rely on this dynamic encoding approach in the main manuscript, and instead adopted a static encoding scheme, which demonstrated superior accuracy in our comparative studies.\nFig. S36. 5 bit short-term memory effect.\nFig. S37. 5-bit dynamic response for MNIST recognition.\n3. In summary, my only question is how different numbers/characters can form a 'time series', and led to the successful classification. Usually, in other literatures, people disassembles a single image into multiple time-series segments (such as binary 4 digit series), to perform the physical reservoir computing.\nResponse:\nThank you very much for your thoughtful question. The 4-bit dynamic response is a widely used encoding strategy for recognizing the MNIST dataset in physical reservoir computing. To implement this approach, the dataset is first binarized, typically by setting all pixel values greater than zero to 1, and the rest to 0. As illustrated in Fig. S38, this binarization allows each image to be decomposed into a sequence of 4-bit segments, effectively forming a time series. This encoding scheme generates 16 possible states, which capture the temporal correlations among different subregions of the image. These temporal correlations enable the reservoir to perform nonlinear transformations, thereby allowing a simple linear readout layer to achieve effective classification. Please see the detailed explanation for our design.\nFig. S39 presents the recognition results of a linear readout layer applied to both binary-encoded and 4-bit encoded images. In the 4-bit encoding scheme, each image is represented by 16 distinct conductivity states, corresponding to the device’s different responses. The results show that binary encoding slightly outperforms 4-bit encoding, as the latter compresses image features and thus reduces their recognizability.\nIn Fig. S40, we captured all four characteristic bits from the device response. By sampling the final value within each pulse duration as a feature, we obtained the mapping diagram shown in Fig. S3a, achieving an accuracy of 87.94%. These results demonstrate that implementing dynamic mechanisms can be an effective strategy for enhancing image recognition accuracy.\nTo further improve accuracy, we introduced a sparse matrix. As shown in Fig. S41, we constructed a (28×28) matrix in which the entries with a value of 1 represent active devices. These devices perform feature mapping through dynamic mechanisms, and their temporal node mappings across columns are aggregated to form higher-level latent space states. By combining these spatial states with the raw data for recognition by the readout layer, we achieved an accuracy of 89.40%. To further enhance performance, we adopted a -1 and +1 strategy, where -1 represents a negative response mechanism (i.e., the positive response multiplied by -1), and +1 represents a positive response. As shown in Fig. S42, this approach increased the accuracy to 89.42%.\nTo further boost accuracy, we converted these mappings into voltage signals and leveraged the device’s bipolarity to achieve an additional nonlinear transformation. This nonlinear mapping, which mimics the tanh function through bipolar modulation, resulted in an accuracy improvement to 91.54%, as illustrated in Fig. S43. The device’s bipolarity thus serves as an effective means to enhance data nonlinearity. Seeking further improvement, we expanded the dimension of the sparse matrix from (28, 28) to (28, 280), which increased the accuracy to 92.29%, as shown in Fig. S44. However, when the sparse matrix was further expanded to (28, 560), the accuracy decreased to 91.29%, as shown in Fig. S45.\nGiven this result, we abandoned the dynamic mechanism and instead flattened the image to construct a sparse matrix of (784, 280) dimensions to assist the bipolar neurons in feature extraction. This approach yielded an accuracy of 89.83%, which was lower than the dynamic mechanism-based result of 92.29%, as shown in Fig. S46.\nHowever, when we further increased the size of the sparse matrix to (784, 800), the accuracy reached 92.30%, nearly matching the result achieved with the dynamic mechanism, as shown in Fig. S47. When the full grayscale values of the MNIST dataset were used, the accuracy improved to 92.81%, as shown in Fig. S48. Finally, by expanding the sparse matrix to 1000 columns and controlling the sparsity at 5%, we achieved our highest accuracy of 94.14%, as shown in Fig. S49.\nFig. S38. 4 bit dynamic mechanism.\nFig. S39. (a) binarized image, (b) confusion matrix of recognition binary numbers, (c) 4bit encoded features (16 states) and (d) confusion matrix of recognition 4-bit encoded features.\nFig. S40. (a) 4-bit encoding (64 states), (b) Confusion matrix.\nFig. S41. (a) Sparse matrix (0 and 1) for higher level features, (b) higher level reservoir state features, (c) binarized images and (d) confusion matrix.\nFig. S42. (a) Sparse matrix (-1, 0 and 1) for higher level features, (b) higher level reservoir state features, (c) binarized images and (d) confusion matrix.\nFig. S43. Identification tasks based on ambipolar neurons and dynamic mechanisms (a) Sparse matrix (-1, 0 and 1) for higher level features, (b) higher level reservoir state features, (c) binarized images and (d) confusion matrix.\nFig. S44. Identification tasks based on ambipolar neurons and dynamic mechanisms (a) binarized images (b) Sparse matrix (-1, 0 and 1) for higher level features (28, 280), (c) higher level reservoir state features and (d) confusion matrix.\nFig. S45. Identification tasks based on ambipolar neurons and dynamic mechanisms (a) binarized images, (b) Sparse matrix (-1, 0 and 1) for higher level features (28,560), (c) higher level reservoir state features and (d) confusion matrix.\nFig. S46. Identification tasks based on ambipolar neurons and without dynamic mechanisms (a) binarized images, (b) Sparse matrix (-1, 0 and 1) for higher level features (784, 280) and confusion matrix.\nFig. S47. Identification tasks based on ambipolar neurons and without dynamic mechanisms (a) binarized images, (b) Sparse matrix (-1, 0 and 1) for higher level features (784, 800) and confusion matrix.\nFig. S48. Identification tasks based on ambipolar neurons and without dynamic mechanisms (a) full gray value images, (b) Sparse matrix (-1, 0 and 1) for higher level features (784, 800) and confusion matrix.\nFig. S49. Identification tasks based on ambipolar neurons and without dynamic mechanisms (a) full gray value images, (b) Sparse matrix (-1, 0 and 1) for higher level features (784, 1000) and confusion matrix.\nIn Figure S35, we demonstrate how dynamic mechanism conversion to static mechanisms is achieved through controlled pulse interval control. Fig. S36 and S37 introduce the 5-bit dynamic mechanism and English letter recognition tasks (Supplementary Fig. S34). Due to the limited features of the dynamic mechanism, the English letters cannot be fully recognized by the layer. Figures S38 to S49 showcase the performance of dynamic mechanisms versus static MNIST digit recognition tasks. Traditional 4-bit digital feature extraction yields lower results than binarized methods due to compressed digital features. By fully utilizing 4-bit response data, accuracy can be improved from 84.92% to 87.94%. This demonstrates that the nonlinearity of dynamic mechanisms enhances binarized image recognition accuracy. Furthermore, we designed sparse matrices to obtain higher-order nonlinear features. The dynamic mechanism achieves 92.29% accuracy, while the static model reaches 94.14% accuracy for full grayscale images. Accuracy may further increase with higher-dimensional sparse matrices (Supplementary Fig. S38 to S49).\nDue to the nonlinear behavior of the device, the linear matrix can be used to enhance the nonlinearity. The matrix multiplication with device performance cannot improve the accuracy and nonlinearity.\nTherefore, we initially developed this work using a time-scale approach and, after multiple rounds of optimization, ultimately retained this temporal framework. During data processing, we systematically extracted features from all 60,000 images in chronological order, with the architecture employing temporal sequences as its coordinate system.\n\nReviewer #3\nWhile I appreciate the authors' efforts in their response, the authors have not convincingly demonstrated the necessity or effectiveness of the proposed reservoir component. And I believe some of the results remain questionable or potentially incorrect.\nFundamentally, the MNIST task is a static image recognition problem as the authors mentioned in their reply. And in their setup, each time step corresponds to a single input image, and no temporal sequence is involved. This effectively eliminates any temporal dynamics, which are the core motivation for using a reservoir. Without sequential input, the reservoir cannot perform dynamic processing, and its role is reduced to that of a static transformation layer. Therefore, applying a reservoir in this context appears conceptually unjustified.\nResponse:\nThank you for your thoughtful comments and for highlighting this important conceptual issue. We have conducted comprehensive computational analyses to directly compare the performance of dynamic and static mechanisms within our system. As detailed in Figures S38 to S49, our results consistently show that static encoding outperforms dynamic encoding for the MNIST image recognition task. This finding is transparently presented and discussed in the manuscript.\nWhile it is true that MNIST is a static image dataset, our study systematically explored both temporal (dynamic) and static encoding schemes to rigorously assess the necessity and effectiveness of the reservoir component. The superior performance of the static approach in our experiments suggests that, for this particular application, the nonlinear transformation capability of the reservoir, rather than its temporal dynamics, is the key factor driving improved accuracy. This insight is significant, as it demonstrates that physical reservoir computing architectures can be flexibly adapted to different types of input data, and that their effectiveness is not solely dependent on temporal processing.\nWe believe these results provide a clear and honest evaluation of the reservoir’s role in our system, and offer valuable guidance for future research on physical reservoir computing in both static and dynamic contexts.\nAs illustrated in Fig. S38, this binarization allows each image to be decomposed into a sequence of 4-bit segments, effectively forming a time series. This encoding scheme generates 16 possible states, which capture the temporal correlations among different subregions of the image. These temporal correlations enable the reservoir to perform nonlinear transformations, thereby allowing a simple linear readout layer to achieve effective classification. Please see the detailed explanation for our design.\nFig. S39 presents the recognition results of a linear readout layer applied to both binary-encoded and 4-bit encoded images. In the 4-bit encoding scheme, each image is represented by 16 distinct conductivity states, corresponding to the device’s different responses. The results show that binary encoding slightly outperforms 4-bit encoding, as the latter compresses image features and thus reduces their recognizability.\nIn Fig. S40, we captured all four characteristic bits from the device response. By sampling the final value within each pulse duration as a feature, we obtained the mapping diagram shown in Fig. S3a, achieving an accuracy of 87.94%. These results demonstrate that implementing dynamic mechanisms can be an effective strategy for enhancing image recognition accuracy.\nTo further improve accuracy, we introduced a sparse matrix. As shown in Fig. S41, we constructed a (28×28) matrix in which the entries with a value of 1 represent active devices. These devices perform feature mapping through dynamic mechanisms, and their temporal node mappings across columns are aggregated to form higher-level latent space states. By combining these spatial states with the raw data for recognition by the readout layer, we achieved an accuracy of 89.40%. To further enhance performance, we adopted a -1 and +1 strategy, where -1 represents a negative response mechanism (i.e., the positive response multiplied by -1), and +1 represents a positive response. As shown in Fig. S42, this approach increased the accuracy to 89.42%.\nTo further boost accuracy, we converted these mappings into voltage signals and leveraged the device’s bipolarity to achieve an additional nonlinear transformation. This nonlinear mapping, which mimics the tanh function through bipolar modulation, resulted in an accuracy improvement to 91.54%, as illustrated in Fig. S43. The device’s bipolarity thus serves as an effective means to enhance data nonlinearity. Seeking further improvement, we expanded the dimension of the sparse matrix from (28, 28) to (28, 280), which increased the accuracy to 92.29%, as shown in Fig. S44. However, when the sparse matrix was further expanded to (28, 560), the accuracy decreased to 91.29%, as shown in Fig. S45.\nGiven this result, we abandoned the dynamic mechanism and instead flattened the image to construct a sparse matrix of (784, 280) dimensions to assist the bipolar neurons in feature extraction. This approach yielded an accuracy of 89.83%, which was lower than the dynamic mechanism-based result of 92.29%, as shown in Fig. S46.\nHowever, when we further increased the size of the sparse matrix to (784, 800), the accuracy reached 92.30%, nearly matching the result achieved with the dynamic mechanism, as shown in Fig. S47. When the full grayscale values of the MNIST dataset were used, the accuracy improved to 92.81%, as shown in Fig. S48. Finally, by expanding the sparse matrix to 1000 columns and controlling the sparsity at 5%, we achieved our highest accuracy of 94.14%, as shown in Fig. S49.\nFig. S38. 4 bit dynamic mechanism.\nFig. S39. (a) binarized image, (b) confusion matrix of recognition binary numbers, (c) 4bit encoded features (16 states) and (d) confusion matrix of recognition 4-bit encoded features.\nFig. S40. (a) 4-bit encoding (64 states), (b) Confusion matrix.\nFig. S41. (a) Sparse matrix (0 and 1) for higher level features, (b) higher level reservoir state features, (c) binarized images and (d) confusion matrix.\nFig. S42. (a) Sparse matrix (-1, 0 and 1) for higher level features, (b) higher level reservoir state features, (c) binarized images and (d) confusion matrix.\nFig. S43. Identification tasks based on ambipolar neurons and dynamic mechanisms (a) Sparse matrix (-1, 0 and 1) for higher level features, (b) higher level reservoir state features, (c) binarized images and (d) confusion matrix.\nFig. S44. Identification tasks based on ambipolar neurons and dynamic mechanisms (a) binarized images (b) Sparse matrix (-1, 0 and 1) for higher level features (28, 280), (c) higher level reservoir state features and (d) confusion matrix.\nFig. S45. Identification tasks based on ambipolar neurons and dynamic mechanisms (a) binarized images, (b) Sparse matrix (-1, 0 and 1) for higher level features (28,560), (c) higher level reservoir state features and (d) confusion matrix.\nFig. S46. Identification tasks based on ambipolar neurons and without dynamic mechanisms (a) binarized images, (b) Sparse matrix (-1, 0 and 1) for higher level features (784, 280) and confusion matrix.\nFig. S47. Identification tasks based on ambipolar neurons and without dynamic mechanisms (a) binarized images, (b) Sparse matrix (-1, 0 and 1) for higher level features (784, 800) and confusion matrix.\nFig. S48. Identification tasks based on ambipolar neurons and without dynamic mechanisms (a) full gray value images, (b) Sparse matrix (-1, 0 and 1) for higher level features (784, 800) and confusion matrix.\nFig. S49. Identification tasks based on ambipolar neurons and without dynamic mechanisms (a) full gray value images, (b) Sparse matrix (-1, 0 and 1) for higher level features (784, 1000) and confusion matrix.\nIn addition, the system architecture is different from the conventional reservoir design. Between the input and output, the authors introduce an additional fully connected layer, creating a multi-layer feedforward structure. This design significantly complicates the interpretation of results, as it becomes unclear whether the reported performance is attributable to the reservoir or to the fully connected layer. Given that a standard fully connected neural network is already known to perform very well on MNIST, it is likely that the performance gain originates primarily from this component rather than from the reservoir itself.\nTo support the role of the reservoir, the authors present results comparing single-layer and two-layer fully connected networks, claiming that their system achieves better accuracy. However, these results are questionable. For one, the reported accuracy of the baseline networks is surprisingly low, far below what is commonly achieved on MNIST. A single-layer fully connected network typically exceeds 90% accuracy, and a two-layer network with sufficient neurons generally surpasses 95%. Furthermore, the authors report identical performance and confusion matrices for both the two-layer and three-layer networks, which is highly implausible given normal variance in training. These inconsistencies suggest that the experiments do not reliably isolate or validate the reservoir’s contribution and may even indicate issues in the experimental procedure itself.\nResponse:\nThank you for raising these concerns regarding our system architecture and the interpretation of our results. We would like to clarify that the additional layer introduced between the input and output is not a conventional fully connected neural network layer, but rather a linear transformation trained using the least squares method. This approach is fundamentally different from standard feedforward neural networks, which typically employ non-linear activation functions and iterative training algorithms.\nOur intention was to isolate and highlight the contribution of the physical reservoir component, rather than to rely on algorithmic enhancements. The linear layers in our system serve only to map the features extracted by the reservoir, without introducing additional nonlinearity. As a result, the observed performance improvements can be attributed to the intrinsic properties of the reservoir devices.\nRegarding the baseline results, we acknowledge that the reported accuracies for single-layer and two-layer linear networks are lower than those typically achieved by fully connected neural networks on MNIST. This is because our comparison is based on strictly linear models, without the benefit of non-linear activations or deep learning techniques. To ensure transparency and reproducibility, we have made our code and data publicly available for independent verification."}, {"round": 4, "reviewer_comments": {"reviewer_1": "The response letter indicates that the the ambipolar transistors are performing static computation. While authors additionally tried utilizing dynamic aspect of the devices in the revision, it seems to be not a main objective of the research. Given that the concept of 'reservoir computing' is proposed for processing dynamical data with lower computational cost, the suggested research is difficult to be categorized as physical reservoir computing, from my understanding. Therefore, I regret to say that the paper is not suitable for publication in Nature Communications.", "reviewer_3": "The authors have addressed my concerns in the latest revision and the manuscript may be accepted."}, "author_response": "Response:\nWe sincerely thank the reviewers for your valuable comments. The issues you raised are very important and have indeed provided directions for further improvement of this study. As you rightly pointed out, our current work primarily focuses on leveraging device characteristics for static computing functions, and exploration in processing dynamic data remains insufficient, which creates a certain gap with the concept of typical physical reservoir computing. In response to this concern, we have made significant additions to the research content during the revision.\nWe propose a novel paradigm that integrates dynamic and static characteristics. This method introduces a dynamic response mechanism to incorporate temporal attributes into data preprocessing, followed by solving through a physically-inspired Echo State Network (ESN) paradigm. This design effectively overcomes the limitations of traditional static approaches, achieving higher-level performance. On the MNIST handwritten digit dataset and the Fashion-MNIST dataset, the proposed paradigm achieves classification accuracies of 96.98% and 86.67%, respectively. In the task of Lorenz attractor prediction, the Normalized Root Mean Square Error (NRMSE) reaches 0.165%, 0.280%, and 0.206%, demonstrating performance comparable to state-of-the-art algorithms. The computational process and relevant explanations of this paradigm are detailed below.\nFig. R1 4-bit encoding with varying pulse intervals. a 5 ms, b 10 ms, c 15 ms, d 20 ms, e 25 ms, f 30 ms.\nWe applied a 5 ms pulse voltage to the gate and subsequently measured the I . ds After the pulse was applied, different time intervals were set to achieve varying levels of 4bit computation. A 4bit encoding scheme was adopted, with pulse intervals set at 5 ms, 10 ms, 15 ms, 20 ms, 25 ms, and 30 ms, respectively. As shown in Fig. R1a, when the pulse interval is 5 ms, the device exhibits clear 4bit encoding characteristics, and its 16 distinguishable conductance states can be used for image encoding. However, as shown in Fig.R1b, when the pulse interval is increased to 10 ms, all conductance states except for the unstimulated code \"0000\" shift toward lower conductance values. This shift introduces distinct feature representations for image processing. When the pulse interval is further increased to 15 ms (Fig. R1c), the conductance states continue to shift toward lower values, indicating that changes in pulse interval alter the device’s response pattern.\nTaking the code \"1110\" as an example for analysis: when the pulse interval is 5 ms, the drain current at the first recorded point (Step 1) is approximately 39.95 nA, while at the last recorded point (Step 4) it is approximately 53.00 nA. When the interval increases to 10 ms, the Step 1 current decreases to 38.70 nA, and Step 4 decreases to 48.64 nA. At 15 ms, Step 1 is 37.45 nA and Step 4 is 39.30 nA. As observed in these three cases, the Step 1 current consistently remains higher than that of Step 4, but the gap between them gradually narrows. This indicates that the first three pulses significantly influence the output of the last pulse, and this influence weakens as the interval increases. This trend reverses when the pulse interval reaches 20 ms (Fig. R1d). Here, the Step 1 current is 35.93 nA, while Step 4 is 35.37 nA, meaning Step 1 begins to fall below Step 4. As the interval further increases to 25 ms (Fig. R1e) and 30 ms (Fig. R1f), the Step 1 currents are 35.28 nA and 34.41 nA, respectively, while Step 4 decreases further to 31.59 nA and 29.70 nA, gradually widening the gap between them.\nThese changes demonstrate that the influence of the first three pulses on the last pulse gradually diminishes as the pulse interval increases. The current response at the final time step transitions from being predominantly shaped by the collective effect of the preceding pulses to being more independently governed by the last pulse itself.\nThis shift in the dynamic response mechanism provides the system with new opportunities to learn and map image features from different perspectives.\nFig. R2 Mapping of the digital number 5 using different pulse intervals.\nThe 4-bit encoding mapping method was adopted, where each row of 28 elements is divided into 7 regions and mapped using 4-bit encoding. The original binarized data exhibits significant differences after mapping. Since the image is divided into 196 subregions, the data acquires a temporal dimension after undergoing the dynamic response transformation of the device. Different pulse intervals lead to distinct mapping patterns, with variations not only in the range of conductance values but also in the differentiated presentation of local features.\nTo facilitate understanding, we take the red-marked region in the figure as an example. The corresponding encoding for this region is \"1000\". It can be observed that at a time interval of 5 ms, the conductance of this region gradually decays over time and stands out clearly against the background. However, when the time interval increases to 10 ms, the features of this region weaken significantly. Further increasing the interval to 15 ms, the I at Step 4 becomes similar to the background, indicating a noticeable degree of decay. This shows that the mapping characteristics at this point differ markedly from those at 5 ms.\nAs the time interval continues to increase, this trend becomes more pronounced. At an interval of 20 ms, the I at Step 3 begins to lose distinctiveness; by 25 ms, the I at Step 2 also becomes blurred; and when the interval reaches 30 ms, the I at Step 2 is already difficult to clearly discern. Thus, the dynamic response mechanism can effectively process image features, enhancing the diversity and hierarchical richness of the data. When applied to image recognition tasks, the system can leverage the features exhibited during the dynamic response process as a reference, thereby improving the accuracy and robustness of image recognition.\nDue to the distinct boundary between the digits and the background, when encoding is used for mapping, the involvement of codes other than pure \"0000\" or \"1111\" tends to yield richer feature variations. The reason is that when only \"0000\" and \"1111\" are involved, the encoding merely reflects the maximum and minimum values of the I . In contrast, when boundaries are present, the system through its dynamic response mechanism can generate a greater diversity of feature expression possibilities.\nFig. R3 Mapping of the shoe using different pulse intervals.\nOn the Fashion dataset, although the overall situation is similar to that described earlier, the elements in this dataset are inherently denser. After binarization, compared to the handwritten digits dataset, it exhibits fewer edge features. Despite the presence of the dynamic response mechanism, due to the relatively similar structure of the elements within the dataset, the expressive differences remain rather limited. For instance, when the pulse interval was set to 5 ms, we applied the same method for mapping and obtained corresponding results. However, since most of the image content occupies the main area of the 28 × 28 pixel grid, the corresponding subset patterns are highly similar. When the time interval was increased to 10 ms, the output only showed a lightening of color, indicating that this is essentially a data scaling operation. Further increasing the interval to 15 ms resulted in continued fading of the color, but no significant new features generated by the dynamic response were observed. Thus, it can be concluded that the dynamic response provides only limited improvement in enhancing the feature representation capability for this dataset.\nFig. R4 Reservoir states of a physical ESN output for different pulse intervals. a digits dataset, b fashion dataset.\nWe utilize ambipolar transistors to construct a sparse matrix and a tanh function, which together form a high-dimensional, sparsely coupled nonlinear dynamical system. In this system, the sparse matrix defines the topology of information propagation and multi-scale memory pathways, while the tanh function introduces saturating nonlinearity, confining the system dynamics within a bounded range. The synergistic interaction between the two generates rich separated states and gradually decaying memory properties, thereby providing an analytically tractable dynamical foundation for temporal information processing.\nExperimental results demonstrate that our proposed physical dynamic preprocessing method effectively achieves high-dimensional nonlinear mapping of data. When the physically preprocessed data is fed into the physically implemented Echo State Network (ESN) system, the reservoir states naturally exhibit a balanced and separable spatial distribution. This effectively avoids common issues such as neuron saturation or degradation often observed in traditionally randomly initialized reservoirs, significantly enhancing feature discriminability and system interpretability. As shown in Fig. R4a and Fig. R4b, which present the spatial distributions of partial reservoir states for the handwritten digit dataset and the Fashion dataset, respectively, the state values are predominantly concentrated within the interval [-0.1, 0.1]. This indicates that the system operates in the approximately linear region of the tanh function, maintaining dynamic stability while retaining sufficient nonlinear expressive capability. Further analysis reveals that the physical ESN essentially projects the input mapping into a higher-dimensional feature space, after which effective classification can be achieved using only a simple linear classifier, with excellent performance. This aligns with the fundamental idea of Cover's theorem—data in high-dimensional spaces are more easily linearly separable—which structurally explains the superior performance of the physical ESN. Moreover, leveraging the dynamic response of ambipolar transistors for image preprocessing essentially integrates the dynamic characteristics of the physical system with nonlinear computational mechanisms, offering a viable pathway toward efficient and low-power neuromorphic computing.\nFig. R5 MNIST datasets accuracy by jointly utilizing six reservoir spatial states. a digits datasets, b fashion datasets.\nWe employ dynamic response systems with different parameters to map input images, forming multiple feature extraction pathways. Distinctive features are extracted from the reservoir states of each pathway, each reflecting the projection characteristics of the input image under different physical dynamic mechanisms. Subsequently, all multi-path features are concatenated and fused, and then fed into a simple classifier for joint recognition.\nThe experimental results demonstrate that this method effectively enhances classification performance. On the handwritten digit dataset and the Fashion dataset, classification accuracies reach 96.98% and 86.67%, respectively (as shown in Figures a and b). In particular, the 4-bit preprocessed data achieved through varying pulse widths, when further processed by a physical ESN, significantly improves classification accuracy. The underlying principle is that multi-parameter preprocessing generates diverse feature mappings. Compared to the original single image, these mappings express richer structural information from multiple dynamic dimensions, thereby enhancing feature discriminability and robustness.\nFig. R6 Lorenz attractor mappings under different pulse intervals along the x, y, and z axes. a x axes, b y axes and c z axes.\nWe model the dynamic response of the transistor as an intrinsic kinetic process. When the chaotic time-series signal of the Lorenz system, such as x(t), is normalized and applied as a gate voltage input, the charge transport and storage mechanisms within the device introduce a memory superposition effect of historical states. As a result, the real-time I reflects not only the current input but also the influence of past signals. In the experiment, pulse sequences with different time intervals (5-25 ms) were used for excitation and readout to probe the dynamic memory characteristics of the device. Fig. R6a presents the measurement results along the x axis: under short-interval pulse excitation (5 ms), the I response exhibits sharp peaks, indicating that the device retains significant memory of recent inputs. When the pulse interval is extended to 10 ms, the amplitude of the Ids response weakens noticeably. As the pulse interval further increases to 25 ms, the response further decays and the features become increasingly smooth. This phenomenon is also observed along the y and z axes, demonstrating that the memory effect of the device is time dependent: under short time intervals, the output is strongly correlated with historical signals; as the interval lengthens, this historical correlation gradually diminishes.\nBased on the above dynamics, the transient responses collected from the device construct a state sequence driven by nonlinear charge-carrier kinetics. This mechanism naturally performs nonlinear embedding and feature expansion of the input Lorenz signal along the temporal dimension, endowing the resulting data with inherent temporal structure and dynamic information. Consequently, this representation not only preserves the time evolution characteristics of the system but also provides a physically rich input space for subsequent nonlinear decoding and modeling.\nThe three-dimensional trajectory data (x, y, z) of the Lorenz attractor were normalized to match the dynamic range of the target hardware system. Subsequently, the normalized signals were linearly converted into corresponding analog voltage signals. Finally, through a transconductance circuit or an equivalent conversion module, the voltage signals were further transformed into current signals with explicit temporal evolution characteristics, serving as time‑series driving sources for subsequent processing. Fig. R6a, b and c respectively show the driven results of the x‑, y‑, and z‑axes under different pulses. These temporal current signals drive the physical device at their respective pulse intervals, inducing distinct nonlinear responses in the system.\nFig. R7 Reservoir spatial states corresponding to different pulse intervals for the x, y, and z axes.\nLeakage rate (forgetting factor) is the core of ESN for time series prediction tasks. When there is no leakage rate (α=1), the reservoir state is completely overwritten by a new nonlinear transformation at each time step, forming rapidly changing discrete dynamical states. By introducing a leakage rate (0<α<1), the state update becomes a combination of linear mixing of old and new information with nonlinear transformation, equivalent to applying low-pass filtering and inertial smoothing to the reservoir dynamics. This results in continuous state trajectories with adjustable memory time constants—smaller α values lead to slowly evolving long-term memory, suitable for handling slowly varying continuous signals.\nBy introducing the forgetting factor earlier, we endowed the system with both short-term and long-term memory behaviors. Consequently, its reservoir state space is continuous. With the introduction of the leakage rate (forgetting factor), the reservoir state update transitions from discrete jumps to continuous evolution, thereby acquiring dynamically tunable short- and long-term memory characteristics. This continuity causes the reservoir state trajectories to manifest as smooth, inertial flows in phase space rather than discontinuous jumps. Specifically, each dimension of the state vector becomes a dynamical variable with exponentially decaying memory, whose time constant is directly controlled by the leakage rate. This continuous dynamic not only better simulates the evolution processes of real physical systems (such as neural membrane potentials or circuit transient responses) but also naturally encodes temporal derivative information of the input signals into the reservoir states.\nBased on our design modifications, we have redefined the role of the physical ESN, moving away from “simulating the evolution of physical systems” and instead emphasizing its core function as a nonlinear transformer. Under this framework, the system’s temporal characteristics are no longer explicitly generated by the physical ESN’s own state dynamics but originate from the preprocessing coupling between the original time-series signals of the Lorenz attractor and the device’s dynamic characteristics.\nSpecifically, the chaotic time-series signals (x(t), y(t), z(t)) of the Lorenz system are first applied as excitation inputs to an ambipolar transistor. The device’s inherent nonlinear transconductance characteristics, carrier recombination dynamics, and the electron-hole competition mechanism regulated by the gate voltage collectively form a complex dynamic preprocessing layer. This process transforms the original deterministic Lorenz trajectories into a set of high-dimensional I response sequences that depend on the physical properties of the device. Although these I data originate from a dynamic system, they constitute a \"device-dynamic-feature-encoded\" static data stream relative to the input of the physical ESN.\nSubsequently, these preprocessed I sequences are fed into a static physical ESN. Here, \"static\" means that the ESN reservoir itself does not possess time-varying state update mechanisms (i.e., no recurrent connections or leaky integration are introduced). Its role is purely to perform spatial nonlinear projection and feature expansion on the input high-dimensional I signals, forming a fixed-dimensional reservoir state space.\nAs shown in the figure, the output reservoir states are structurally static high-dimensional point sets, where each point corresponds to a feature vector at a specific time after multiple nonlinear transformations.\nAlthough the final reservoir states are static and do not contain explicit temporal recurrence, the critical time dependencies have been effectively encoded into the statistical structure and sequential relationships of the I sequences during the front-end \"device dynamic preprocessing\" stage. Therefore, this static high-dimensional state space still carries implicit dynamical information about the evolution of the Lorenz system. By training a simple linear readout layer or a lightweight nonlinear decoder, the system can successfully reconstruct and predict the future states of the Lorenz attractor from these static yet information-rich features.\nThis demonstrates that even with a static reservoir architecture, as long as the front-end preprocessing effectively injects temporal dynamic information, the entire system can still accomplish complex chaotic time-series prediction tasks.\nFig. R8 Prediction results of the Lorenz attractor by the physical ESN under different pulse intervals.\nBy preprocessing raw time-series data through dynamic characteristics and then solving it via static nonlinear transformations, excellent predictive performance can still be achieved. The dynamic preprocessing stage encodes the temporal evolution characteristics and physical constraints of the system into a high-dimensional feature space, while the subsequent static nonlinear transformation focuses on extracting separable patterns from these features. Specifically, the dynamic processing converts temporal dependencies into spatially distributed statistical features through the transient responses of devices or analog circuits. Subsequently, the static nonlinear layer adjusts the dimensions and separates the patterns of these features, ultimately enabling the reconstruction of system dynamics through linear readout. This \"dynamic encoding–static decoding\" architecture not only reduces reliance on recurrent computational resources but also demonstrates that temporal information can be transformed into spatial features through preprocessing, thereby enabling the modeling and prediction of dynamic systems within a static framework.\nAs shown in the figure, despite the clear reconstruction of the Lorenz attractor’s topological structure across different data collection time intervals, the predicted trajectories still exhibit subtle distortions and local deviations. This phenomenon can be attributed to the modeling strategy of \"predicting the state at t+2 using the state at t\". While the longer prediction step validates the model’s generalization ability, it also amplifies the extreme sensitivity of chaotic systems to initial conditions (i.e., the butterfly effect). In the physical reservoir composed of ambipolar transistors, this sensitivity is further modulated by device non-idealities, leading to cumulative prediction errors manifested as trajectory jitter or slight morphological distortions of the attractor in phase space.\nThese results are comparable to those obtained by introducing a leakage factor, at least indicating the effectiveness of our approach. Furthermore, when reservoir features extracted from five different pulse intervals are used for prediction, the results become significantly smoother.\nFig. R9 NRMSE under different pulse intervals.\nThe NRMSE of the predictions across the three Lorenz system dimensions reveals distinct characteristics regarding the influence of pulse interval duration. On the x‑axis, the total prediction error (i.e., the combined error across all intervals) is the highest, with an NRMSE value of 0.23%, indicating that the x‑dimension dynamics are the most challenging to reconstruct under the current modeling setup. On the y‑axis, predictions for the 5 ms, 20 ms, and 25 ms intervals exhibit relatively higher errors, with NRMSE values of 0.745%, 0.846%, and 0.862% respectively, suggesting that the y‑dynamics are particularly sensitive to both very short and long pulse intervals. In the z‑dimension, the highest NRMSE occurs at the shortest pulse interval (5 ms) with a value of 0.799%, reflecting that the faster sampling may inadequately capture the slower-scale variations characteristic of the z‑variable in the Lorenz system. Overall, these results highlight that prediction accuracy is not only dimension-dependent but also closely tied to the temporal sampling strategy, with excessively short or long intervals leading to increased reconstruction errors in specific dynamical components.\nTo better align with the reviewers’ feedback and emphasize the core contribution, we have revised the manuscript title to “Physical echo state network based on the nonlinearity and dynamic response of ambipolar heterostructure transistors”. We believe this computing paradigm is novel and offers greater suitability and potential for time-series and image processing tasks compared to conventional reservoir computing approaches.\nRevision in the manuscript:\nFig. 3. Dynamic-mechanism-based physical ESN for MNIST handwritten-digit recognition. a Four-bit encoding built on the dynamic response. b Realisation of “0,0,0,0” and “1,0,1,0” codes under different pulse intervals. c Mapping of the MNIST dataset achieved by varying pulse intervals. d Corresponding reservoir-state space generated by the physical ESN for these mappings and e Confusion matrix obtained with five combined pulses, showing the final classification performance.\nThe MNIST datasets are used for classification39-41. These images are all 28 × 28 size and are fed into the physical ESN in the same way as above. 60,000 images are used as the training set, while the remaining 10,000 images serve as the test set. Thanks to the inherently excellent nonlinearity of ambipolar transistors, even when performing direct image straightening, the data can be mapped to a high-dimensional space. This allows the readout layer to easily interpret these reservoir states, enabling a single-step solution with high accuracy. The ESN achieves an accuracy of 95.27% on the grayscale MNIST handwritten digit dataset and 86.10% on the MNIST fashion dataset(Supplementary Fig. S19-S23). This is not the ultimate limit of the architecture. To achieve higher accuracy, scaling up the sparse matrix size is required. Given its already high accuracy, further improvements would necessitate additional devices, introducing inherent limitations. Another method involves prepossessing the datasets with the dynamic response of the device and then using the prepossessing datasets to complete the recognition task. Since the dynamic response can provide features for the datasets, effectively adding dynamic nonlinearity, it can also yield favorable results. We further binarize the MNIST dataset and use 4-bit encoding to map the images, followed by processing with a physical ESN for recognition. Different pulse intervals in the 4-bit encoding are employed to impart more features (Supplementary Fig. S24 and S25). Information encoding is first established through a four-bit scheme (Fig. 3a). Different types of nonlinearity are achieved by adjusting the interval pulses. For example, when the interval is increased, the responses produced by the \"0,0,0,0\" and \"1,0,1,0\" encodings are different (Fig. 3b). This gives us further opportunities for exploration. When images are processed using encodings with different pulse intervals, different image mappings are obtained (Fig. 3c). When these mappings are processed by the physical ESN, they produce different reservoir space states (Fig. 3d). Both the digital dataset and the fashion dataset exhibit discrete reservoir states (Supplementary Fig. S26). However, their time attributes are incorporated due to preprocessing with different pulse intervals. This results in the potential for creating better performance opportunities. When the mappings converted from six pulse intervals are read out by the physical ESN, the accuracy reaches 96.98% (Fig. 3e). Similarly, on the fashion dataset, after dynamic processing, the accuracy also increases relatively to 86.67% (Supplementary Fig. S27), higher than 86.10%. Therefore, introducing nonlinearity and dynamic mechanism through physical methods in MNIST dataset processing may serve as an effective approach to enhance performance while reducing computational costs.\nFig. 4. Lorenz-attractor prediction task realized by dynamic-response-based physical reservoir computing. a Schematic of the approach: device dynamics encode and forecast the Lorenz attractor. b-d Nonlinear mapping of Lorenz-attractor trajectories along x, y and z axes using the intrinsic dynamics of the single device (V =1V, the Lorentz attractor curve is normalized and then converted to a 0 to +10V V signal g using DAC). e Prediction NRMSE versus pulse-interval variation. f Forecasted Lorenz-attractor time traces under different operating modes.\nFor neuromorphic computing systems, the capability to forecast temporal and spatial sequences is essential for ESNs42–44. To evaluate this, we apply our ESN to predict the Lorenz attractor, using data at time step t as input and at t+2 as target. This task represents a second key challenge in time‑series prediction, an area where hardware‑based architectures are actively studied. We explore two training strategies. The first adopts the forgetting‑factor method commonly used in ESNs (Supplementary Fig. S28), where the factor α substitutes for explicit dynamic mechanisms. By tuning α, the physical ESN achieves normalized root‑mean‑square errors (NRMSE) of 0.165%, 0.28%, and 0.206% along the x‑, y‑, and z‑axes, respectively. The forgetting factor ensures continuity in reservoir states across time (Supplementary Figs. S29–S31). The second approach directly embeds temporal attributes by preprocessing the data through the dynamic response of the physical device before training (Fig. 4a). Here, training proceeds with a forgetting factor of 1, equivalent to a static mapping setup. The device’s dynamic response performs nonlinear transformation on the Lorenz attractor, and the resulting time‑varying Ids serves as input to the ESN (Fig. 4b–d). By varying the sampling intervals during preprocessing, different dynamic features can be extracted for prediction. Though the reservoir states under this method appear discrete, they retain encoded temporal information (Supplementary Fig. S32). This preprocessing strategy yields NRMSE values of 0.165%, 0.28%, and 0.206% on the x‑, y‑, and z‑axes, respectively (Fig. 4e)—performance that closely matches, and in some dimensions slightly exceeds, that of the forgetting‑factor method. Consequently, the system achieves low prediction error and generates smooth Lorenz attractor trajectories (Fig. 4f).\nRevision in the SI:\nFig. S24 4-bit encoding with varying pulse intervals. a 5 ms. b 10 ms. c 15 ms. d 20 ms. e 25 ms. f 30 ms.\nFig. S25 Mapping of the shoe using different pulse intervals.\nFig. S26 Reservoir states of a physical ESN output for different pulse intervals. a digits dataset. b fashion dataset.\nFig. S27 fashion datasets accuracy by jointly utilizing six reservoir spatial states.\nFig. S32 Reservoir spatial states corresponding to different pulse intervals for the x, y, and z axes."}, {"round": 5, "reviewer_comments": {"reviewer_1": "My concern especially regarding the processing of temporal information using the suggested device system is well addressed. The manuscript can be accepted.\n(Remarks on code availability)", "reviewer_3": "The authors have made substantial modifications, transforming the previous static input into a temporal one, thereby enabling the reservoir to exploit its temporal properties. I believe this significantly improves the paper, which may be accepted as is now.\n(Remarks on code availability)"}, "author_response": "Response: We sincerely appreciate your review comments, which have significantly enriched our content."}]

Decision letter