Ferroelectricity in graphene nanoribbon devices enabled by collective water molecule dynamics

Domain: physics
Published: 2025-12-09
Editorial outcome: minor_revision

Abstract
Water is omnipresent in nanoscale systems, yet its collective dynamics and impact on emerging electronics remain poorly understood. Here, we investigate the role of water molecule dynamics in the ferroelectric response of graphene nanoribbon devices. Our findings demonstrate that the collective dynamics of water molecules stabilize the ferroelectric effect. We find that a minimum bi-layer thickness is required for the temperature stability of the ferroelectric effect. In contrast, mono-layer ribbons show a 70% shrinkage of the hysteresis window between 120 and 400 K. Using a combination of electrical transport measurements and molecular dynamics simulations, we conclude that water molecules bridging between graphene nanoribbon layers stabilize the formation of water clusters via intermolecular Coulomb interactions, driving a robust ferroelectric behavior and remnant polarization observed at the device level. This work lays the foundations for exploiting water dynamics in next-generation ferroelectric heterostructures, with direct implications for neuromorphic computing and memory devices.

Introduction
Water plays essential roles in a wide range of phenomena, from materials and electronics to life itself1,2. Understanding the interaction of water with low-dimensional materials has proven essential for nanofluidics3, energy storage4,5, water splitting6, purification7,8, and water-assisted ferroelectricity9–11. Water plays a critical role in tuning material properties at the nanoscale. It turns a challenge into opportunities for nanoelectronics by modifying electronic transport, charge distribution, and doping, making it a critical factor for applications in neuromorphic computing and molecular-level devices, especially when considering the impact of the edge-adsorbed species on graphene nanoribbons (Gr NRs)12–16.

Graphene nanoribbons with a high edge-to-surface ratio provide an excellent platform to study the influence of water dynamics at the edge-water molecule interface. Previous work showed that water molecules adsorbed at graphene edges in single-layer microdevices encapsulated by hexagonal boron nitride (hBN) can induce ferroelectricity10, a phenomenon also observed in nanoribbon-based devices11. This phenomenon originates from the dipolar nature of water that can result in switching between two polarization states under external electric fields. However, the exact mechanism behind water-induced polarization in Gr NR devices, the impact of field-related dynamics of edge-adsorbed water, the influence of graphene thickness on water anchoring, temperature-dependent stability, and the interplay of these factors on ferroelectricity remain poorly understood. These knowledge gaps limit the development of next-generation devices based on water-induced ferroelectricity, including applications in radio frequency switches, neuromorphic computing, and memcapacitors17–21. This study builds up on the existing knowledge of water-induced ferroelectricity in graphene edges10,11. It brings forward new experimental results demonstrating how collective water dynamics are critical to this effect through graphene layer number and temperature dependence. These results are supported by molecular dynamics simulations describing the collective behavior of water at graphene edges. Our twofold experimental and computational approach shows that this ferroelectric effect at the device level in Gr NRs integrated into field-effect transistors (FETs) strongly depends on the graphene layer number, requiring a minimum bi-layer thickness to sustain temperature-stable ferroelectricity. We elucidate the physical mechanism behind these results as driven by the collective behavior of water that induces edge anchoring and interlayer bridging by molecular water necks that enhance intermolecular Coulomb interactions. By investigating the interplay between external fields, temperature, layer number, and adsorbed water at graphene edges, we open the door to future nanoelectronics and sensor applications leveraging the ubiquitous nature of water confined at the nanoscale22–29.

Results
Observation of the hysteretic response in oxygen-terminated graphene nanoribbon devices
Figure. 1a depicts the schematic of a Gr NR FET on hBN. Para-hexaphenyl (6P) organic nanoneedles self-assembled on hBN-Gr-hBN heterostructures serve as self-aligned masks11,24,30,31. Reactive ion etching (RIE) selectively removes exposed graphene to produce Gr nanoribbons, as illustrated in Fig. 1b, c; detailed in the Methods section and Supplementary Fig. S1. The resulting nanoribbons have an average width of 20 nm, as estimated from high-resolution atomic force microscopy (AFM) images (Supplementary Fig. S2). Alternative to organic nanostructures, recent studies have shown that directed self-assembly of block copolymers can produce complex arrays of highly ordered sub-10 nm wide nanopatterns32–34, that can be used to fabricate 2D material nanoribbon FETs35.

The top and bottom encapsulation with hBN insulates graphene from potential charge-trap sources, ensuring that external interactions primarily occur at the nanoribbon’s edges (Supplementary Fig. S3). Organic nanostructure growth is done near thermodynamic equilibrium conditions, meaning that even minor changes in the substrate can significantly impact self-assembly. The interfacial bubbles introduce bi-axial tensile strain in the top hBN layer, and the organic nanostructure self-assembly is hindered in such strained regions. As a result, the nanoribbon network effectively “avoids” these bubbles (Supplementary Fig. S3d–g).

The nanoribbon edge termination is controlled by the choice of a precursor gas used for RIE. Using oxygen (O2) or sulphur hexafluoride (SF6) results in oxygen- or fluorine-terminated edges, respectively10. Electrical transfer results show that devices with oxygen-terminated edges exhibit pronounced and stable hysteresis, with switchable p- and n-doping characteristics depending on the gate bias (VSG) sweeping direction (Fig. 1d). The electrical response comparison between the nanoribbons and the starting flake FETs is shown in Supplementary Fig. S4. The hysteresis window (VH) is defined as the difference between two charge neutrality points (CNPs), which correspond to the maxima in the total two-terminal device resistance Rtot(VSG). In contrast, F-terminated ribbons show negligible hysteresis due to their hydrophobic nature (Fig. 1d, blue plots), as previously observed for graphene flakes10. Given that oxygen-termination enables water molecule adsorption, our focus will be on oxygen-terminated graphene nanoribbons and their water-induced ferroelectricity10,36,37. To ensure water molecule adsorption on the oxygenated ribbon edges, Gr NR devices were exposed to ambient conditions (relative humidity of 22–26% at 298 K) for at least 30 min following the RIE step. This procedure was also applied to the SF6-etched devices to maintain comparability between the two datasets.

We have used a total of 22 devices, comprising 16 O2 RIE-treated and 6 SF6 RIE-treated devices. To illustrate sample-to-sample variation in performance parameters, we present peak apparent linear mobility (µlin) statistics in Supplementary Fig. S5. Furthermore, µlin(VSG) curves for both SF6 and O2 etched devices were compared before and after the RIE step, see Supplementary Fig. S6. The assembly of the nanostructure network varies between samples, impacting the scattering of ribbon widths. Consequently, the VH values show a scatter of about ±25% among different devices under comparable operating conditions. This scattering of VH values is illustrated in Supplementary Fig. S7, as a function of the layer number and the total sweeping range of the gate voltage (VSG-tot). Further discussion regarding the relationship between ribbon width and VH is provided later in the text from the perspective of molecular dynamics (MD) simulations.

The details of the mechanism behind the hysteretic characteristics presented in Fig. 1d are known, as Caridad et al. proposed that single molecules switch between the two states10. The external gate fields disrupt the equal probability of water arrangements while inducing a torque on the water adsorbed at the edges. The torque results from Coulomb forces acting on the water dipole, allowing it to cross the energy barrier between the two states ca. 25–40 meV10,38,39.

Interestingly, Rtot(VSG) hysteresis and its thermal stability exhibit profound differences depending on whether the oxygenated edges of the ribbons are mono-layered or multi-layered. Fig. 1e, f compares the temperature dependence of the normalized hysteresis window for mono-layer and multi-layer (2–5 layers) Gr NR FETs over a temperature range of 120–400 K. The VSG sweeping range was fixed between −40 V and +40 V for each VH value, with a constant sweeping rate of 2 Vs−1.

In mono-layer Gr NR FETs, the hysteresis window decreases as temperature increases, although significant values persist across the whole temperature range. In contrast, multi-layer devices with two or more graphene layers (L ≥ 2) show a ferroelectric effect with negligible variations across the investigated temperature range (see Fig. 1f). For multi-layer devices (2L to 5L), we did not observe any significant differences in the response of the Gr NR FETs with different layer thickness. However, there is a trend indicating that both VH and µlin tend to decrease as the number of layers increases (Supplementary Figs. S5 and S7).

To estimate the height of the barrier observed in the experiments, we fitted the VH(T) curves (illustrated by solid black lines in Fig. 1e, f) using an analytical model based on the Arrhenius equation. This model considers the activation energy (Ea) or the energy barrier between the two polarization states, expressed as the fitting constant in: VH(T) = VH0 / (1 - exp(-Ea/(kBT))). Here, kB is the Boltzmann constant, T is the temperature, and VH0 scales the fit to the value obtained at the lowest temperature. For the monolayers, we obtained an activation energy of Ea = (17.45 ± 11.60) meV. The fits for the individual VH(T) curves of the mono-layer GrNRN FETs can be found in Supplementary Fig. S8. In contrast, the data for the multilayer devices in Fig. 1f shows an energy barrier higher than 150 meV.

These observations led us to hypothesize that a single water molecule switching between the sides of the graphene plane is responsible for the ferroelectric response in mono-layer devices, leading to a temperature-dependent behavior. When the thermal energy of the molecule exceeds the switching energy barrier, the molecules can thermally transition between states, which explains the temperature-dependent results in Fig. 1e. Additional results regarding the dependence of the hysteresis window on the VSG sweeping range can be found in Supplementary Fig. S9.

Stability of the hysteretic response: excluding charge trap origins
The direction of the observed hysteresis could originate from the charge traps in the gate dielectric layer. However, the same hysteresis direction could also result from the edge dipoles, as illustrated in Supplementary Fig. S10. To further verify the nature of the observed hysteresis, we performed several experiments to rule out the influence of extrinsic effects such as impurities and charge traps. Since trapping can occur on the timescale of seconds40,41, we used varying sweep rates (from 0.5 Vs−1 to 20 Vs−1) to evidence its effects on VH. Figure. 2a shows that VH as a function of the sweeping rate remains largely unperturbed at 298 K and 77 K under forward and backward gate sweeps, suggesting that the hysteresis is not dominated by charge traps. Furthermore, measurements at 4.2 K in Fig. 2b show that the hysteresis persists, ruling out contributions from thermally activated traps42,43. Interestingly, near the charge neutrality point, the resistance exhibits aperiodic fluctuations with high reproducibility, indicating the opening of a transport gap. While such fluctuations could be attributed to quantum interference phenomena44,45, Chen et al.46 showed that such insulating states can be attributed to the adsorption of water molecules on the nanoribbon edges, generating strong electric fields that modulate the NRs’ bandgap.

To further demonstrate the impact of water molecule adsorption, the devices were vacuum annealed at 520 K for 60 min at a pressure of 8·10−3 mbar. Subsequent electrical measurements were performed at 300 K, avoiding any ambient exposure. As shown in Fig. 2c, d, the hysteresis window VH drastically decreases after annealing, consistent with the desorption of water molecules from the nanoribbon edges. Remarkably, the hysteresis is restored upon exposure to ambient conditions (relative humidity of 22–26% at 298 K), providing strong evidence for the critical role of water adsorption in the ferroelectric response. In addition, the unintentional p–type doping evidenced in Fig. 2c implies that only positive VSG values are required to compensate for built–in charge transfer doping before the total field direction can reverse and switch the water dipoles orientation (Supplementary Fig. S11). The dependence of VH on the vacuum conditions is detailed in Supplementary Fig. S12, where we observed a 40% increase in VH when going from vacuum conditions to ambient pressure at 300 K under controlled relative humidity of 40%. Furthermore, VH remained constant at pressures below 101 mbar.

Molecular dynamics simulations of the oxygenated and water-terminated graphene nanoribbons
We performed MD simulations to understand better the role of water molecules in inducing electric fields and to elucidate the mechanism behind the temperature-dependent hysteresis differences between mono-layer and multi-layer Gr NR devices. The tendency of the water molecules to bind at oxygenated edges and form clusters39 suggests that a collective behavior stabilizes the molecules, resulting in a temperature-independent effect. Two conditions must be met to achieve such stable ferroelectric ordering of water molecules on Gr edges: (i) a low thermal barrier that allows a single polar molecule to switch between states, and (ii) a cluster large enough to be stabilized by intermolecular Coulomb interactions. While the first condition is well established, we have designed molecular dynamics simulations to investigate the second condition. The simulations implicitly consider the hydrogen bonds within the dynamically calculated local charges between water molecules and the oxidized graphene edge. Moreover, we also calculate the electric field, which is a direct consequence of the molecule’s polarization, as shown later. Our calculations agree with the model presented in ref. 10, and provide new insights into the importance of the collective behavior of water at multi-layered edges.

To conceptualize the mechanistic difference between the two systems, we can think of both mono- and bi-layers as a bi-stable system with two states and a barrier between them. The collective interaction that stabilizes the cluster in one state could essentially be viewed as the deepening of one minimum, depending on how many molecules are in that minimum. This introduces effectively a discrepancy in the rate of returning to the lower occupied state by thermal fluctuations. The process is illustrated in Supplementary Fig. S13.

Therefore, the existence of water interaction and the size of the cluster play a key role in achieving the temperature-independent hysteretic effect. Naturally, there is a limit to how large a cluster can be while remaining anchored at the edge. As the cluster size increases, the induced dipole field also grows. At the same time, the energy required to detach the molecules furthest from the edge decreases (Supplementary Fig. S14). These two effects compete. In the presented simulations, we have limited the amount of water molecules per unit length of the edge to below one-half of the limit determined by the cluster stability in strong external fields at 300 K (more details in the “Methods” section).

The MD simulations of mono-layer graphene confirm the predicted energy barrier of 25−40 meV for switching a single water molecule10,38. In contrast, the multi-layer system shows a fundamentally different behavior. The layer separation in graphite allows water molecules to form a bridge between two oxygen-terminated edges. These bridging molecules remain stable and do not react to external electric fields, yet they promote the formation of a water cluster around them. MD simulations show that anchored water molecules within the cluster can sustain high electric fields (up to 5 V nm−1) and temperatures up to 500 K before desorbing from the edges. The cluster surrounding the bridging molecules is mobile in the electric field, and its collective behavior helps stabilize the structure, resulting in a hysteresis window that is essentially temperature-independent. This collective self-stabilizing behavior of polar objects resembles that previously observed in colloidal systems and molecular motors47,48. Our calculations do not rule out the formation of water clusters in mono-layers. In mono-layers, flipping of individual molecules across the basal plane of the ribbons is frequently observed, even in the presence of clusters. However, such transitions of individual molecules are essentially absent in the bi-layer case. The distribution of water molecules along the edges in both the 1L and 2L cases is presented in Supplementary Fig. S15.

Expanding on our findings, we investigated how varying electric fields affect water clusters using MD simulations. An external homogeneous electric field generated by infinite planes was applied to 1L and 2L Gr NRs, as described in the Methods section. The fraction of polarized molecules above (n+) and below (n-) the graphene layer quantifies the ferroelectric effect arising from molecular switching. We observed charge bi-stability and switching between two states under a cyclical change of the external field. The resulting difference between n+ and n- per unit length of the NR edge is presented in Fig. 3a, b, comparing the hysteresis loops obtained for mono- and bi-layer Gr NRs. It is important to note that in the experiment, the electric field source was generated between the ribbon and the flat surface (gate electrode). In contrast, in the simulations, the electro-neutral system with the ribbon and water molecules was placed between two uniformly charged planar electrodes, generating a homogeneous external out-of-plane field. Further, in the simulations, the electric field applied was over an order of magnitude higher than in the experiment when considering a parallel capacitance model. This was done to reduce simulation time since the electric field strength exponentially prolongs the switching time of the model system. Field enhancement is expected in our experiments due to fringing capacitance effects10. An estimate of the enhancement and the field profiles based on a finite elements model of the device structure is provided in the Supporting Information, Supplementary Fig. S16. We anticipate about one order of magnitude enhancement of the out–of–plane field component due to the fringing capacitance effects, reaching up to 2 V nm−1 within the experimentally achievable biasing range.

For a mono-layer ribbon, a weak hysteresis in n+-n- is observed as a function of the external electric field. Moreover, this hysteresis decreases with increasing temperature (as shown by the orange circles in Fig. 3c), indicating a low energy barrier between the states. This phenomenon can be attributed to the absence of a bridging water molecule found in multi-layer systems. Our experiments with mono-layer Gr NRs showed a more robust hysteresis, especially at elevated temperatures. The discrepancy between the modeled and experimentally observed mono-layer systems may be attributed to more complex interfaces present in the actual experiments that are not captured in the simulations. The model considers a free-standing graphene ribbon, while the experiment involves hBN encapsulation, which likely creates complex electrostatic environments that could also partially enhance the collective behaviour of the adsorbed water clusters. In addition, the model system has a considerably shorter temporal evolution (4 ns) in response to the external electrostatic perturbations compared to the experimental timescale of over 100 ms. This shorter simulation timescale may prevent the model from adequately relaxing within the achievable 4 ns window, especially at temperatures below 250 K. In contrast, the experimental system has more time to achieve polarization and to experienced ecoherence of the polarized state once the external field is turned off. Consequently, direct quantitative comparisons of the remnant field values between the experiment and the model are not feasible. Nevertheless, the model effectively captures the trends in temperature and layer number observed in the experiments, allowing us to elucidate the physical mechanisms behind the ferroelectric behavior of Gr NRs devices.

For a bi-layer NR, a prominent hysteresis was noted in the MD simulations (Fig. 3b), and the opening of the hysteresis loop was found to be temperature-independent in the probing range from 250–450 K, as shown by purple squares in Fig. 3c. The field lines and field strength acting on a 2L nanoribbon generated via MD simulations are presented in Supplementary Fig. S17. As the electric field decreases, some molecules remain stationary until the field polarization switches, causing these molecules to migrate to the other side with respect to the nanoribbon’s basal plane.

The water cluster at the edge of a 2L NR is stable in the external electric field up to 4.1 V nm−1, after which the molecules were observed to be removed from the cluster by the external field at 300 K. n+-n- hysteresis presented in Fig. 3d, uses the external field range up to the limit of the desorption of the molecules from the cluster. In this case, a clear indication of saturation is observed. Such high fields probed in these MD simulations are not reachable in the experiment due to the SiO2 and hBN field breakdown. To visualize the ordering of the water cluster, we follow its evolution along a down-sweep of the field in the hysteresis shown in Fig. 3d. Figure 3e–g represents snapshots of the edge segments for a 2L Gr NR system with the corresponding charge distribution and induced electric field lines (below). The system consists of both externally applied fields (not shown) and induced fields with opposite polarities in the z-direction. Figure 3e presents the case of saturation with an external field of 4 V nm−1 applied in the z-direction. An induced charge redistribution is observed due to the molecular polarization at the edges, which in turn creates an effective electric dipole moment. This changes the local field over the whole graphene bi-layer profile. The effective electric field is, therefore, a sum of the two fields, and as the external field decreases, they could cancel each other in parts of the graphene plane. Fig. 3f displays the configuration where the induced electric dipole points roughly along the basal plane of Gr NR, i.e., representing one of the coercivity points of the out-of-plane dipole field hysteresis. Consequently, the induced electric field component perpendicular to the ribbon will be negligible. Finally, for −4 V nm−1, we observe a complete reversal of the induced electric field (Fig. 3g). As the field strength increases, the migration and molecular polarisation repeat in the opposite order, generating a hysteresis loop. The bonding between the polar molecular ensemble and the graphene edge, together with the intermolecular Coulomb interactions, should be strong enough to prevent the external electric field from tearing off molecules from the cluster (Supplementary Fig. S14b). The details of the large-scale atomistic model are presented in Supporting Information Supplementary Fig. S18.

The width of the ribbons significantly impacts the strength of the dipole-induced fields. Through MD simulations of polarized dipoles, for both 1L and 2L ribbons, we have extracted the dipole-induced fields perpendicular to the basal plane of graphene and their dependence along the direction of the ribbon width. As shown in Supplementary Fig. S19, the field induced by a single polarized edge follows the theoretically expected x-2 dependence. We hypothesize that an optimal width of about 10 nm maximizes the intensity of the hysteresis effect, arising from two competing processes. For ribbons that are too narrow, the field induced by one row of polarized edge dipoles begins to affect the dipoles at the opposite edge, resulting in an antiparallel state. For widths below 10 nm, the field from one polarized edge is expected to impact the other edge, overcoming external fields for widths below 5 nm. In the other extreme, as the ribbon width increases, the electrical response begins to show more bulk-like properties of graphene, resulting in a reduction of hysteresis. Our MD simulations suggest that for widths above 100 nm, the contribution from the edges to the electrical response becomes negligible.

Probing the remanence of the water dipole-induced fields in graphene nanoribbons
To confirm the pronounced hysteresis predicted by the MD simulation, we tracked the drain current (ID) evolution of Gr NR FETs without an external gate field after pre-biasing the devices into the n+ or the n- state. Figure. 4a presents a hysteretic transfer curve of the water-terminated Gr NR FET, where the range of the gate voltages was chosen so that one of the charge neutrality points is at VSG = 0 V49. The asymmetric VSG sweep approach was used for CNP position adjustment (Supplementary Fig. S20). This way, the difference in ID for the n+ and the n- states is maximized. The two ID(VSG = 0 V) states can be defined as “0” and “1”, corresponding to the low or the CNP state and the high current state, schematically illustrated for a bi-layer in Fig. 4b.

By sweeping from VSG = 0 V to one end of the range and returning to VSG = 0 V, the system is set either to the ID “0” or “1” state. From this point on, the ID is recorded as a function of time without an externally applied gate field. The dipole-induced remanent field was observed in both mono- and multi-layer Gr NR FETs. The results comparing a mono- and a multi-layer Gr NR FET are shown in Fig. 4c, where the normalized difference between ID1(t) and ID0(t) was tracked for 600 s at 300 K and under low vacuum. To characterize the decay of the field, the curves were fitted to C∙exp(-t τ-1) where C stands for a scaling constant and τ for a time constant.

The remanence of the field for a mono-layer is less pronounced, and the current difference between the two states decays faster than in the multi-layer Gr NR FETs (Fig. 4c). The time constants indicate that the 1L Gr NR devices lose 90% of the initial ID difference between the “0” and “1” states in about 30 minutes and multi-layer Gr NR FETs in about 90 minutes. Further, a faster decaying component with τ = (64 ± 32) s, responsible for about 15–20% of the initial drop of the ID difference, was noticed in all devices and attributed to weakly bound molecules. Gr NR FETs with L ≥ 2 require several hours of relaxation to have the initial transfer curve sweep starting from the depolarized state of the water molecules, which is consistent with the extracted time constants.

Methods
2D material heterostructures
Crystals of 2D materials were mechanically exfoliated from bulk material and transferred onto a 290–300 nm SiO2/Si substrate using commercially available Nitto tape. The substrates were chemically cleaned by sonication in acetone and IPA, as well as plasma cleaned in a mixture of oxygen and argon. Mono-layer and multi-layer thick flakes were identified via optical contrast and Raman measurements. Mono-layer graphene crystals were transferred onto hBN substrates using deterministic dry transfer methods. The transfer was achieved using metal-coated SiNx membranes rather than the conventional polymer stamp transfer techniques50. Graphite electrodes and top hBN capping were transferred via a PDMS-based deterministic dry transfer method.

Organic nanostructures
Organic molecules were grown by hot-wall epitaxy11. As a source, parahexaphenyl (6P) was used. The molecules were sublimated at 493.15 K under 2·10−6 mbar pressure in a vertically mounted semi-closed quartz tube. The walls of the tube were kept at 513.15 K. The reactor tube is closed by the sample holder that is kept at a temperature of 383.15 K. To ensure minimal temperature variation between the subsequent growths, the system was left to stabilize for 45 min after all three growth temperatures were reached, yielding under ±1 K variation during the growth. The samples were exposed to the 6P vapor for 4–6 min, resulting in an equivalent of 0.8 ML (mono-layer) – considering layer-by-layer growth of 6P on the surrounding SiO2/Si. An example of the nanoribbon network (AFM topography map) after the etching step is presented in Supplementary Fig. S2.

Reactive ion etching
The etching parameters were adopted from the reference11. Reactive ion etching (Oxford Instruments PlasmaLab 100) was done in two steps; first employing SF6 plasma (50 sccm, 80 W, 8-10 s) to remove the top hBN layer, followed by oxygen plasma (50 sccm, 80 W, 4–10 s) to remove the graphene layers and to oxygenate the nanoribbon edges. Alternatively, prolonged SF6 etching (20 s) was used to etch through the whole stack, followed by brief oxygen etching to oxygenate graphene nanoribbon edges.

To ensure that the water molecules can populate the oxygenated ribbon edges, upon etching, all devices were exposed the ambient air, with regulated temperature and relative humidity (298 K and 22–26%). The waiting time of at least 30 minutes was applied under these conditions for each device in between the RIE step and the electrical measurements under vacuum. Further, no significant difference was observed between the freshly etched devices (with the 30 min holding time) and the devices that were tested, even several days after the RIE step, also stored under ambient conditions.

Electrical characterization
Electrical characterization of the flake and Gr NR FETs was done using a Keithley 2636A Source-Meter attached to the Instec vacuum probe station. The samples were contacted with Au-coated Ti electrical cantilever microprobes. Variable temperature (77 - 425 K) electrical measurements were performed using a liquid nitrogen flow cryostat with the devices mounted on a silver plate for thermal uniformity. The temperatures were monitored via Instec’s mK2000 temperature controller connected to the probe station with a resolution of 0.01 K. To ensure stable and reproducible electrical measurements, upon reaching the targeted temperature, the system was left to thermalize for 10 min prior to the electrical sweeps. This ensures that during the electrical measurements, the variation of the temperature is below ±0.5 K. Cryo (4.2 K) measurements were performed using a closed-cycle He cryostat and Lake Shore cryotronics temperature controller. All electrical measurements carried out in the Instec probe station (temperature range of 77–673 K) were performed at (2.2 ± 0.3)·10−2 mbar or as specified in the figure captions. The electrical measurements carried out at 4.2 K were done under the pressure of (4.0 ± 0.1)·10−6 mbar. To ensure stable low temperature, a holding time of 60 minutes was applied between reaching the targeted temperature and the start of the measurements. This ensured that the fluctuations were below ±0.1 K at 4.2 K. Sample annealing was done in the Instec probe station under (8.5 ± 0.5)·10−3 mbar, for 60 minutes and at 520 K. The annealing parameters were set to remove the leftovers from the organic mask successfully, and not to damage the ribbons or to compromise the gate dielectric. The resistance values are presented as the total device resistance and were measured in a two-terminal configuration. The presented resistance values are neither corrected for the lead resistance nor for the contact resistance between the nanoribbons and graphite electrodes.

Electrical characterization under varied vacuum conditions (Supplementary Fig. S12) was performed by dosing air with a controlled relative humidity level of 40% using a mass flow meter from Omega, with a flow precision of 0.1 sccm (standard cubic centimeters) and a range of 0.1–200 sccm. Vacuum levels were maintained to below 1% variation during the electrical sweeps. To achieve lower pressure levels, the flow of the humid air was set against the continuous pumping. To achieve the pressure levels in the range of 102 mbar, the chamber was disconnected from the pumps, and a small dose of humid air was released to the chamber.

AFM measurements
AFM measurements were performed using the Horiba/AIST-NTO megascope AFM system. Nanosensor probes were used (resonant frequency ~300 kHz, tip radius of 2–4 nm). AFM images were processed in the open-source software Gwyddion v2.56. For AFM image processing, zero-order line filtering was applied, leveling of the base plane, and tip deconvolution were considered, using a model tip with a radius of 4 nm.

MD simulations
The local configuration of H2O molecules and the binding states at the oxygenated graphene nanoribbon (Gr NR) edges were investigated using a ReaxFF atomistic model incorporated with molecular dynamics (MD) simulations51–53. The model was applied to mon- and bi-layer Gr nanoribbons and includes the elements of the experimental system. The graphene NR is 20 nm wide in the x-direction, terminated with oxygen, and 34 nm long in the y-direction, with applied periodic boundary conditions. To achieve complete (1:1) edge coverage, we chose 4.7 and 9.4 water molecules per nm length of the edge. To test the influence of one edge cluster on the other, the NR width was varied between 5 nm and 30 nm, considering only a mono-layer Gr nanoribbon. We calculated both in-plane and out-of-plane electric fields generated by the edges for various widths of the graphene band. Based on these calculations, we found that for ribbons wider than 5 nm, the electric fields generated within the cluster become two orders of magnitude larger than the average effective field exerted by one cluster on the other.

The MD simulations were performed using time steps of 1 fs and an NVT thermostat in LAMMPS, a commonly used distributed classical MD code51. The simulated water temperature ranged from 250 K to 450 K, depending on the simulation, while the positions of carbon atoms in the ribbon were fixed to maintain a planar structure, similar to the experimental setup. The interatomic forces within the graphene-water system were derived using the appropriate ReaxFF potential by Zhang et al.52,53. In addition, a ReaxFF model by Chenoweth et al.54 was tested, with no qualitative differences observed. Charge equilibration (Qeq) is used in molecular dynamics simulations to calculate the distribution of partial charges on the graphene ribbon edge and water molecules. This distribution changes with time to match changes in the local environment55.

The homogeneous electric field was varied cyclically between −2 V nm−1 and 2 V nm−1 in the simulations used to evaluate the temperature dependence of the hysteresis width, as these field strengths were estimated to be achievable in the experiment, considering the fringing capacitance effects (see also Supplementary Fig. S16). To achieve the saturation, the field was varied cyclically between −4.1 V nm−1 and 4.1 V nm−1. In all cases, the averages over five cycles are shown.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "The manuscript by Aslam et al. reports on the role played by water on the ferroelectric effect of a graphene nanoribbons\ndevice. I think the manuscript is well-written and clearly tackles the water molecule dynamics, although some points should\nbe investigated in more detail. Hence, I recommend revision before publication addressing the following comments.\n- I do not understand the Authors’ choice to use organic molecules like parahexaphenyl (6P) as masks for their nanoribbons\nin both this manuscript and their previous work (ref. 11), when more precise self-assembled nanopatterns can be exploited\n(see for instance Kim et al., Adv. Mater. 2023, 35, 2207338). Is there a specific reason? Do not the Authors believe that a\n“more ordered” pattern would provide more reliable and comparable devices?\n- As per previous entry, a more controlled nanoribbons fabrication would allow to control their width (improving the statistics\nreported in Figure S2) and the pitch, i.e. the distance between two adjacent nanoribbons. I believe that such a higher control\nwill, in turn, enable to answer questions like: can the width of the nanoribbons assumed to be an additional factor among\nthose considered by the Authors (number of graphene NR layers, temperature…)? Is there a trade-off between the width of\nthe nanoribbon and “bulk-like” effects for which the adsorption of water is negligible?\n- The main flaw of the manuscript, in this Reviewer’s opinion, is that no description of the ambient conditions is provided. For\ninstance, caption of Figure 1 states that the transfer characteristics for 2L GNR FETs are measured under low vacuum\n(2x10^-2 mbar), but it is not reported how the water adsorption has occurred in ambient conditions (humidity?). The only\nhumidity value I found was in the caption of Figure 2. However, due to the importance of this factor in the discussion, I think\nthe Authors should clearly add this information throughout the manuscript.\n- Yet again in Figure 1, it not clear the differences (if any) among the devices. In particular the thickness range for multilayer\ndevices is 2-5 layer, cannot the Authors be more specific on the number of devices with 2, 3, … layers? Moreover in the\nmanuscript devices with L>2 are poorly commented, why?\n- The temperature at which the devices were annealed in Figure 2c,d is missing, and again the ambient humidity upon\nexposure (line 129, page 4) is not reported as well.\n- Can the Authors comment and compare the hysteresis window for 1L GNR FET measured at 4.2 K (Figure 2b) and 120-\nsmaller than that in high vacuum (15 V? at 4.2 K in Figure 2b), so do the Authors confirm the observed trend -shrinkage of\nthe hysteresis window- even in vacuum?\n- I do not understand why clustering of water molecules is ruled out in single-layer nanoribbons.\n- The Authors stated that two conditions mut be met to stabilize the observed ferroelectric effect (page 4, lines 142-145). How\ndo the two conditions are competitive to each other? Provided that the clustering of water molecules results into a\ntemperature-independent ferroelectric behavior, can the Authors comment on the marked differences between mono and\nmulti-layered nanoribbons to get the same macroscopic effect? Did the Authors investigate the size of the water molecules\ncluster on the ferroelectric effect?", "reviewer_2": "In this manuscript, the authors report that graphene nanoribbon (Gr NR) devices with oxygen-terminated edges exhibit stable\nhysteresis and suggest that water molecules at the graphene edges induce ferroelectricity. They also note differences in the\nstability of the hysteresis between mono-layer and multi-layer Gr NR FETs, attributing these differences to collective water\nmolecule dynamics, as simulated using molecular dynamics (MD). I find their work interesting and believe it could have\nsignificant implications for future research and applications. However, I am not entirely convinced that the transport\nmeasurements specifically reflect the role of water molecules at the edges, while safely excluding the influence of other\nfactors.\nBelow are some concerns I would like to raise:\n1. Could the authors provide the mobility values for their mono-layer and multi-layer Gr NR devices encapsulated by hBN?\nGiven the high mobility typically observed in hBN-encapsulated devices, it would be useful to determine whether the\nobserved hysteresis is primarily related to the Gr NR edges. Additionally, could the authors explain why the mobility of the\nSF6 RIE device is significantly lower than that of the O2 RIE device in Figure 1(d)?\n2. The authors mention an energy barrier of 25-40 meV for switching a water molecule. Can authors estimate the size of this\nbarrier from Figure 1(e)? Is the size of barrier consistent with their observations? The hysteresis shown in Figure 1(e)\nremains robust up to 420 K, which also appears to differ from their MD model calculations. The authors attribute this\ndiscrepancy to complex interfaces in the experiment, but further clarification would be helpful.\n3. The authors dismiss the role of charge traps because the hysteresis persists at 4.2 K. However, can the authors rule out\nthe possibility of water-air bubbles at the van der Waals interfaces as a contributing factor?\n4. The authors could consider conducting experiments on wider graphene devices encapsulated by hBN. Could they\ncomment on how the width of the graphene affects the hysteresis characteristics?", "reviewer_3": "This manuscript reports on the role of water molecule dynamics in the ferroelectric response of graphene nanoribbon (Gr\nNR) devices. Through an analysis of the dependence on the number of graphene layers and temperature, it reveals the\nstabilizing effect of the collective behavior of water molecules on the ferroelectric effect. To enhance the quality of the\nmanuscript, minor revisions should be considered. All recommendations are listed below.\n1.The manuscript involves the preparation of multilayer graphene nanoribbons, water molecule adsorption, and temperature-\ndependent measurements. However, it does not clarify how experimental consistency was ensured, particularly regarding\nthe reproducibility of experiments under varying temperature and electric field conditions. The authors should demonstrate\nthe consistency of experimental results across different devices.\n2.The experimental results in the manuscript show good consistency with the molecular dynamics simulation results for\nmultilayer graphene nanoribbons (Gr NRs), but there are certain discrepancies in the case of monolayer graphene\nnanoribbons. For example, the opening of the hysteresis window shown in Fig. 3 (a) is much less clear as compared to the\nexperimental results (Figure S4). The authors explained this by introducing the possible more complex interfaces due to\nhBN encapsulation, please explain in more detail.\n3.The time scale of the molecular dynamics simulations (4 ns) in the manuscript significantly differs from the experimental\ntime scale (100 ms). In what way could this substantial discrepancy in time scales affect the accuracy of the simulation\nresults?\n4.The ferroelectric effect in multilayer graphene nanoribbons (2-5 layers) remains nearly unaffected at temperatures up to\nWhile the authors attribute this phenomenon to the absence of collective water molecule behavior in monolayer graphene\nnanoribbons, they do not explore the mechanistic differences in temperature stability between graphene nanoribbons of\nvarying layer numbers.\n5.It would be valuable if the authors could further clarify the potential benefits of designing water-switch-based devices for in-\nmemory computing and synaptic circuit applications.\n6.There are minor typographical errors in the manuscript, such as \"VH\" on line 87 not being italicized. Please ensure that all\ncontent adheres to standard formatting conventions."}, "author_response": "Point-by-point response\nto the reviewers’ comments and questions\n“Ferroelectricity in Graphene Nanoribbon Devices Enabled by Collective Water Molecule Dynamics”\n(NCOMMS-23-18909A-Z)\nThe comments and questions from the reviewers have been separated and presented in regular font, while our responses are\nprovided in blue font, followed by a brief description of the changes made in the revised text, and highlighted in the marked-up\nversion of the revised manuscript and revised SI files, using red font.\nAuthors’ remarks: We sincerely thank all three reviewers for their time and effort in assessing our manuscript, particularly for their\npositive and encouraging feedback. We have carefully addressed all their comments and questions, which has allowed us to\nsignificantly improve the quality and level of our manuscript as detailed in this point-by-point response letter.\n--- REVIEWER COMMENTS ---\nReviewer #1 (Remarks to the Author)\nThe manuscript by Aslam et al. reports on the role played by water on the ferroelectric effect of a graphene nanoribbons device. I\nthink the manuscript is well-written and clearly tackles the water molecule dynamics, although some points should be investigated\nin more detail. Hence, I recommend revision before publication addressing the following comments.\nAuthors’ response: We thank Reviewer #1 for the positive evaluation of our work. We agree with the reviewer that some points\nrequired additional investigations, which, thanks to these constructive comments, helped improve the level and clarity of this work.\nThese revisions include new data and analysis as described below.\nR1Q1: - I do not understand the Authors’ choice to use organic molecules like parahexaphenyl (6P) as masks for their nanoribbons\nin both this manuscript and their previous work (ref. 11), when more precise self-assembled nanopatterns can be exploited (see for\ninstance Kim et al., Adv. Mater. 2023, 35, 2207338). Is there a specific reason? Do not the Authors believe that a “more ordered”\npattern would provide more reliable and comparable devices?\nAuthors’ response: We thank the referee for discussing the BCP nanopatterning process for 2D materials. While other nanoribbon\npatterning methods like BCP provide a higher degree of nanostructure order for device optimization, in addition to our experience\nwith 6P masks, we strategically chose these molecules for our fundamental study due to the simplicity of the fabrication process\nand formation of a clean interface with minimal defects:\n1. Simplicity: The masks leave a bare 2D material surface, as the molecules do not form a film over the 2D material, but\nrather strongly de-wet the substrate (resulting in ribbon formation). Therefore, this method requires only one etching step\nto both pattern the structures and functionalize their edges. This simplicity avoids potential damage to the basal plane and\nline-edge/width roughness associated with multiple RIE steps.\n2. Clean interface: Parahexaphenyl (6P) interacts weakly with 2D materials. This method minimizes unintended doping\nand charge transfer due to mask residues, producing clean interfaces critical to revealing new fundamental phenomena,\nas demonstrated in this work.\n3. Minimizing interfacial defects: We use hetero-stacks of 2D materials (hBN/graphene/hBN), which are prone to the\nformation of interfacial bubbles, even when silicon nitride membranes are used for transfer. Since epitaxial self-assembly\nof organic molecules is highly sensitive to any perturbations of the substrate’s lattice, the nanostructures “avoid”\nassembling on the bubbles as the commensuration changes due to strain. Consequently, interfacial “bubbles” as a\nsource of potential charge traps are largely avoided in nanoribbon networks. Examples of this are provided within the\nanswer to Question 3 posed by Reviewer 2.\nWe agree that the precise control offered by BCP is absolutely needed to bring the observed effect to applications.\nChanges made: We now discuss in the beginning of the Results section the advances in BCP self-assembly and patterning of 2D\nmaterials-based nanoribbons via this method. We also discuss the advantages of BCP for integrating 2D materials-based\nnanoribbons into future technologies.\nreported in Figure S2) and the pitch, i.e. the distance between two adjacent nanoribbons. I believe that such a higher control will, in\nturn, enable to answer questions like: can the width of the nanoribbons assumed to be an additional factor among those considered\nby the Authors (number of graphene NR layers, temperature…)? Is there a trade-off between the width of the nanoribbon and\n“bulk-like” effects for which the adsorption of water is negligible?\nAuthors’ response: This is an excellent point. A narrower width distribution, in addition to a way to control the pitch between the\nadjacent ribbons, could allow for tuning the hysteresis width and optimize performance with respect to the applied gate bias.\nHowever, note that these parameters should not impact the nature of the observed effect. Our conclusions regarding the layer\nnumber and temperature dependence of the hysteretic effect should remain unchanged, even with a narrower ribbon distribution.\nWe expect that there is an optimal width that maximizes the hysteresis effect. However, performing systematic experiments to\nelucidate the role of ribbon width on V goes beyond the scope of this study.\nH\nFollowing the reviewer’s advice, we have estimated the width limit and discussed this issue in more detail in the revised manuscript.\nFrom MD simulations of polarized dipoles, both for 1L and 2L ribbons, we have extracted the dipole-induced fields perpendicular to\nthe basal plane of graphene and their dependence along the x direction of the ribbon width. The fields induced by only one polarized\nedge are illustrated to show that, from the MD model, we obtain the theoretically expected 1/x2 dependence (added in Figure S19).\nThere are two scenarios to consider:\n1. As the ribbon width reduces, at some point, the field induced by one row of the polarized edge dipoles would start\nto affect the dipoles at the other edge. The dipoles from the two edges of a ribbon would effectively try to flip each\nother, as the edge dipole-induced field has the opposite sign of the external field by the gate that introduces their\npolarization. From the MD simulations, we estimate that below a 10 nm width, the field induced by one dipole line\n(edge) starts to affect the polarization of the opposing edge. Below 5 nm, the dipole field would be sufficient to\novercome the external field and would likely result in an antiparallel state. Consequently, we expect a reduction\nin the hysteresis effect for widths below 10 nm.\n2. Looking into the other extreme, as the ribbon width increases, more “bulk” properties of graphene would be\nobserved in the electrical response, and the hysteresis would reduce. Again, from the MD simulations, we estimate\nthat the widths above 100 nm GrNRs would give instead a negligible contribution from the edges to the\nelectrical response. In general, the wider the ribbons, the lower the impact of the edge-induced fields.\nJudging from the dipole-induced field distribution from our MD simulations, we expect that the nanoribbon width in the range of\n10 nm is optimal for maximizing the hysteresis window with respect to the total range of the field that can be applied by the\ngate. Additionally, optimizations with high-k dielectrics, reduced gate thickness, and the use of dual-gate geometries would enable\nfurther enhancement of the water switch's effectiveness.\nChanges made: The impact of the width of the ribbons is discussed in the main text, and is referenced to the added section in the\nSI Figure S7 and Figure S19 where, respectively, the statistics of V for all the devices (see also the reply to question 6) and\nH\nadditional MD results on the field decay along the ribbon width direction are provided. Estimates of the expected optimal ribbon\nwidth, along with other potential improvements of the device architecture, are also discussed in the main text.\n(E = 2 V/nm in the positive y-direction, 300 K) considering only one ribbon edge, in the case of 1L ribbon (a-c) and 2L ribbon (d-\next\nf). (a,d) Contour maps presenting the surface charge density. Electric field streamlines are superimposed to indicate the direction\nof electrostatic influence. The mirror-image placement of charges ensures lateral symmetry, and the color-coded regions emphasize\nstrong confinement effects in the near-surface zone. (b,e) corresponding E field strength map (semi-log scale). (c,f) Intensity of the\ny\ndipole-induced fields (solid blue lines) along the x-direction at 0.2 nm above the graphene plane, in log-log scale. A reference x-2\ntrend line (dashed orange) illustrates the asymptotic decay expected from a 1D line of dipoles, confirming consistency with analytical\npredictions.\nR1Q3: The main flaw of the manuscript, in this Reviewer’s opinion, is that no description of the ambient conditions is provided. For\ninstance, caption of Figure 1 states that the transfer characteristics for 2L GNR FETs are measured under low vacuum\n(2×10-2 mbar), but it is not reported how the water adsorption has occurred in ambient conditions (humidity?). The only\nhumidity value I found was in the caption of Figure 2. However, due to the importance of this factor in the discussion, I think\nthe Authors should clearly add this information throughout the manuscript.\nAuthors’ response: We agree with Reviewer 1, this is indeed an important omission in the description of the fabrication process\nthat has now been corrected.\nIn brief, our RIE system (Oxford’s PlasmaLab 100) is exposed to ambient air, with a regulated temperature and relative humidity\n(298 K and 22 - 26 % RH). A waiting time of at least 30 minutes was applied in these ambient conditions and humidity for\neach device between the RIE step and the electrical measurements under vacuum. Further, no significant difference was observed\nbetween the freshly etched devices (with the 30 minutes holding time) and the devices that were tested even several days after the\nRIE step, also stored under the ambient conditions.\nChanges made: This is now directly addressed in the main text, and more details are provided in the SI methods.\nR1Q4: Yet again in Figure 1, its not clear the differences (if any) among the devices. In particular the thickness range for multilayer\ndevices is 2-5 layer, cannot the Authors be more specific on the number of devices with 2, 3, … layers? Moreover, in the manuscript\ndevices with L>2 are poorly commented, why?\nAuthors’ response: We are thankful to Reviewer 1 for this question, which allowed us to clarify our idea for multi-layer devices.\nWe tested thicker devices (up to 5 L) and have not observed any systematic difference between the response of the L > 2 FETs.\nThe limit of 5 L was chosen as an empirical estimate to prevent longer RIE steps from damaging the ribbons. During our\ninvestigation, we observed that all monolayer devices exhibited a strong temperature dependence of the V , whereas all multilayer\nH\ndevices displayed consistent temperature-independent hysteresis. Therefore, we have grouped the devices to emphasize the\nfundamental differences between monolayer and multilayer systems. The goal of Figure 1e,f is to present the evidence of this\nclaim.\nrevised SI, we specify the exact number of devices used in our study and include statistical graphs comparing the apparent peak\nmobility and V values across all devices (see also the reply to question 6 and reviewer 2 question 1).\nH\nChanges made: We have edited the main text and Figures 1 and 4 to clearly indicate the number of layers. We also ensure that\nwithin the SI the layer number is clearly communicated for all presented results.\nR1Q5: The temperature at which the devices were annealed in Figure 2c,d is missing, and again the ambient humidity upon\nexposure (line 129, page 4) is not reported as well.\nAuthors’ response: We thank the reviewer for noticing this omission. The annealing was performed at 250°C for 60 minutes\nunder 8 × 10-3 mbar vacuum. Re-exposure to ambient relative humidity (22-26 %) is now also specified.\nChanges made: The humidity and the annealing parameters are now clearly stated in the main text, in the SI Methods, the caption\nof Fig. 2, and in the SI Fig. S11.\nR1Q6: Can the Authors comment and compare the hysteresis window for 1L GNR FET measured at 4.2 K (Figure 2b) and 120-400 K\n(Figure 1e)? Moreover, it seems that the hysteresis window under low vacuum (11.2 V at 300 K? in Figure S8) is smaller than that\nin high vacuum (15 V? at 4.2 K in Figure 2b), so do the Authors confirm the observed trend -shrinkage of the hysteresis window-\neven in vacuum?\nAuthors’ response: We thank the reviewer for this insightful point. Our results show that qualitatively, at 4.2 K, the device shows\na robust hysteresis consistent with the trend in Figure 1e, taken as the low-temperature limit where the thermal energy is insignificant\nto cause depolarization. A direct quantitative comparison between these two experiments, shown in Figure 2b and the SI Figure S12\n(originally Fig. S8), is not feasible since the data were obtained from two different devices in two separate measurement systems.\nThe need for different devices arises from the varying environmental conditions required for each experiment. Further, V variations\nH\nare experienced between devices, primarily due to differences in ribbon width. Therefore, we have refrained from making direct\ncomparisons of the absolute values between the different devices.\nTo demonstrate this and provide a clear overview of the VH scattering, we have added a new statistical summary of the hysteresis\nwindow for the entire set of 22 devices used in this study (16 O etched and 6 SF etched). These additions are now discussed in\n2 6\nthe revised manuscript.\nThe data presents the hysteresis window over the total gate sweeping range (V /V ) as a function of the number of layers\nH SG-tot\n(Fig. S7a), V /V as a function of V (Fig. S7b), and also unscaled V as a function of V (Fig. S7c).\nH SG-tot SG-tot H SG-tot\nTo answer the question about the vacuum trend, we have carried out a new experiment where we investigate the dependence of\nthe V on a single device from the vacuum level (8×10-4 mbar) to ambient conditions (40% RH). These new results, presented in\nH\nFigure S12, show that the V reduces by about 40% when going from the ambient pressure down to about 101 mbar, and we did\nH\nnot observe any further reduction of the V as the pressure was decreased further down to a high vacuum. This definitely confirms\nH\nthat the hysteresis window shrinks when the ribbons are exposed to the vacuum conditions and that it remains stable\nbelow 101 mbar. The apparent trend noted by the reviewer is likely due to device-to-device variability, not a continuous\npressure-dependent effect. This question is now discussed in the revised manuscript text.\nChanges made: Along with the discussion on the potential ribbon width effects on the V , in the revised manuscript we also refer\nH\nto the V variations within our device set (added SI Figure S7), we discuss that the scattering likely originates from the variations in\nH\nthe ribbon widths, and also argue why only relative V values are compared between different devices.\nH\nFigure S12 (original version SI Figure S8) has been replaced with more detailed experiments that show how the vacuum level\nimpacts the V in a range from 8×10-4 mbar up to the ambient.\nH\nH 2 6\nFETs). Red circles represent O etched devices and purple diamonds represent SF etched ones. Shaded areas serve as a guide\n2 6\nto the eye to indicate the trends and the device-to-device scattering of the parameters. The data are extracted from the electrical\ntransfer curves measured at 300 K and under 2·10-2 mbar. Each data point with its error bar is extracted as a mean and standard\ndeviation from five transfer curves. (a,b) V is scaled by the total V sweeping range (V ), as a function of the number of\nH SG SG-tot\ngraphene layers (a) and as a function of V (b). (c) presents V (V ) without scaling. Likely, the main source for the\nSG-tot H SG-tot\ndevice-to-device scatter is the variation in the width of the nanoribbons in the self-assembled networks.\nFigure S12: (a) electrical transfer curves (total device resistance; two-terminal) of a 1L Gr NR FET measured at 300 K in ambient\npressure and 40 % relative humidity (RH), and immediately after under a vacuum of 8.3·10-4 mbar. (b) for the same device, the\ndependence of the V on the vacuum conditions of the electrical measurement chamber. For each data point, V was extracted\nH H\nfrom the electrical transfer curves measured at 300 K with V sweeping between -40 V and +40 V at a 2.5 Vs-1 rate. Vacuum level\nSG\nwas controlled by adjusting the inlet of 40% RH air using a 0.1-200 sccm flow controller. The measurements were conducted in two\nways: starting from the ambient and lowering the pressure, and starting from the lowest pressure and increasing to the ambient.\nError bars indicate a standard deviation of the V from three subsequent sweeps. Going from 10-3 mbar to the ambient pressure\nH\nand RH 40% the V opening increases by about 40 %, and the noticeable increase in V starts only after 101 mbar pressure.\nH H\nAuthors’ response: We thank the reviewer for pointing out that this issue was not clearly discussed in the original manuscript. We\ndid not intend to rule out clustering in mono-layers, but rather to highlight the fundamental differences in the nature of the dynamics\nthat give rise to clustering. To clarify, although less frequently than in bi-layers, water clusters do form in mono-layers; however,\nthe complex water dynamics involve motion on the graphene basal plane and flips of individual molecules, which make water\nclusters less stable. In contrast, in bi-layers, the motion of individual molecules across the basal plane is less frequent, while\nindividual molecule flipping was not observed, resulting in larger and more stable clusters. The reasons for these differences are\nfurther discussed in our response to the following comment below (R1Q8).\nChanges made: This point is clarified in the revised manuscript, now stating that the formation of clusters and the collective motion\nof water molecules are not excluded for the mono-layers. Additional data on the cluster size statistics per unit length of the edge\nare provided in the SI Figure S15.\nFigure S15. Probability of water cluster density for (a) mono-layer ribbon, and (b) bi-layer ribbon. The data is extracted from a\nrelaxed configuration at 300 K and without applying the external field. The results demonstrate two points: first, the formation of\nwater clusters cannot be ruled out for the mono-layer ribbons, and second, individual water molecules adsorbed at the ribbon edges\nwere essentially not found in the case of the bi-layers. Please, mind that the bars below 4 nm-1 correspond to the individually\nadsorbed molecules, and that for the bi-layers only the mobile water molecules are counted, not the “bridge” molecules between\nthe layers.\ndo the two conditions are competitive to each other? Provided that the clustering of water molecules results into a\ntemperature-independent ferroelectric behavior, can the Authors comment on the marked differences between mono and\nmulti-layered nanoribbons to get the same macroscopic effect? Did the Authors investigate the size of the water molecules\ncluster on the ferroelectric effect?\nAuthors’ response: These two conditions are not directly competitive but regard two characteristics of the ferroelectric effect. In a\nsimplified picture, we can think of both mono- and bi-layer cases as being a bi-stable system with two states and a barrier between\nthem. The barrier to cross between the two sides of the ribbon (the two energy states) is rather similar for both cases, that is the 1st\ncondition – low thermal barrier. This condition allows the water dipoles to be switched by the gate at the device level. The 2nd\ncondition – collective self-stabilizing cluster – could be seen as a deepening of one minimum, depending on how many\nmolecules are in that minimum state. This introduces a discrepancy in the rate of returning to the lower occupied minima by thermal\nfluctuations (see added Fig. S13), and explains how the dipoles remain in a polarized state after switching, resulting in a memory\neffect.\nTherefore, the existence of water interaction and the size of the cluster play a key role in achieving the temperature-independent\nhysteretic effect. There is, of course, a limit on how large the cluster can be and remain anchored at the edge. The larger the cluster,\nthe larger is the dipole induced field (results added to Fig. S14a), and at the same time, the lower the energy needed to remove the\nindividual molecules at the cluster’s periphery since they are less strongly bound (results added to Fig. S14b). These two effects\nare competitive with each other.\nRegarding the differences between mono- and multi-layers, the critical point is how each system supports water clusters. In\nmulti-layers, the water clusters are stabilized by strong interlayer bridging water molecules acting as a stabilizing scaffold. This\ninterlayer anchoring effect enables the formation of larger and more robust water clusters, resulting in the collective\ntemperature-independent behavior observed in our experiments. The interlayer anchoring effect is absent in monolayers (there are\nno adjacent layers to build the stabilizing scaffold), which limits the size and stability of the water clusters that can be formed.\nTherefore, in contrast to multi-layers, monolayers are dominated by the dynamics of individual water molecules, resulting in\ntemperature-dependent effects.\nOur MD calculations presented in Figure 3 of the manuscript used the largest fully stable clusters with respect to the entire\ntemperature and the external field ranges (up to 4.1 V/nm). Under these conditions, the largest stable clusters at the ribbon edge\ncorrespond to 15 water molecules/nm for mono-layers and 20 molecules/nm for bi-layers. Under lower fields and lower\ntemperatures, larger clusters are expected.\nChanges made: This issue is now clearly addressed in the revised manuscript. Two figures have been added to the SI. Fig. S13\nillustrates the impact of the two conditions on the switching process and the barrier height. Fig. S14 provides additional calculation\ndetails on the relation between the amount of water per unit length of the ribbon edge and the stability of water molecules within the\ncluster at high external fields.\np\nsystem in equilibrium (equal distribution of the particles between the two states), with the exchange rate between the two states\ndescribed by k(T), and the energy barrier between the two states (E ). (b) If the external field (E ) is applied, the potential profile is\na ext\nperturbed and one minimum is favored. (c) Upon canceling the E , the capability of the system to stay polarized depends on the\next\nratio between the barrier height and the thermal energy of the particles. In (a-c), we do not consider interaction between the particles,\nand even if one state is more occupied, the rate of crossing between the states is the same. In contrast, (d-f) considers that the\ninteraction between the particles favors the state with higher occupation, i.e., reduces the exchange rate. This would be the case\nwhen collective dynamics self-stabilizes the system in a polarized state. Let’s assume that by thermal fluctuations one particle\ncrosses over the barrier (d) and introduces non-equal occupation of the states. The interaction between the particles and the\ndiscrepancy in the occupation between the two states introduces the energy difference between the two minima E and more\npopulated state is slightly more favored (e). Lastly, if we would polarize such an interacting system in an external field (as illustrated\nin (b)) and switch the external field off, being significantly more populated, one state would be significantly more favored (f). In this\ncase, the exchange rates between the two states would significantly differ, and to de-polarize the system, a much larger effective\nbarrier would be observed.\nwater molecules present in the simulation and expressed per unit length of the edge. The model is relaxed in the external field of\nsimulation in a high external field of 4 V/nm, shown as a function of the number of water molecules per unit length of the edge.\nCircles and crosses represent the calculation results for 1L and 2L ribbons, respectively. For both sets of calculations, the\ntemperature was set to 300 K.\nReviewer #2 comments: In this manuscript, the authors report that graphene nanoribbon (Gr NR) devices with oxygen-terminated\nedges exhibit stable hysteresis and suggest that water molecules at the graphene edges induce ferroelectricity. They also note\ndifferences in the stability of the hysteresis between mono-layer and multi-layer Gr NR FETs, attributing these differences to\ncollective water molecule dynamics, as simulated using molecular dynamics (MD). I find their work interesting and believe it could\nhave significant implications for future research and applications. However, I am not entirely convinced that the transport\nmeasurements specifically reflect the role of water molecules at the edges, while safely excluding the influence of other factors.\nAuthors’ response: We are thankful to Reviewer 2’s constructive comments; we agree that excluding the influence of other factors\nis critical to this work. To address this crucial point, we highlight four key experimental observations that collectively demonstrate\nthe role of edge-adsorbed water in the ferroelectric behavior of Gr NR FETs.\n1. The hysteresis effect appears only in devices with hydrophilic, oxygen-terminated edges. Identical devices with\nhydrophobic edges obtained with SF RIE, with all the other parameters kept the same (fabrication steps, graphene\n6\nthickness, hBN capping, stacking methods, graphite contacts, exposure to ambient humidity), show negligible hysteresis.\nThis points to an effect related to the edge’s ability to adsorb water.\n2. Vacuum annealing at 520 K completely removes the hysteresis behavior in oxygen-terminated devices. However, when\nthe same device is re-exposed to ambient humidity and measured again under vacuum, we see the reappearance of the\nhysteretic behavior, proving the direct role of water adsorption from the air.\n3. Potential contributions from trapped charges due to the formation of interfacial bubbles are avoided in our fabrication\nmethod with self-assembled masks. This is confirmed by new AFM results in Figure S3.\n4. Other contributions due to environmental exposure of the nanoribbon basal planes are minimized by hBN encapsulation\nof the top and bottom sides of the graphene channel. This device configuration isolates the device from environmental\ncharge traps, ensuring that the edges and not the basal planes are selectively exposed for interaction with the environment.\nThese independent observations provide strong evidence to support our conclusion that the ferroelectric effect observed originates\nfrom the dynamics of water molecules at the Gr NR edges instead of other potential artifacts.\nR2: Below are some concerns I would like to raise:\nR2Q1: Could the authors provide the mobility values for their mono-layer and multi-layer Gr NR devices encapsulated by hBN?\nGiven the high mobility typically observed in hBN-encapsulated devices, it would be useful to determine whether the\nobserved hysteresis is primarily related to the Gr NR edges. Additionally, could the authors explain why the mobility of the\nSF6 RIE device is significantly lower than that of the O2 RIE device in Figure 1(d)?\nAuthors’ response: Following the reviewer’s advice, we performed a comprehensive mobility analysis for all 22 devices used in\nthis study. The data are now included in the revised SI. The following results have been added:\n Statistics of the apparent electron mobility (not excluding potential contributions of the contact resistance) at 300 K and\nfor the total of 22 devices used in the study (6 SF etched and 16 O etched) are now added to the SI Figure S5. The\n6 2\ngraph provides the peak mobility values from µ(V ) calculated using the parallel capacitance model. The model is chosen\nG\nto provide a straightforward comparison between starting flake FETs and Gr NR FETs. If the fringing capacitance model\nis used (1D wire over a plane), about 1 order of magnitude lower mobility values are obtained for the Gr NR FETs. More\ndetails on the fringing capacitance model are available in [10.1103/PhysRevLett.106.256801] [10.1038/s41699-022-\n00356-y].\n Several representative µ(V ) graphs have been added in the SI Figure S6, two O-terminated and two F-terminated\nG\ndevices, presenting the data for 1L and 2L Gr NRs in each case. µ(V ) curves are provided before and after the RIE step\nG\nfor the same device.\nRegarding the SF RIE devices and mobility, we conducted a negative experiment with six SF RIE devices. Since the comparison\n6 6\nbetween O and SF RIE was previously reported for the edges of graphene flakes [10.1021/acs.nanolett.8b00797], a larger data-set\n2 6\non SF RIE devices was unnecessary.\nto O2-etched ones. Since these devices have been used only as a negative experiment, we preselected the lower-performing\ndevices for the SF RIE process. These results confirm that the main difference between the two device types is the presence or\n6\nabsence of water-induced hysteresis, O vs. SF -etched, and not their carrier mobility.\n2 6\nRegarding the O RIE devices and mobility: Considering the comparisons between flake- and NR-FETs, and using a parallel\n2\ncapacitance model, µ(V ) values were found to be in a similar range as for the original flakes. An increase in peak µ(V ) values is\nG G\ncommonly observed for O RIE-FETs. Further, µ(V ) curves were found frequently to be asymmetric for electrons and holes around\nCNP+/CNP-, and also exhibit a sweep direction dependency. We attribute this to the interplay between the gate and the induced\ndipole fields, which can enhance the apparent mobility on one side of the CNP depending on the sweep direction, correlating the\nelectrical characteristics to the dipole switching.\nChanges made: Mobility values for the entire device set used in this study and gate-dependent mobility curves are added to the SI\nFigure S5 and Figure S6. These new results and the difference between the parallel plate and fringing capacitance models are\ndiscussed in the revised text.\nlin\nvalues are extracted as peak mobilities from the electrical transfer curves recorded at 300 K (2·10-2 mbar). The total of 22 devices\nthat were used in the study are presented. Hexagons correspond to the flake-FETs before the RIE step (step 5 in the fabrication\nscheme – Figure S1), diamonds represent Gr NR-FETs etched using SF RIE, and circles using O RIE. Mobilities are extracted in\n6 2\nall cases using a linear regime approximation and a parallel capacitance model. Fringing capacitance model for the NR-FETs results\nin about one order of magnitude lower values [S2]. Length and width of each device were determined via a combination of the\noptical micrographs and AFM topography images. The larger uncertainty for the Gr NR-FETs originates from the uncertainty of the\nNR widths measured by the AFM.\nFigure S6: Linear apparent electron mobility as a function of gate voltage. Each sub-panel compares the µ (V ) curves of devices\nlin SG\nbefore and after the RIE step (labeled as flake- and NR-FETs, respectively). The x-axis is compensated for the V values for\nCNP\nbetter comparison. The curves are obtained from the electrical transfer data recorded at 300 K (2·10-2 mbar) using the linear regime\nparallel capacitance mobility model. (a) and (b) presents two devices with O RIE etching of 1 and 2 graphene layers, respectively.\n2\n(c) and (d) present comparable devices to (a) and (b), only with SF RIE step.\n6\nbarrier from Figure 1(e)? Is the size of barrier consistent with their observations? The hysteresis shown in Figure 1(e) remains robust\nup to 420 K, which also appears to differ from their MD model calculations. The authors attribute this discrepancy to complex\ninterfaces in the experiment, but further clarification would be helpful.\nAuthors’ response: We thank the reviewer for this multi-part question that led us to perform a new analysis that significantly\nstrengthened our conclusions. Following the reviewer’s advice, we estimated the activation energy (E) by fitting our experimental\na\ndata to an Arrhenius model: . Here, k is the Boltzmann constant, T is the temperature, and V simply\nB H0\nscales the fit to the value obtained at the lowest temperature (considering the experimental range). From the V (T) experiments\nH\npresented in Figure 1e of the𝑉H o⁄r𝑉ig H i 0 na=l 1m−an𝑒u𝑥s𝑝c(r−ip𝐸t, a f⁄o𝑘r B t𝑇h)e mono-layers, we have obtained Ea = (17.45 ± 11.60) meV. The data for\nthe multi-layers (Figure 1h) yields a barrier higher than 150 meV.\nConcerning the discrepancy with the MD model and clarification of “complex interfaces”, the reviewer correctly noticed that the\nexperimentally deduced hysteresis values seem to be more robust than the obtained ones in our idealized MD model, showing a\nweaker hysteresis. We attribute this discrepancy to several factors present in the experiments (our “complex interfaces” case) that\nare not introduced in the MD simulations. First, the hBN encapsulation in our experiments can provide an additional electrostatic\npotential that contributes to the water cluster stabilization. In contrast, this effect is not considered in the simulations that are based\non an isolated Gr NR. Second, in comparison to the nanosecond time scale in which the simulations were run, the experimental\nvalues were averaged over a much longer time scale (>100 ms), allowing the water clusters more time to relax and to settle into a\nstable polarized state, resulting in a more pronounced hysteresis. Following the reviewer’s advice, this point is now discussed in the\nrevised text.\nChanges made: The estimation of the energy barrier is briefly discussed in the text, the fits and the E values are given in\na\nFigure 1e,f. Further, Figure S8 is added to the SI, showing separately the fits for each of the mono-layer devices.\nFigure 1: …. (e,f) Temperature dependence of the hysteresis windows (V ) normalized to the value recorded at 120 K, presenting\nH\nthe data for the mono-layer and multi-layer (2 to 5 layers, with thickness labelled at the right side for each curve in (f)) Gr NR FETs.\nDifferent lines represent different devices, and error bars depict V variation in the subsequent measurements.\nH\nShaded areas serve as a guide to the eye. Solid black lines are the Arrhenius function-based analytical fits used to extract the energy barrier (E).\na\nDashed black lines in (e) correspond to one sigma deviations of the E estimate for mono-layers.\na\nFigure 1e of the main manuscript. The fitting is carried out by the root mean square minimization, and the uncertainties imply a\nstandard deviation.\nR2Q3: The authors dismiss the role of charge traps because the hysteresis persists at 4.2 K. However, can the authors rule out the\npossibility of water-air bubbles at the van der Waals interfaces as a contributing factor?\nAuthors’ response: We thank the reviewer for noticing this critical point that was not carefully addressed in the original manuscript.\nIndeed, we can rule out contributions from interfacial air bubbles, which is one of the benefits of our fabrication process.\nThe organic nanostructure self-assembly masks obtained by hot-wall epitaxy (HWE) growth near thermodynamic equilibrium are\nhighly sensitive to substrate strain. Interfacial bubbles introduce strain in the top hBN layer, which hinders this growth. Consequently,\nthe organic nanostructure masks, and therefore the final graphene nanoribbon network, preferentially form over pristine unstrained\nregions, avoiding the bubbles. Hence, the Gr NR network is formed away from these potential charge traps.\nChanges made: This issue is now clearly discussed in the main text, explaining how our method allows minimizing the influence\nfrom interfacial air bubbles. The SI has been updated with a new Figure S3, providing direct evidence by AFM topography imaging\nof the nanostructure’s assembly, avoiding the bubbles.\nFigure S3. … (d-g) AFM topography images of the nanoribbon networks before the RIE step (scale bar 500 nm, z scale 30 nm).\nObserved bubbles are characteristic interfacial traps of water-air mixture within the layers of the van der Waals hetero-stack\n(hBN/graphene/hBN). Bi-axial strain of the top layer induced by the bubbles hinders the self-assembly of the organic nanostructures.\nConsequently, the organic nanostructures tend to form between the bubbles or tend to terminate on the bubbles.\non how the width of the graphene affects the hysteresis characteristics?\nAuthors’ response: We thank the reviewer for bringing up this excellent point. In flakes (width commonly above 2 µm), this effect\nis not observed, unless the current is passed along the flake’s edge [10.1021/acs.nanolett.8b00797]. Otherwise, an alternative\napproach to investigate the impact of the ribbon width on the V would be to use EVU lithography or self-assembly of block\nH\ncopolymers. Although such a systematic experimental investigation is beyond the scope of this work, we performed molecular\ndynamics simulations to elucidate how NR’s width affects hysteresis.\nChanges in nanoribbon width should affect the magnitude of the hysteresis window, but neither the nature of the observed effect\nnor the conclusions regarding the layer number and temperature dependence.\nReviewer 1 raised similar concerns in their questions 1 and 2. Therefore, the following part of the answer is repeated from\nReviewer 1 Question 2:\nWe have estimated the width limit and discussed this issue in more detail in the revised manuscript. From MD simulations of\npolarized dipoles, both for 1L and 2L ribbons, we have extracted the dipole-induced fields perpendicular to the basal plane of\ngraphene and their dependence along the x direction of the ribbon width. The fields induced by only one polarized edge are\nillustrated to show that, from the MD model, we obtain the theoretically expected 1/x2 dependence (added in Figure S19).\nThere are two scenarios to consider:\n As the ribbon width reduces, at some point, the field induced by one row of the polarized edge dipoles would start\nto affect the dipoles at the other edge. The dipoles from the two edges of a ribbon would effectively try to flip each\nother, as the edge dipole-induced field has the opposite sign of the external field by the gate that introduces their\npolarization. From the MD simulations, we estimate that below a 10 nm width, the field induced by one dipole line\n(edge) starts to affect the polarization of the opposing edge. Below 5 nm, the dipole field would be sufficient to\novercome the external field and would likely result in an antiparallel state. Consequently, we expect a reduction\nin the hysteresis effect for widths below 10 nm.\n Looking into the other extreme, as the ribbon width increases, more “bulk” properties of graphene would be\nobserved in the electrical response, and the hysteresis would reduce. Again, from the MD simulations, we estimate\nthat the widths above 100 nm GrNRs would give rather a negligible contribution from the edges to the\nelectrical response. In general, the wider the ribbons, the lower the impact of the edge-induced fields.\nJudging from the dipole-induced field distribution from our MD simulations, we expect that the nanoribbon width in the range of\n10 nm is optimal for maximizing the hysteresis window with respect to the total range of the field that can be applied by the\ngate. Additionally, optimizations with high-k dielectrics, reduced gate thickness, and the use of dual-gate geometries would enable\nfurther enhancement of the water switch's effectiveness.\nChanges made: The impact of the width of the ribbons is discussed in the main text, and is referenced to the added section in the\nSI Figure S7 and Figure S19 where, respectively, the statistics of V for all the devices (see also the reply to reviewer 1 question 6)\nH\nand additional MD results on the field decay along the ribbon width direction are provided. Estimates of the expected optimal ribbon\nwidth, along with other potential improvements of the device architecture, are also discussed in the main text.\n(E = 2 V/nm in the positive y-direction, 300 K) considering only one ribbon edge, in the case of 1L ribbon (a-c) and 2L ribbon (d-f).\next\n(a,d) Contour maps presenting the surface charge density. Electric field streamlines are superimposed to indicate the direction of\nelectrostatic influence. The mirror-image placement of charges ensures lateral symmetry, and the color-coded regions emphasize\nstrong confinement effects in the near-surface zone. (b,e) corresponding E field strength map (semi-log scale). (c,f) Intensity of the\ny\ndipole induced fields (solid blue lines) along the x-direction at 0.2 nm above the graphene plane, in log-log scale. A reference x-2\ntrend line (dashed orange) illustrates the asymptotic decay expected from a 1D line of dipoles, confirming consistency with analytical\npredictions.\nThis manuscript reports on the role of water molecule dynamics in the ferroelectric response of graphene nanoribbon (Gr NR)\ndevices. Through an analysis of the dependence on the number of graphene layers and temperature, it reveals the stabilizing effect\nof the collective behavior of water molecules on the ferroelectric effect. To enhance the quality of the manuscript, minor revisions\nshould be considered. All recommendations are listed below.\nAuthors’ response: We thank the reviewer for the positive evaluation and constructive recommendations. We have addressed all\nthe reviewers’ comments as described below.\nR3Q1: The manuscript involves the preparation of multilayer graphene nanoribbons, water molecule adsorption, and\ntemperature-dependent measurements. However, it does not clarify how experimental consistency was ensured, particularly\nregarding the reproducibility of experiments under varying temperature and electric field conditions. The authors should demonstrate\nthe consistency of experimental results across different devices.\nAuthors’ response: We thank the referee for bringing up this important point. Following the reviewer’s advice, we took care of this\nquestion, addressing it in two ways:\n1) We conducted and included a comprehensive statistical analysis of all 22 Gr NR devices used in this study, comprising\ngraphene thicknesses from 1L to 5L. These results are included in the revised SI, new Figures S5 and S7. Although there\nis some expected scatter in the absolute mobility and hysteresis values due to the random nature of the nanostructure\nself-assembly process, the statistical results show that the trends and performance reported are consistent across all\ndevices.\n2) To clarify how we maintained experimental consistency, we expanded the SI Methods section to include all fabrication\nconditions. These conditions adhered to the standards commonly found in nanofabrication labs, with each device’s\nfabrication parameters closely monitored. All electrical measurements were conducted using systematically defined\nmeasurement routines.\nChanges made: The SI Methods now state additional information on the experimental procedures, especially focusing on ensuring\ncomparable experimental conditions between the different devices.\nThese include:\n- Holding times in the HWE setup before the growth to ensure stable growth temperatures (45 min).\n- Variations of the temperatures for the HWE experiments are also stated.\n- Holding times in the electrical setup to ensure stable temperature during variable temperature electrical measurements;\n10 min for the flow LN2 cryostat and the 77 K to 425 K range, and 60 min for the closed-cycle LHe cryostat and 4.2 K\nmeasurements.\n- Temperature variations are specified for the electrical setups.\n- Vacuum variations are specified for the electrical setups.\n- Device annealing time and temperature are stated (see also reviewer 1 question 5)\n- Relative humidity exposure after RIE and after vacuum annealing is now clearly defined (humidity values and exposure\ntime) - (see also response to reviewer 1 questions 3 and 5)\nSI Figures S5 and S7 are added to provide an overview of the peak apparent electron mobility and a summary of the hysteresis\nwindow (V ). The addition is also briefly discussed in the main text (see also reviewer 2 question 1, and reviewer 1 question 6).\nH\nlin\nThe values are extracted as peak mobilities from the electrical transfer curves recorded at 300 K (2·10-2 mbar). The total of 22 devices\nthat were used in the study is presented. Hexagons correspond to the flake-FETs prior to the RIE step (step 5 in the fabrication\nscheme – Figure S1), diamonds represent Gr NR-FETs etched using SF RIE, and circles using O RIE. Mobilities are extracted in\n6 2\nall cases using a linear regime approximation and a parallel capacitance model. Fringing capacitance model for the NR-FETs results\nin about one order of magnitude lower values [S2]. Length and width of each device were determined via a combination of the\noptical micrographs and AFM topography images. The larger uncertainty for the Gr NR-FETs originates from the uncertainty of the\nNR widths measured by the AFM.\nFigure S7. Statistics of the hysteresis window (V ) considering the entire set of devices used in the study (16 O RIE and 6 SF RIE\nH 2 6\nFETs). Red circles present O etched devices and purple diamonds present SF etched ones. Shaded areas serve as a guide\n2 6\nto the eye to indicate the trends and the device-to-device scattering of the parameters. The data is extracted from the electrical transfer\ncurves measured at 300 K and under 2·10-2 mbar. Each data point with its error bar is extracted as a mean and standard deviation\nfrom five transfer curves. (a,b) V is scaled by the total V sweeping range (V ), as a function of the number of graphene layers\nH SG SG-tot\n(a) and as a function of V (b). (c) presents V (V ) without scaling. Likely, the main source for the device-to-device scatter is\nSG-tot H SG-tot\nthe variation in the width of the nanoribbons in the self-assembled networks.\nmultilayer graphene nanoribbons (Gr NRs), but there are certain discrepancies in the case of monolayer graphene nanoribbons.\nFor example, the opening of the hysteresis window shown in Fig. 3 (a) is much less clear as compared to the experimental\nresults (Figure S4). The authors explained this by introducing the possible more complex interfaces due to hBN encapsulation,\nplease explain in more detail.\nAuthors’ response: We thank the reviewer for this question, and agree that this issue was not clearly explained in the original\nmanuscript. After the reviewer’s comment, we have now provided more details about the discrepancy between mono-layers and\nmulti-layers in the revised manuscript, clarifying that the collective motion of water and cluster formation is not excluded for the\nmono-layers (please see also our response to reviewer 1’s question 7). We attribute these differences to the stabilizing effect of the\nhBN encapsulation, which was not included in the MD simulations due to its computational complexity. In short, the hBN\nencapsulation with the polar B-N bonds provides an additional electrostatic potential near the mono-layer’s edge that favors the\nwater molecules’ collective behavior. Another consequence of the hBN encapsulation is the spatial confinement it creates, hindering\nthe motion of water molecules. Since these effects were not captured in the MD simulations, there is a discrepancy with the actual\nexperimental results.\nChanges made: The expected influence of the hBN caps on the V difference between the experiment and the MD simulation is\nH\nnow discussed in more detail.\nR3Q3: The time scale of the molecular dynamics simulations (4 ns) in the manuscript significantly differs from the experimental time\nscale (100 ms). In what way could this substantial discrepancy in time scales affect the accuracy of the simulation results?\nAuthors: We thank the referee for pointing out this issue that was insufficiently addressed in the original manuscript. The different\ntime scales are a known challenge when comparing MD simulations with actual experiments. In our case, the system has more time\nto relax, meaning that upon the application of the bias, stronger polarization of the dipoles would be expected in our experiments.\nTo compensate for this in the simulations, we applied higher fields than achievable by the experiments to speed up the system\ndynamics and probe the stability of the water clusters and their collective response. From this perspective, we had to limit the lower\ntemperature for the model to 250 K, as at lower T, a much longer time would be needed for the relaxation. Therefore, and considering\nthat the MD model is idealized and a much smaller system than the one probed in the experiment, we relied on these simulations\nto reveal the physical mechanism, but not for a direct quantitative match with the experimental data. The MD simulations proved to\nbe a valuable tool for identifying trends that support our main experimental findings, particularly the temperature dependence in\nmonolayers versus the independence in bilayers, thus confirming the mechanism we proposed.\nChanges made: The discussion on the impact of the differing time scales between the experiments and the MD simulations is now\nexpanded. We clearly state that, due to these discrepancies, the absolute values for the remnant fields and the hysteresis opening\ncannot be directly compared to the experimental results in a quantitative manner. Instead, the temperature trends are effectively\ncaptured and used for the comparison with the experiment, validating the physical mechanism we proposed.\nR3Q4: The ferroelectric effect in multilayer graphene nanoribbons (2-5 layers) remains nearly unaffected at temperatures up to\nauthors attribute this phenomenon to the absence of collective water molecule behavior in monolayer graphene nanoribbons,\nthey do not explore the mechanistic differences in temperature stability between graphene nanoribbons of varying layer numbers.\nAuthors’ response: Following the reviewer’s advice, we have added a new Figure S13 to the revised SI to better illustrate the\ndifference in temperature stability, showing the self-stabilizing collective state in multi-layers. This is a critical mechanistic difference\nthat is based on the formation of water bridges in multi-layers (L > 2), which is physically impossible in mono-layers. These bridges\nact as anchoring sites for water molecules, enabling the formation of larger and more robust water clusters compared to the case\nof mono-layers. Consequently, the strong intermolecular interactions mediated by these water bridges give rise to a self-stabilizing\nstate depicted in Figure S13, which can be seen as a deepening of the potential energy well of the polarized state. This provides\nan additional energy barrier for depolarization in multi-layers, making the system resistant to thermal energy and thereby providing\ntemperature-independent behavior to the hysteresis, as observed in the V (T) trend for all multi-layered devices we investigated.\nH\nIn contrast, in mono-layers that cannot support the formation of these water bridges, the stabilizing anchoring effect is missing,\nmaking the water clusters at the edges less stable and their polarization more easily affected by thermal energy, which explains the\nobserved temperature dependence.\nRegarding the mechanistic differences for a varying number of layers, in Figure 1f of the revised manuscript, we now specify the\nnumber of layers for each curve. The critical observation is that these plots vs. temperature are flat and overlapping, with no\ndiscernible trend, in contrast to the results for mono-layers in Figure 1e, which show temperature-dependent behavior. This is why\nthe mechanistic difference we focus on is between the mono-layer case and the multi-layer L>2 case.\nChanges made: We now discuss in more detail these mechanistic differences in temperature stability. A conceptualization of the\nbi-stable system and the differences between such a system with and without the effects of the inter-particle interactions and\nself-stabilization are mentioned in the main text and are discussed in more detail in the SI Figure S13.\np\nsystem in equilibrium (equal distribution of the particles between the two states), with the exchange rate between the two states\ndescribed by k(T), and the energy barrier between the two states (E ). (b) If the external field (E ) is applied, the potential profile is\na ext\nperturbed and one minimum is favored. (c) Upon canceling of the E the capability of the system to stay polarized depends on the\next\nratio between the barrier height and the thermal energy of the particles. In (a-c) we do not consider interaction between the particles\nand even if one state is more occupied the rate of crossing between the states is the same. In contrast (d-f) considers that the\ninteraction between the particles favors the state with higher occupation, i.e. reduces the exchange rate. This would be the case\nwhen collective dynamics self-stabilizes the system in a polarized state. Let’s assume that by thermal fluctuations one particle\ncrosses over the barrier (d) and introduces non-equal occupation of the states. The interaction between the particles and the\ndiscrepancy in the occupation between the two states introduces the energy difference between the two minima E and more\npopulated state is slightly more favored (e). Lastly, if we would polarize such interacting system in an external field (as illustrated\nin (b)) and switch the external field off, being significantly more populated one state would be significantly more favored (f). In this\ncase the exchange rates between the two states would significantly differ, and to de-polarize the system a much larger effective\nbarrier would be observed.\nmemory computing and synaptic circuit applications.\nAuthors’ response: We thank the reviewer for raising this valuable point, which is actually one of the directions that our team is\ncurrently working on. We envision the non-volatile field-programmable nature of the edge adsorbed water switch as promising for\ncomputing applications in two configurations:\nFirst, as robust memristors with non-volatile synaptic weights. This could be achieved by maximizing the hysteresis effect in a dual\ngate (top and bottom) geometry architecture with high-k gate dielectrics to enhance the V range while reducing operating voltages\nH\nV . Further engineering of narrower and pitch-controlled ribbons would enhance the on/off ratio, allowing memristive switching as\nSG\nhighly tunable synaptic elements controlled by the edge water dipoles.\nTop gate\nWater-terminated Gr-NR\nsource drain\nBottom gate\nExample structure of a dual-gate Gr-NR FET.\nSecond, one could exploit the nanoribbon network as a programmable ferroelectric electrode in a vertical field effect transistor\n(V-FET) configuration with an ambipolar 2D semiconductor as a channel. In this case, the dipolar field of the water molecules would\nnon-volatily dope the 2D semiconductor channel, allowing switching between a p-type, n-type, and the off state, making possible\na device that acts as a multi-state memory cell and a logic gate. A network of these could make a truly reconfigurable circuit, which\nis a key goal for in-memory computing.\nsource\ndrain\n2D semiconductor channel\nWater-terminated Gr-NR electrode\nGate dielectric\ngate\nExample structure of a V-FET with Gr-NR electrode\nChanges made: The possible geometries for future integration of the water switch into computing-in-memory elements are briefly\naddressed in the Discussion section of the revised manuscript.\nR3Q6: There are minor typographical errors in the manuscript, such as \"VH\" on line 87 not being italicized. Please ensure that all\ncontent adheres to standard formatting conventions.\nAuthors’ response: We thank the reviewer for the careful reading of our work and for noticing this omission. We have performed\na thorough review of the entire manuscript and SI to correct all typos and ensure consistent formatting of all variables, constants,\nand abbreviation definitions.\nChanges made: Corrected typographical errors are highlighted in the marked-up versions of the revised manuscript and SI by red\nfont color.\nReviewer #1 (Remarks to the Author):\nThe Authors exhaustively replied to all raised concerns and added significant new data and extra\ndiscussion in the revised version of the manuscript. I therefore recommend the publication of\nthis manuscript as is.\nAuthors’ response: We’re thankful to Reviewer #1 for their time and effort in evaluating this\nwork and for the positive feedback following the revisions.\nReviewer #2 (Remarks to the Author):\nBecause transport experiments probe the overall average electrical properties of a sample, it is\ninherently challenging to disentangle contributions that arise solely from the edge. As a result,\narguments by authors cannot yet be regarded as completely conclusive. However, it is also true\nthat the authors provided reasonable answers and persuasive arguments to the various\nconcerns I raised.\nIn my view, the manuscript is suitable for publication at this stage, allowing the broader\nreadership to assess its merits.\nAuthors’ response: We’re thankful to Reviewer #2 for their critical feedback and evaluation of\nour work, which helped make it significantly better.\nReviewer #3 (Remarks to the Author):\nThe authors have addressed my comments clearly, I have no more questions.\nAuthors’ response: We’re thankful to Reviewer #3 for their time and effort in evaluating this\nwork and contributing to its increased quality and clarity in the revised version."}, {"round": 2, "reviewer_comments": {"reviewer_1": "The Authors exhaustively replied to all raised concerns and added significant new data and extra discussion in the revised version of the manuscript. I therefore recommend the publication of this manuscript as is.", "reviewer_2": "Because transport experiments probe the overall average electrical properties of a sample, it is inherently challenging to disentangle contributions that arise solely from the edge. As a result, arguments by authors cannot yet be regarded as completely conclusive. However, it is also true that the authors provided reasonable answers and persuasive arguments to the various concerns I raised.\nIn my view, the manuscript is suitable for publication at this stage, allowing the broader readership to assess its merits.", "reviewer_3": "The authors have addressed my comments clearly, I have no more questions."}, "author_response": "Authors’ response: We’re thankful to Reviewer #1 for their time and effort in evaluating this work and for the positive feedback following the revisions.\nAuthors’ response: We’re thankful to Reviewer #2 for their critical feedback and evaluation of our work, which helped make it significantly better.\nAuthors’ response: We’re thankful to Reviewer #3 for their time and effort in evaluating this work and contributing to its increased quality and clarity in the revised version."}]

Decision letter