# Fileset

[s41467-025-65922-6.pdf](https://mdr.nims.go.jp/filesets/14b9ce3e-ca2b-47ec-94c7-81326217d9eb/download)

## Creator

Muhammad Awais Aslam, Igor Stanković, Gennadiy Murastov, Amy Carl, Muhammad Zubair Khan, Zehao Song, [Kenji Watanabe](https://orcid.org/0000-0003-3701-8119), [Takashi Taniguchi](https://orcid.org/0000-0002-1467-3105), Alois Lugstein, Christian Teichert, Roman Gorbachev, Raul D. Rodriguez, Aleksandar Matković

## Rights

[Creative Commons BY Attribution 4.0 International](https://creativecommons.org/licenses/by/4.0/)

## Other metadata

[Ferroelectricity in graphene nanoribbon devices enabled by collective water molecule dynamics](https://mdr.nims.go.jp/datasets/b70acf65-5de4-45e7-8ab6-47f03284df87)

## Fulltext

Ferroelectricity in graphene nanoribbon devices enabled by collective water molecule dynamicsArticle https://doi.org/10.1038/s41467-025-65922-6Ferroelectricity in graphene nanoribbondevices enabled by collective watermoleculedynamicsMuhammad Awais Aslam1,9, Igor Stanković 1,2,9 , Gennadiy Murastov 1,Amy Carl 3,4, Muhammad Zubair Khan1, Zehao Song 5, Kenji Watanabe 6,Takashi Taniguchi 7, Alois Lugstein 5, Christian Teichert1,Roman Gorbachev 3, Raul D. Rodriguez 8 & Aleksandar Matković 1Water is omnipresent in nanoscale systems, yet its collective dynamics andimpact on emerging electronics remain poorly understood. Here, we investi-gate the role of water molecule dynamics in the ferroelectric response ofgraphene nanoribbon devices. Our findings demonstrate that the collectivedynamics of water molecules stabilize the ferroelectric effect. We find that aminimum bi-layer thickness is required for the temperature stability of theferroelectric effect. In contrast, mono-layer ribbons show a 70% shrinkage ofthe hysteresis window between 120 and 400K. Using a combination of elec-trical transport measurements and molecular dynamics simulations, we con-clude that water molecules bridging between graphene nanoribbon layersstabilize the formation of water clusters via intermolecular Coulomb interac-tions, driving a robust ferroelectric behavior and remnant polarizationobserved at the device level. This work lays the foundations for exploitingwater dynamics in next-generation ferroelectric heterostructures, with directimplications for neuromorphic computing and memory devices.Water plays essential roles in a wide range of phenomena, frommaterials and electronics to life itself1,2. Understanding the interactionof water with low-dimensional materials has proven essential fornanofluidics3, energy storage4,5, water splitting6, purification7,8, andwater-assisted ferroelectricity9–11. Water plays a critical role in tuningmaterial properties at the nanoscale. It turns a challenge into oppor-tunities for nanoelectronics bymodifying electronic transport, chargedistribution, and doping, making it a critical factor for applications inneuromorphic computing and molecular-level devices, especiallywhen considering the impact of the edge-adsorbed species on gra-phene nanoribbons (Gr NRs)12–16.Graphene nanoribbons with a high edge-to-surface ratio providean excellent platform to study the influence of water dynamics at theedge-water molecule interface. Previous work showed that watermolecules adsorbed at graphene edges in single-layer microdevicesencapsulated by hexagonal boron nitride (hBN) can induceferroelectricity10, a phenomenon also observed in nanoribbon-baseddevices11. This phenomenon originates from the dipolar nature ofwater that can result in switching between two polarization statesunder external electric fields. However, the exact mechanism behindwater-induced polarization in Gr NR devices, the impact of field-related dynamics of edge-adsorbed water, the influence of grapheneReceived: 2 May 2023Accepted: 27 October 2025Check for updates1Chair of Physics, Department Physics, Mechanics, and Electrical Engineering, Montanuniversität Leoben, Leoben, Austria. 2Scientific Computing Laboratory,Center for the Study of Complex Systems, Institute of Physics Belgrade, University of Belgrade, Belgrade, Serbia. 3Department of Physics and Astronomy andNational Graphene Institute, University of Manchester, Manchester, United Kingdom. 4Department of Materials and National Graphene Institute, University ofManchester, Manchester, United Kingdom. 5Institute of Solid State Electronics, TU Wien, Vienna, Austria. 6Research Center for Electronic and OpticalMaterials, National Institute forMaterials Science, Tsukuba, Japan. 7ResearchCenter forMaterials Nanoarchitectonics, National Institute forMaterials Science,Tsukuba, Japan. 8Tomsk Polytechnic University, Tomsk, Russia. 9These authors contributed equally: Muhammad Awais Aslam, Igor Stanković.e-mail: igor@ipb.ac.rs; raul@tpu.ru; aleksandar.matkovic@unileoben.ac.atNature Communications |        (2025) 16:10982 11234567890():,;1234567890():,;http://orcid.org/0000-0001-5756-7196http://orcid.org/0000-0001-5756-7196http://orcid.org/0000-0001-5756-7196http://orcid.org/0000-0001-5756-7196http://orcid.org/0000-0001-5756-7196http://orcid.org/0000-0002-7681-3118http://orcid.org/0000-0002-7681-3118http://orcid.org/0000-0002-7681-3118http://orcid.org/0000-0002-7681-3118http://orcid.org/0000-0002-7681-3118http://orcid.org/0000-0003-3601-1636http://orcid.org/0000-0003-3601-1636http://orcid.org/0000-0003-3601-1636http://orcid.org/0000-0003-3601-1636http://orcid.org/0000-0003-3601-1636http://orcid.org/0000-0003-3276-7516http://orcid.org/0000-0003-3276-7516http://orcid.org/0000-0003-3276-7516http://orcid.org/0000-0003-3276-7516http://orcid.org/0000-0003-3276-7516http://orcid.org/0000-0003-3701-8119http://orcid.org/0000-0003-3701-8119http://orcid.org/0000-0003-3701-8119http://orcid.org/0000-0003-3701-8119http://orcid.org/0000-0003-3701-8119http://orcid.org/0000-0002-1467-3105http://orcid.org/0000-0002-1467-3105http://orcid.org/0000-0002-1467-3105http://orcid.org/0000-0002-1467-3105http://orcid.org/0000-0002-1467-3105http://orcid.org/0000-0001-5693-4775http://orcid.org/0000-0001-5693-4775http://orcid.org/0000-0001-5693-4775http://orcid.org/0000-0001-5693-4775http://orcid.org/0000-0001-5693-4775http://orcid.org/0000-0003-3604-5617http://orcid.org/0000-0003-3604-5617http://orcid.org/0000-0003-3604-5617http://orcid.org/0000-0003-3604-5617http://orcid.org/0000-0003-3604-5617http://orcid.org/0000-0003-4016-1469http://orcid.org/0000-0003-4016-1469http://orcid.org/0000-0003-4016-1469http://orcid.org/0000-0003-4016-1469http://orcid.org/0000-0003-4016-1469http://orcid.org/0000-0001-8072-6220http://orcid.org/0000-0001-8072-6220http://orcid.org/0000-0001-8072-6220http://orcid.org/0000-0001-8072-6220http://orcid.org/0000-0001-8072-6220http://crossmark.crossref.org/dialog/?doi=10.1038/s41467-025-65922-6&domain=pdfhttp://crossmark.crossref.org/dialog/?doi=10.1038/s41467-025-65922-6&domain=pdfhttp://crossmark.crossref.org/dialog/?doi=10.1038/s41467-025-65922-6&domain=pdfhttp://crossmark.crossref.org/dialog/?doi=10.1038/s41467-025-65922-6&domain=pdfmailto:igor@ipb.ac.rsmailto:raul@tpu.rumailto:aleksandar.matkovic@unileoben.ac.atwww.nature.com/naturecommunicationsthickness on water anchoring, temperature-dependent stability, andthe interplay of these factors on ferroelectricity remain poorlyunderstood. These knowledge gaps limit the development of next-generation devices based on water-induced ferroelectricity, includ-ing applications in radio frequency switches, neuromorphic com-puting, and memcapacitors17–21. This study builds upon the existingknowledge of water-induced ferroelectricity in graphene edges10,11. Itbrings forward new experimental results demonstrating how col-lective water dynamics are critical to this effect through graphenelayer number and temperature dependence. These results are sup-ported by molecular dynamics simulations describing the collectivebehavior of water at graphene edges. Our twofold experimental andcomputational approach shows that this ferroelectric effect at thedevice level in Gr NRs integrated into field-effect transistors (FETs)strongly depends on the graphene layer number, requiring a mini-mum bi-layer thickness to sustain temperature-stable ferroelec-tricity. We elucidate the physical mechanism behind these results asdriven by the collective behavior of water that induces edgeanchoring and interlayer bridging by molecular water necks thatenhance intermolecular Coulomb interactions. By investigating theinterplay between external fields, temperature, layer number, andadsorbed water at graphene edges, we open the door to futurenanoelectronics and sensor applications leveraging the ubiquitousnature of water confined at the nanoscale22–29.ResultsObservation of the hysteretic response in oxygen-terminatedgraphene nanoribbon devicesFigure. 1a depicts the schematic of a Gr NR FET on hBN. Para-hexaphenyl (6 P) organic nanoneedles self-assembled on hBN-Gr-hBNheterostructures serve as self-aligned masks11,24,30,31. Reactive ionetching (RIE) selectively removes exposed graphene to produce Grnanoribbons, as illustrated in Fig. 1b, c; detailed in theMethods sectionand Supplementary Fig. S1. The resulting nanoribbons have an averagewidth of 20 nm, as estimated from high-resolution atomic forcemicroscopy (AFM) images (Supplementary Fig. S2). Alternative toorganic nanostructures, recent studies have shown that directed self-assembly of block copolymers can produce complex arrays of highlyordered sub-10 nm wide nanopatterns32–34, that can be used to fabri-cate 2D material nanoribbon FETs35.The top and bottom encapsulation with hBN insulates graphenefrompotential charge-trap sources, ensuring that external interactionsprimarily occur at the nanoribbon’s edges (Supplementary Fig. S3).Organic nanostructure growth is done near thermodynamicfixed bridgewater moleculesmobileclusterH: O: C: 2L3L5L4Lgraphenenanoribbons(a)drainsourcehBNSiO2/Siglobal back gate(b)(c)(d)graphenegraphenenanoribbonnanostructureorganichBNhBN(W/L). Rtot   [k�]0.00.20.40.60.81.01.21.41.61.8 hysteresiswindow300 K-30 -20 -10 0 10 20 30 40VSG - VCNP [V](e) (f)SF6 etchO2 etch150 200 250 300 350 400T [K]mono-layers0.20.40.60.81.01.2(VH)   (arb. units)normalized hysteresis window 150 200 250 300 350 400T [K]0.20.40.60.81.01.2(VH)   (arb. units)multi-layersEa=(17.45+11.60)meVEa>150meVFig. 1 | Device structure and the observation of the hysteretic response.a Schematic representation of a nanoribbon field-effect transistor. b, c Opticalimages and corresponding sketches of stacked hBN/Gr/hBN flakes and patternednanoribbons. d Width and length (W L1) scaled total device resistance(two–terminal) transfer characteristics for 2 LGr NR FETs (measured under 2·10−2mbar). Arrows indicate the sweepdirection. The VSG range is shiftedwith respect tothe mean VCNP values between both sweeping directions (considering the shifts of3 V and 28V for oxygenated and fluorinated ribbons, respectively).e, f Temperature dependence of the hysteresis windows (VH) normalized to thevalue recorded at 120K, presenting thedata for themono-layer andmulti-layer (2–5layers, with thickness labelled at the right side for each curve in (f)) Gr NR FETs.Different lines represent different devices, and error bars depict VH variation in thesubsequent measurements. Shaded areas serve as a guide to the eye. Solid blacklines are the Arrhenius function-based analytical fits used to extract the energybarrier (Ea). Dashed black lines in (e) correspond to one sigma deviations of the Eaestimate for monolayers. Insets in (e) and (f) present MD simulation models exhi-biting typical local configurations of water clusters in the absence of the externalelectric field.Article https://doi.org/10.1038/s41467-025-65922-6Nature Communications |        (2025) 16:10982 2www.nature.com/naturecommunicationsequilibrium conditions, meaning that even minor changes in the sub-strate can significantly impact self-assembly. The interfacial bubblesintroduce bi-axial tensile strain in the top hBN layer, and the organicnanostructure self-assembly is hindered in such strained regions. As aresult, the nanoribbon network effectively “avoids” these bubbles(Supplementary Fig. S3d–g).The nanoribbon edge termination is controlled by the choice ofprecursor 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-terminatededges exhibit pronounced and stable hysteresis, with switchable p andn doping characteristics depending on the gate bias (VSG) sweepingdirection (Fig. 1d). The electrical response comparison between thenanoribbons and the starting flake FETs is shown in SupplementaryFig. S4. The hysteresis window (VH) is defined as the differencebetween two charge neutrality points (CNPs), which correspond to themaxima in the total two-terminal device resistance Rtot(VSG). In con-trast, F-terminated ribbons show negligible hysteresis due to theirhydrophobic nature (Fig. 1d, blue plots), as previously observed forgraphene flakes10. Given that oxygen-termination enables watermolecule adsorption, our focus will be on oxygen-terminated gra-phene nanoribbons and their water-induced ferroelectricity10,36,37. Toensurewatermolecule adsorption on the oxygenated ribbon edges, GrNR devices were exposed to ambient conditions (relative humidity of22−26 % at 298K) for at least 30min following the RIE step. This pro-cedure was also applied to the SF6-etched devices to maintain com-parability between the two datasets.We have used a total of 22 devices, comprising 16 O2 RIE-treatedand 6 SF6 RIE-treated devices. To illustrate sample-to-sample varia-tion in performance parameters, we present peak apparent linearmobility (µlin) statistics in Supplementary Fig. S5. Furthermore,µlin(VSG) curves for both SF6 and O2 etched devices were comparedbefore and after the RIE step, see Supplementary Fig. S6. Theassembly of the nanostructure network varies between samples,impacting the scattering of ribbon widths. Consequently, the VHvalues show a scatter of about ± 25 % among different devices undercomparable operating conditions. This scattering of VH values isillustrated in Supplementary Fig. S7, as a function of the layer num-ber and the total sweeping range of the gate voltage (VSG-tot). Furtherdiscussion regarding the relationship between ribbonwidth and VH isprovided later in the text from the perspective of moleculardynamics (MD) simulations.The details of the mechanism behind the hysteretic character-istics presented in Fig. 1d are known, as Caridad et al. proposed thatsingle molecules switch between the two states10. The external gatefields disrupt the equal probability of water arrangements whileinducing a torque on the water adsorbed at the edges. The torqueresults from Coulomb forces acting on the water dipole, allowing itto cross the energy barrier between the two states ca.25–40meV10,38,39.Interestingly, Rtot(VSG) hysteresis and its thermal stability exhibitprofound differences depending on whether the oxygenated edgesof the ribbons are mono-layered or multi-layered. Fig. 1e, f comparesthe temperature dependence of the normalized hysteresis windowfor mono-layer and multi-layer (2–5 layers) Gr NR FETs over a tem-perature range of 120–400K. The VSG sweeping range was fixedbetween − 40 V and + 40 V for each VH value, with a constantsweeping rate of 2 V s−1.In mono-layer Gr NR FETs, the hysteresis window decreases astemperature increases, although significant values persist across thewhole temperature range. In contrast, multi-layer devices with two ormore graphene layers (L ≥ 2) show a ferroelectric effect with negligiblevariations across the investigated temperature range (see Fig. 1f). Formulti-layer devices (2 L to 5 L), we did not observe any significant dif-ferences in the response of the Gr NR FETs with different layerthickness. However, there is a trend indicating that both VH and µlintend to decrease as the number of layers increases (SupplementaryFigs. 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 analyticalmodel based on the Arrhenius equation. Thismodelconsiders the activation energy (Ea) or the energy barrier between thetwopolarization states, expressed as thefitting constant in:VH(T) =VH0/ (1 - exp(-Ea/(kBT)). Here, kB is the Boltzmann constant, T is the tem-perature, and VH0 scales the fit to the value obtained at the lowesttemperature. For the monolayers, we obtained an activation energy ofEa = (17.45 ± 11.60) meV. The fits for the individual VH(T) curves of themono-layer GrNRN FETs can be found in Supplementary Fig. S8. Incontrast, the data for the multilayer devices in Fig. 1f shows an energybarrier higher than 150meV.These observations led us to hypothesize that a single watermolecule switching between the sides of the graphene plane isresponsible for the ferroelectric response in mono-layer devices,leading to a temperature-dependent behavior. When the thermalenergy of the molecule exceeds the switching energy barrier, themolecules can thermally transition between states, which explains thetemperature-dependent results in Fig. 1e. Additional results regardingthe dependence of the hysteresis window on the VSG sweeping rangecan be found in Supplementary Fig. S9.Stability of the hysteretic response: excluding charge traporiginsThe direction of the observed hysteresis could originate from thecharge traps in the gate dielectric layer. However, the same hysteresisdirection could also result from the edge dipoles, as illustrated inSupplementary Fig. S10. To further verify the nature of the observedhysteresis, we performed several experiments to rule out the influenceof extrinsic effects such as impurities and charge traps. Since trappingcanoccur on the timescale of seconds40,41, we used varying sweep rates(from0.5 V s−1 to 20V s−1) to evidence its effects onVH. Figure. 2a showsthat VH as a function of the sweeping rate remains largely unperturbedat 298K and 77 K under forward and backward gate sweeps, suggest-ing 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. Interest-ingly, near the charge neutrality point, the resistance exhibits aper-iodic fluctuationswith high reproducibility, indicating the opening of atransport gap. While such fluctuations could be attributed to quantuminterference phenomena44,45, Chen et al.46 showed that such insulatingstates can be attributed to the adsorption of water molecules on thenanoribbon edges, generating strong electric fields that modulate theNRs’ band gap.To further demonstrate the impact of watermolecule adsorption,the devices were vacuum annealed at 520K for 60min at a pressure of8·10−3mbar. Subsequent electrical measurements were performed at300K, avoiding any ambient exposure. As shown in Fig. 2c, d, thehysteresis window VH drastically decreases after annealing, consistentwith the desorption of water molecules from the nanoribbon edges.Remarkably, the hysteresis is restored upon exposure to ambientconditions (relative humidity of 22–26 % at 298K), providing strongevidence for the critical role of water adsorption in the ferroelectricresponse. In addition, the unintentional p–type doping evidenced inFig. 2c implies thatonly positiveVSG values are required to compensatefor built–in charge transfer doping before the total field direction canreverse and switch the water dipoles orientation (SupplementaryFig. S11). The dependence ofVHon the vacuumconditions is detailed inSupplementary Fig. S12, where we observed a 40% increase in VH whengoing from vacuum conditions to ambient pressure at 300K undercontrolled relative humidity of 40%. Furthermore, VH remained con-stant at pressures below 101mbar.Article https://doi.org/10.1038/s41467-025-65922-6Nature Communications |        (2025) 16:10982 3www.nature.com/naturecommunicationsMolecular dynamics simulations of the oxygenated and water-terminated graphene nanoribbonsWe performed MD simulations to understand better the role of watermolecules in inducing electric fields and to elucidate the mechanismbehind the temperature-dependent hysteresis differences betweenmono-layer and multi-layer Gr NR devices. The tendency of the watermolecules to bind at oxygenated edges and form clusters39 suggeststhat a collective behavior stabilizes the molecules, resulting in atemperature-independent effect. Two conditions must be met toachieve such stable ferroelectric ordering of water molecules on Gredges: (i) a low thermal barrier that allows a single polar molecule toswitch between states, and (ii) a cluster large enough to be stabilizedby intermolecular Coulomb interactions. While the first condition iswell established, we have designedmolecular dynamics simulations toinvestigate the second condition. The simulations implicitly considerthe hydrogen bonds within the dynamically calculated local chargesbetween water molecules and the oxidized graphene edge. Moreover,we also calculate the electric field, which is a direct consequence of themolecule’s polarization, as shown later.Our calculations agreewith themodel presented in ref. 10, and provide new insights into the impor-tance of the collective behavior of water at multi-layered edges.To conceptualize the mechanistic difference between the twosystems, we can think of bothmono- and bi-layers as a bi-stable systemwith two states and a barrier between them. The collective interactionthat stabilizes the cluster in one state could essentially beviewed as thedeepening of oneminimum, depending on howmanymolecules are inthat minimum. This introduces effectively a discrepancy in the rate ofreturning to the lower occupied state by thermal fluctuations. Theprocess is illustrated in Supplementary Fig. S13.Therefore, the existence of water interaction and the size of thecluster play a key role in achieving the temperature-independenthysteretic effect. Naturally, there is a limit to how large a cluster can bewhile remaining anchored at the edge. As the cluster size increases, theinduced dipole field also grows. At the same time, the energy requiredto detach the molecules furthest from the edge decreases (Supple-mentary Fig. S14). These two effects compete. In the presented simu-lations, we have limited the amount of watermolecules per unit lengthof the edge to below one-half of the limit determined by the clusterstability in strong external fields at 300K (more details in the “Meth-ods” section).The MD simulations of mono-layer graphene confirm the pre-dicted energy barrier of 25−40meV for switching a single watermolecule10,38. In contrast, the multi-layer system shows a fundamen-tally different behavior. The layer separation in graphite allows watermolecules to form a bridge between two oxygen-terminated edges.These bridging molecules remain stable and do not react to externalelectric fields, yet they promote the formation of a water clusteraround them. MD simulations show that anchored water moleculeswithin the cluster can sustain high electric fields (up to 5 V nm−1) andtemperatures up to 500K before desorbing from the edges. Thecluster surrounding the bridging molecules is mobile in the electricfield, and its collective behavior helps stabilize the structure, resultingin a hysteresis window that is essentially temperature-independent.This collective self-stabilizing behavior of polar objects resembles thatpreviously observed in colloidal systems and molecular motors47,48.Our calculations do not rule out the formation of water clusters inmono-layers. In mono-layers, flipping of individual molecules acrossthe basal plane of the ribbons is frequently observed, even in thepresenceof clusters.However, such transitions of individualmoleculesare essentially absent in the bi-layer case. The distribution of watermolecules along the edges in both the 1 L and 2 L cases is presented inSupplementary Fig. S15.Expanding on our findings, we investigated how varying electricfields affect water clusters using MD simulations. An externalhomogeneous electric field generated by infinite planes was appliedto 1 L and 2 LGr NRs, as described in the Methods section. Thefraction of polarized molecules above (n+) and below (n-) the gra-phene layer quantifies the ferroelectric effect arising frommolecularswitching. We observed charge bi-stability and switching betweentwo states under a cyclical change of the external field. The resultingdifference between n+ and n- per unit length of the NR edge is pre-sented in Fig. 3a, b, comparing the hysteresis loops obtained formono- and bi-layer Gr NRs. It is important to note that in theexperiment, the electric field source was generated between theribbon and the flat surface (gate electrode). In contrast, in thesimulations, the electro-neutral system with the ribbon and watermolecules was placed between two uniformly charged planar elec-trodes, generating a homogeneous external out-of-plane field. Fur-ther, in the simulations, the electric field applied was over an order ofmagnitude higher than in the experiment when considering a parallelcapacitance model. This was done to reduce simulation time sincethe electric field strength exponentially prolongs the switching timeof the model system. Field enhancement is expected in our experi-ments due to fringing capacitance effects10. An estimate of theenhancement and the field profiles based on a finite elements modelof the device structure is provided in the Supporting Information,Supplementary Fig. S16. We anticipate about one order of magnitudeenhancement of the out–of–plane field component due to theFig. 2 | Exclusion of the charge-trap origin of the hysteretic response. a Changein the VH with respect to varying VSG sweeping rate (2 L device, measured under2·10−2 mbar). b Total device resistance (two–terminal; semi–log curves) transfercurves demonstrating that the hysteresis prevails at 4.2 K with the emergence of atransport gap, which is highlighted in red color (1 L device, measured under 4·10−6mbar). c, d Total device resistance (two–terminal) transfer curves of the samedevice before and after annealing at 520K under vacuum (1 L device, measuredunder 2·10−2mbar, 300K). Red lines next to the x–axis in (b–d) indicate estimates ofthe unpolarized CNP values. Arrows indicate the sweep direction.Article https://doi.org/10.1038/s41467-025-65922-6Nature Communications |        (2025) 16:10982 4www.nature.com/naturecommunicationsfringing capacitance effects, reaching up to 2 V nm−1 within theexperimentally achievable biasing range.For amono-layer ribbon, a weak hysteresis in n+ - n- is observed asa function of the external electric field. Moreover, this hysteresisdecreases with increasing temperature (as shown by the orange circlesin Fig. 3c), indicating a low energy barrier between the states. Thisphenomenon can be attributed to the absence of a bridging watermolecule found in multi-layer systems. Our experiments with mono-layer Gr NRs showed a more robust hysteresis, especially at elevatedtemperatures. The discrepancy between the modeled and experi-mentally observed mono-layer systems may be attributed to morecomplex interfaces present in the actual experiments that are notcaptured in the simulations. The model considers a free-standinggraphene ribbon, while the experiment involves hBN encapsulation,which likely creates complex electrostatic environments that couldalso partially enhance the collective behaviour of the adsorbed waterclusters. In addition, the model system has a considerably shortertemporal evolution (4 ns) in response to the external electrostaticperturbations compared to the experimental timescale of over 100ms.This shorter simulation time scale may prevent the model from ade-quately relaxing within the achievable 4 ns window, especially attemperatures below 250K. In contrast, the experimental system hasmore time to achieve polarization and to experience decoherence ofthe polarized state once the external field is turned off. Consequently,direct quantitative comparisons of the remnant field values betweenthe experiment and the model are not feasible. Nevertheless, themodel effectively captures the trends in temperature and layer num-ber observed in the experiments, allowing us to elucidate the physicalmechanisms behind the ferroelectric behavior of Gr NRs devices.For a bi-layer NR, a prominent hysteresis was noted in the MDsimulations (Fig. 3b), and the opening of the hysteresis loopwas foundto be temperature-independent in the probing range from 250–450K,as shown by purple squares in Fig. 3c. The field lines and field strengthacting on a 2 L nanoribbon generated via MD simulations are pre-sented in Supplementary Fig. S17. As the electric field decreases, somemolecules remain stationary until the field polarization switches,causing these molecules to migrate to the other side with respect tothe nanoribbon’s basal plane.The water cluster at the edge of a 2 L NR is stable in the externalelectric field up to 4.1 V nm−1, after which themolecules were observedto be removed from the cluster by the external field at 300K. n+ - n-hysteresis presented in Fig. 3d, uses the external field range up to theFig. 3 | Molecular dynamics calculation results. a, b Obtained from the MDsimulations, respectively, for mono- and bi-layer nanoribbon, the count of thedifference of the water dipoles above and below the graphene basal plane (n+ - n-)per unit of the edge length, as a function of the applied external electric field. Solidred and dashedblack hysteresis loopspresent the responses of the systems relaxedat 400K and 300K, respectively. c Temperature dependence of the opening of then+ - n- hysteresis as a function of the temperature, considering the external fieldsweeping range of ± 2 V nm-1. Circles indicate mono-layer and squares bi-layernanoribbons. Solid lines present the mean linear fits, and the dashed area serves asa guide to the eye. Error bars present the standard deviation from five differentcurves. d Example of the bi-layer nanoribbon n+ - n- hysteresis with increasedexternal field sweeping range up to ± 4.1 V nm−1. Points labeled (e–g) in (d) have thecross-sections of the resulting model structure presented at the top of the sub-panels (e–g). The respective density of charges, as well as electric field lines, areshown to scale below. The field in (e–g) follows a downwards sweep direction from4.1 V nm−1 to 4.1 V nm−1, and depicts the collective ordering of water at (e)Eext = 4 V nm−1, (f) Eext = 0.8 V nm−1, and (g) Eext = −4 V nm−1. f presents the coercivitypoint, in which the nanoribbon experiences approximately a net zero out-of-planedipolar field from the water clusters adsorbed at the edges.Article https://doi.org/10.1038/s41467-025-65922-6Nature Communications |        (2025) 16:10982 5www.nature.com/naturecommunicationslimit of the desorption of themolecules from the cluster. In this case, aclear indication of saturation is observed. Such high fields probed inthese MD simulations are not reachable in the experiment due to theSiO2 and hBN field breakdown. To visualize the ordering of the watercluster, we follow its evolution along a down-sweep of the field in thehysteresis shown in Fig. 3d. Figure 3e–g represents snapshots of theedge segments for a 2 LGr NR system with the corresponding chargedistribution and induced electric field lines (below). The system con-sists of both externally applied fields (not shown) and induced fieldswith opposite polarities in the z-direction. Figure 3e presents the caseof saturation with an external field of 4 V nm−1 applied in the z-direc-tion. An induced charge redistribution is observed due to the mole-cular polarization at the edges, which in turn creates an effectiveelectric dipole moment. This changes the local field over the wholegraphene bi-layer profile. The effective electric field is, therefore, asum of the two fields, and as the external field decreases, they couldcancel each other in parts of the graphene plane. Fig. 3f displays theconfiguration where the induced electric dipole points roughly alongthe basal plane of Gr NR, i.e., representing one of the coercivity pointsof the out-of-plane dipole field hysteresis. Consequently, the inducedelectricfield component perpendicular to the ribbonwill be negligible.Finally, for − 4 V nm−1, we observe a complete reversal of the inducedelectricfield (Fig. 3g). As thefield strength increases, themigration andmolecular polarisation repeat in the opposite order, generating ahysteresis loop. The bonding between the polar molecular ensembleand the graphene edge, together with the intermolecular Coulombinteractions, should be strong enough to prevent the external electricfield from tearing off molecules from the cluster (SupplementaryFig. S14b). The details of the large-scale atomistic model are presentedin Supporting Information Supplementary Fig. S18.The width of the ribbons significantly impacts the strength of thedipole-induced fields. Through MD simulations of polarized dipoles,for both 1 L and 2 L ribbons, we have extracted the dipole-inducedfields perpendicular to the basal plane of graphene and their depen-dence along the direction of the ribbon width. As shown in Supple-mentary Fig. S19, the field induced by a single polarized edge followsthe theoretically expected x-2 dependence. We hypothesize that anoptimalwidth of about 10 nmmaximizes the intensity of the hysteresiseffect, arising from two competing processes. For ribbons that are toonarrow, the field induced by one row of polarized edge dipoles beginsto affect the dipoles at the opposite edge, resulting in an antiparallelstate. For widths below 10 nm, the field from one polarized edge isexpected to impact the other edge, overcoming external fields forwidths below5 nm. In the other extreme, as the ribbonwidth increases,the electrical response begins to show more bulk-like properties ofgraphene, resulting in a reduction of hysteresis. Our MD simulationssuggest that forwidths above 100nm, the contribution from the edgesto the electrical response becomes negligible.Probing the remanence of the water dipole-induced fields ingraphene nanoribbonsTo confirm the pronounced hysteresis predicted by the MD simula-tion, we tracked the drain current (ID) evolution of Gr NR FETs withoutan 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-terminatedGrNRFET,where the range of the gate voltageswas chosenso that one of the charge neutrality points is at VSG = 0V49. The asym-metric VSG sweep approach was used for CNP position adjustment(Supplementary Fig. S20). This way, the difference in ID for the n+ andthe n- states is maximized. The two ID(VSG = 0V) states can be definedas “0” and “1”, corresponding to the low or the CNP state and the highcurrent state, schematically illustrated for a bi-layer in Fig. 4b.By sweeping from VSG = 0V to one end of the range and returningto VSG = 0 V, the system is set either to the ID “0” or “1” state. From thispoint on, the ID is recorded as a function of time without an externallyapplied gate field. The dipole-induced remanent field was observed inboth mono- and multi-layer Gr NR FETs. The results comparing amono- and a multi-layer Gr NR FET are shown in Fig. 4c, where thenormalizeddifferencebetween ID1(t) and ID0(t)was tracked for 600 s at300K 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 scalingconstant and τ for a time constant.The remanence of the field for a mono-layer is less pronounced,and the current differencebetween the two states decays faster than inthe multi-layer Gr NR FETs (Fig. 4c). The time constants indicate thatthe 1 LGr NR devices lose 90 % of the initial ID difference between the“0” and “1” states in about 30minutes and multi-layer Gr NR FETs inabout 90minutes. Further, a faster decaying component withτ = (64 ± 32) s, responsible for about 15–20% of the initial drop of the IDdifference, was noticed in all devices and attributed to weakly boundmolecules. Gr NR FETs with L ≥ 2 require several hours of relaxation tohave the initial transfer curve sweep starting from the depolarizedstate of the water molecules, which is consistent with the extractedtime constants.DiscussionOur study demonstrates the critical role of water dynamics in stabi-lizing ferroelectric effects in graphene nanoribbon FETs. A key finding(a) (b) (c)energyenergyzz“0” “1”-40 -20 0 20 40“0”“1”VSG  [V]VSD  = 10 mVset “1” set “0”406080200I D  [nA]( ID1(t) - ID0(t) )/( I D1(0) - I D0(0) )0 100 200 300 400 500 600t [s]1 L Gr-NRGr-NRVSG = 0 V0.40.60.81.04 LWater-induceddipole switchFig. 4 | Probing remanence of the hysteretic response. a The hysteresis curve for1 LGr NR FET after the application of a bias stress (measured under 2·10−2 mbar;unpolarizedpositionof theCNP∼ 2 V).bThe schematic representationof retentionstates for the water clusters adsorbed at the edges. c The source-drain currentdeviation in time with respect to the 1 L τ = (1041 ± 315) s, and 4 L τ = (3020 ± 125) sGrNR devices (measured under 2·10−2 mbar at 300K).Article https://doi.org/10.1038/s41467-025-65922-6Nature Communications |        (2025) 16:10982 6www.nature.com/naturecommunicationsis the difference in temperature dependency of the hysteresis windowbetweenwater-terminatedmono-layer andmulti-layer nanoribbons. Inbi-layer or thicker NRs, the ferroelectric response remains essentiallytemperature-independent up to 450K due to a collective self-stabilizing effect. This behavior requires a large water cluster, allow-ing molecules to remain bound in one state through intermolecularCoulomb interactions. Molecular dynamics simulations reveal themechanismbehind thedifferent ferroelectric behavior betweenmono-and multi-layer NRs. In multi-layer nanoribbons, the presence of twoadjacent layers enables the formation of a water molecule bridge that,along with the oxygen-terminated edges, promotes the formation of awater cluster.Ourmodel predicts the collective behavior of this clusterand the bistability in external electric fields, resulting in hysteresis ofthe induced dipole fields acting on the Gr NRs. The robusttemperature-independent hysteresis response, characterized by highremnant fields in water-terminated graphene nanoribbons, makesthem promising candidates for the creation of heterostructures withother 2Dmaterials. Improvements using dual gate architectures, high-k dielectrics, and fabricating narrower ribbons could lead to thedevelopment of water-switch memristors. In addition, using ribbonnetworks as electrodes in vertical FET configurations, combined withdipole fields to modulate an atomically thin 2D semiconductor placedover the ribbon network, would enable the creation of reconfigurablememory transistors. These innovations could significantly enhancecore functionalities for future computing in memory and synapticcircuits based on water switches.Methods2D material heterostructuresCrystals of 2D materials were mechanically exfoliated from bulkmaterial and transferred onto a 290–300nm SiO2/Si substrate usingcommercially available Nitto tape. The substrates were chemicallycleaned by sonication in acetone and IPA, aswell as plasma cleaned in amixture of oxygen and argon. Mono-layer and multi-layer thick flakeswere identified via optical contrast and Raman measurements. Mono-layer graphene crystals were transferred onto hBN substrates usingdeterministic dry transfer methods. The transfer was achieved usingmetal-coated SiNx membranes rather than the conventional polymerstamp transfer techniques50. Graphite electrodes and top hBN cappingwere transferred via a PDMS-based deterministic dry transfer method.Organic nanostructuresOrganic molecules were grown by hot-wall epitaxy11. As a source,parahexaphenyl (6 P) was used. The molecules were sublimated at493.15 K under 2·10−6mbar pressure in a vertically mounted semi-closed quartz tube. The walls of the tube were kept at 513.15 K. Thereactor tube is closed by the sample holder that is kept at a tempera-ture of 383.15 K. To ensureminimal temperature variation between thesubsequent growths, the systemwas left to stabilize for 45min after allthree growth temperatures were reached, yielding under ± 1 K varia-tion during the growth. The samples were exposed to the 6 P vapor for4–6min, resulting in an equivalent of 0.8ML (mono-layer) – con-sidering layer-by-layer growth of 6 P on the surrounding SiO2/Si. Anexample of the nanoribbon network (AFM topography map) after theetching step is presented in Supplementary Fig. S2.Reactive ion etchingThe etching parameters were adopted from the reference11. Reactiveion etching (Oxford Instruments PlasmaLab 100) was done in twosteps;first employing SF6plasma (50 sccm,80W, 8-10 s) to remove thetop hBN layer, followed by oxygen plasma (50 sccm, 80W, 4–10 s) toremove the graphene layers and to oxygenate the nanoribbon edges.Alternatively, prolonged SF6 etching (20 s) was used to etch throughthe whole stack, followed by brief oxygen etching to oxygenate gra-phene nanoribbon edges.To ensure that the water molecules can populate the oxygenatedribbon 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 30minutes was applied under these con-ditions for each device in between the RIE step and the electricalmeasurements under vacuum. Further, no significant difference wasobserved between the freshly etched devices (with the 30min holdingtime) and the devices that were tested, even several days after the RIEstep, also stored under ambient conditions.Electrical characterizationElectrical characterization of the flake and Gr NR FETs was doneusing a Keithley 2636A Source-Meter attached to the Instec vacuumprobe station. The samples were contacted with Au-coated Ti elec-trical cantilever microprobes. Variable temperature (77 - 425 K)electrical measurements were performed using a liquid nitrogen flowcryostat with the devices mounted on a silver plate for thermal uni-formity. The temperatures were monitored via Instec’s mK2000temperature controller connected to the probe station with a reso-lution of 0.01 K. To ensure stable and reproducible electrical mea-surements, upon reaching the targeted temperature, the system wasleft to thermalize for 10min prior to the electrical sweeps. Thisensures that during the electrical measurements, the variation of thetemperature is below ± 0.5 K. Cryo (4.2 K) measurements were per-formed using a closed-cycle He cryostat and Lake Shore cryotronicstemperature controller. All electrical measurements carried out inthe Instec probe station (temperature range of 77–673 K) were per-formed at (2.2 ± 0.3)·10−2mbar or as specified in the figure captions.The electrical measurements carried out at 4.2 K were done underthe pressure of (4.0 ± 0.1)·10−6mbar. To ensure stable low tempera-ture, a holding time of 60minutes was applied between reaching thetargeted temperature and the start of the measurements. Thisensured that the fluctuations were below ±0.1 K at 4.2 K. Sampleannealing was done in the Instec probe station under(8.5 ± 0.5)·10−3mbar, for 60minutes and at 520K. The annealingparameters were set to remove the leftovers from the organic masksuccessfully, and not to damage the ribbons or to compromise thegate dielectric. The resistance values are presented as the totaldevice resistance and were measured in a two-terminal configura-tion. The presented resistance values are neither corrected for thelead resistance nor for the contact resistance between the nanor-ibbons and graphite electrodes.Electrical characterization under varied vacuum conditions(Supplementary Fig. S12) was performed by dosing air with a con-trolled relative humidity level of 40% using a mass flow meter fromOmega, with a flow precision of 0.1 sccm (standard cubic centimeters)and a range of 0.1–200 sccm. Vacuum levelsweremaintained to below1% variation during the electrical sweeps. To achieve lower pressurelevels, the flow of the humid air was set against the continuouspumping. To achieve the pressure levels in the range of 102mbar, thechamberwasdisconnected from thepumps, and a small doseof humidair was released to the chamber.AFM measurementsAFM measurements were performed using the Horiba/AIST-NT Ome-gascope AFM system. Nanosensor probes were used (resonant fre-quency ~ 300 kHz, tip radius of 2–4 nm). AFM images were processedin the open-source software Gwyddion v2.56. For AFM image proces-sing, zero-order line filtering was applied, leveling of the base plane,and tip deconvolution were considered, using a model tip with aradius of 4 nm.MD simulationsThe local configuration of H2Omolecules and the binding states at theoxygenated graphene nanoribbon (Gr NR) edges were investigatedArticle https://doi.org/10.1038/s41467-025-65922-6Nature Communications |        (2025) 16:10982 7www.nature.com/naturecommunicationsusing aReaxFF atomisticmodel incorporatedwithmoleculardynamics(MD) simulations51–53. The model was applied to mon- and bi-layer Grnanoribbons and includes the elements of the experimental system.The graphene NR is 20 nm wide in the x-direction, terminated withoxygen, and 34 nm long in the y-direction, with applied periodicboundary conditions. To achieve complete (1:1) edge coverage, wechose 4.7 and 9.4 water molecules per nm length of the edge. To testthe influence of one edge cluster on the other, the NRwidthwas variedbetween 5 nm and 30 nm, considering only a mono-layer Gr nanor-ibbon. We calculated both in-plane and out-of-plane electric fieldsgeneratedby the edges for variouswidths of the grapheneband. Basedon these calculations, we found that for ribbons wider than 5 nm, theelectric fields generated within the cluster become two orders ofmagnitude larger than the average effective field exerted by onecluster on the other.The MD simulations were performed using time steps of 1 fs andan NVT thermostat in LAMMPS, a commonly used distributed classicalMD code51. The simulated water temperature ranged from 250K to450K, depending on the simulation, while the positions of carbonatoms in the ribbonwerefixed tomaintain aplanar structure, similar tothe experimental setup. The interatomic forces within the graphene-water system were derived using the appropriate ReaxFF potential byZhang et al.52,53. In addition, a ReaxFF model by Chenoweth et al.54 wastested, with no qualitative differences observed. Charge equilibration(Qeq) is used in molecular dynamics simulations to calculate the dis-tribution of partial charges on the graphene ribbon edge and watermolecules. This distribution changeswith time tomatch changes in thelocal 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 tem-perature dependence of the hysteresis width, as these field strengthswere estimated to be achievable in the experiment, considering thefringing capacitance effects (see also Supplementary Fig. S16). Toachieve 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 cyclesare shown.Data availabilitySource Data: The complete dataset generated in this study has beendeposited in the Zenodo database56 Source data are provided inthis paper.References1. Garcia, R. Interfacial liquid water on graphite, graphene, and 2Dmaterials. ACS Nano 17, 51–69 (2022).2. Maréchal, Y. The Hydrogen Bond and the Water Molecule: ThePhysics and Chemistry of Water, Aqueous and Bio-Media. (Else-vier, 2006).3. Emmerich, T. et al. Enhanced nanofluidic transport in activatedcarbon nanoconduits. Nat. Mater. 21, 696–702 (2022).4. Mendoza-Sánchez, B. & Gogotsi, Y. Synthesis of two-dimensionalmaterials for capacitive energy storage. Adv. Mater. 28,6104–6135 (2016).5. Artemov, V. et al. Bulk electricity storage in 1-nm water channels.Preprint at https://arxiv.org/abs/2410.11983v3 (2025).6. Cai, J. et al. Wien effect in interfacial water dissociation throughproton-permeable graphene electrodes. Nat. Commun. 13,5776 (2022).7. Dervin, S., Dionysiou, D. D. & Pillai, S. C. 2D nanostructures for waterpurification: graphene and beyond. Nanoscale 8,15115–15131 (2016).8. Surwade, S. P. et al. Water desalination using nanoporous single-layer graphene. Nat. Nanotechnol. 10, 459–464 (2015).9. Chin, H. T. et al. Ferroelectric 2D ice under graphene confinement.Nat. Commun. 12, 6291 (2021).10. Caridad, J. M. et al. A graphene-edge ferroelectric molecularswitch. Nano Lett. 18, 4675–4683 (2018).11. Aslam,M.A. et al. Single-crystalline nanoribbonnetworkfield effecttransistors from arbitrary two-dimensional materials. NPJ 2D Mater.Appl. 6, 76 (2022).12. Wang, W. & Li, Z. Potential barrier of graphene edges. J. Appl. Phys.109, 114308 (2011).13. Kvashnin, D.G., Sorokin, P. B., Brüning, J.W. &Chernozatonskii, L. A.The impact of edges and dopants on thework function of graphenenanostructures: The way to high electronic emission from purecarbon medium. Appl. Phys. Lett. 102, 183112 (2013).14. Jiang, L., Wang, J., Liu, P., Song, W. & He, B. Study of wateradsorption on graphene edges. RSC Adv. 8, 11216–11221 (2018).15. Tan, Y. Z. et al. Atomically precise edge chlorination of nano-graphenes and its application in graphene nanoribbons. Nat.Commun. 4, 2646 (2013).16. Wagner, P. et al. Bandgapengineering via edge-functionalization ofgraphene nanoribbons. J. Phys. Chem. C. 117, 26790–26796 (2013).17. Liu, X., Liu, Y., Chen,W., Li, J. & Liao, L. Ferroelectric memory basedon nanostructures. Nanoscale Res. Lett. 7, 285 (2012).18. Liu, C. et al. Two-dimensional materials for next-generation com-puting technologies. Nat. Nanotechnol. 15, 545–557 (2020).19. Guan, Z. et al. Recent progress in two-dimensional ferroelectricmaterials. Adv. Electron. Mater. 6, 1900818 (2020).20. Li, Y. et al. Enhanced bulk photovoltaic effect in two-dimensionalferroelectric CuInP2S6. Nat. Commun. 12, 5896 (2021).21. Krems,M., Pershin, Y. V. &Di Ventra,M. Ionicmemcapacitive effectsin nanopores. Nano Lett. 10, 2674–2678 (2010).22. Ding, Y. et al. Transition Metal Dichalcogenides HeterostructuresNanoribbons. ACS Mater. Lett. 5, 1781–1786 (2023).23. Canton-Vitoria, R., Hotta, T., Xue, M., Zhang, S. & Kitaura, R.Synthesis and characterization of transition metal dichalcogenidenanoribbons based on a controllable O2 etching. J. Am. Chem. Soc.Au 3, 775–784 (2023).24. Murastov, G. et al. Photoinduced edge-specific nanoparticle dec-oration of two-dimensional tungsten diselenide nanoribbons.Commun. Chem. 6, 166 (2023).25. Ma, Z. et al. Lattice-guided growth of dense arrays of aligned tran-sition metal dichalcogenide nanoribbons with high catalytic reac-tivity. Sci. Adv. 11, eadr8046 (2025).26. Plačkić, A. et al. Electrochemistry at the Edge of a van der WaalsHeterostructure. Small 20, 2306361 (2024).27. Huang, Y. et al. Ultrasensitive quantum capacitance detector at theedge of graphene. Mater. Today 73, 38–46 (2024).28. Kempt, R., Kuc, A., Brumme, T. & Heine, T. Edge Conductivity inPtSe2 Nanostructures. Small Struct. 5, 2300222 (2024).29. Gabbett, C. et al. Understanding how junction resistances impactthe conduction mechanism in nano-networks. Nat. Commun. 15,4517 (2024).30. Matković, A. et al. Light-AssistedChargePropagation inNetworksofOrganic Semiconductor Crystallites on Hexagonal Boron Nitride.Adv. Funct. Mater. 29, 1903816 (2019).31. Kratzer, M. et al. Effects of polymethylmethacrylate-transfer resi-dues on the growth of organic semiconductor molecules on che-mical vapor deposited graphene. Appl. Phys. Lett. 106,103101 (2015).32. Chang, T. H. et al. Directed self-assembly of block copolymer filmson atomically-thin graphene chemical patterns. Sci. Rep. 6,31407 (2016).33. Tavakkoli, K. G. et al. Multilayer block copolymer meshes byorthogonal self-assembly. Nat. Commun. 7, 10518 (2016).Article https://doi.org/10.1038/s41467-025-65922-6Nature Communications |        (2025) 16:10982 8https://arxiv.org/abs/2410.11983v3www.nature.com/naturecommunications34. Tung, C. H. et al. Directed self-assembly of polystyrene-block-polyhedral oligomeric silsesquioxane monolayer by nano-trenchfor nanopatterning. Small 20, 2403581 (2024).35. Kim, J. H. et al. Atomically flat, 2D edge-directed self-assembly ofblock copolymers. Adv. Mater. 35, 2207338 (2023).36. Berashevich, J. &Chakraborty, T. Doping graphene by adsorption ofpolarmolecules at theoxidized zigzagedges.Phys. Rev. BCondens.Matter Mater. Phys. 81, 205431 (2010).37. Lin, X., Ni, J. & Fang, C. Adsorption capacity of H2O, NH3, CO, andNO2 on the pristine graphene. J. Appl. Phys. 113, 034306 (2013).38. Winarto, Yamamoto, E. & Yasuoka, K. Water molecules in a carbonnanotube under an applied electric field at various temperaturesand pressures. Water 9, 473 (2017).39. Offei-Danso, A.,Morzan, U. N., Rodriguez, A., Hassanali, A. & Jelic, A.The collective burst mechanism of angular jumps in liquid water.Nat. Commun. 14, 1345 (2023).40. Wang,H.,Wu, Y., Cong,C., Shang, J. & Yu, T.Hysteresis of electronictransport in graphene transistors. ACS Nano 4, 7221–7228 (2010).41. Aslam, M. A. et al. All van der waals semiconducting PtSe2 fieldeffect transistors with low contact resistance graphite electrodes.Nano Lett. 24, 6529–6537 (2024).42. Singh, A. K. &Gupta, A. K. Reversible control of doping in graphene-on-SiO2 by cooling under gate-voltage. J. Appl. Phys. 122,195305 (2017).43. Park, Y., Baac, H. W., Heo, J. & Yoo, G. Thermally activated trapcharges responsible for hysteresis in multilayer MoS2 field-effecttransistors. Appl. Phys. Lett. 108, 083102 (2016).44. Oostinga, J. B., Heersche, H. B., Liu, X., Morpurgo, A. F. & Vander-sypen, L. M. Gate-induced insulating state in bilayer graphenedevices. Nat. Mater. 7, 151–157 (2008).45. Nam, Y., Ki, D. K., Koshino, M., McCann, E. & Morpurgo, A. F.Interaction-induced insulating state in thick multilayer graphene.2D Mater. 3, 045014 (2016).46. Chen, C. et al. Water-induced bandgap engineering in nanoribbonsof hexagonal boron nitride. Adv. Mater. 35, 2303198 (2023).47. Klapp, S. H. Collective dynamics of dipolar andmultipolar colloids:Frompassive to active systems.Curr. Opin. Colloid Interface Sci. 21,76–85 (2016).48. Klumpp, S. & Lipowsky, R. Cooperative cargo transport by severalmolecular motors. Proc. Natl. Acad. Sci. USA 102,17284–17289 (2005).49. Bagheri, M. H., Loibl, R. T. & Schiffres, S. N. Control of wateradsorption via electrically doped graphene: effect of fermi level onuptake and H2O orientation. Adv. Mater. Interfaces 8,2100445 (2021).50. Wang, W. Polymer-free assembly of ultraclean van der Waals het-erostructures. Nat. Rev. Phys. 4, 504–508 (2022).51. Thompson, A. P. et al. LAMMPS-a flexible simulation tool forparticle-based materials modeling at the atomic, meso, and con-tinuum scales. Comput. Phys. Commun. 271, 108171 (2022).52. Aktulga, H. M., Fogarty, J. C., Pandit, S. A. & Grama, A. Y. Parallelreactive molecular dynamics: Numerical methods and algorithmictechniques. Parallel Comput. 38, 245–259 (2012).53. Zhang, W. & van Duin, A. C. T. Improvement of the ReaxFFDescription for Functionalized Hydrocarbon/Water Weak Interac-tions in the Condensed Phase. J. Phys. Chem. B 122,4083–4092 (2018).54. Chenoweth, K., Van Duin, A. C. & Goddard, W. A. ReaxFF reactiveforce field for molecular dynamics simulations of hydrocarbonoxidation. J. Phys. Chem. A 112, 1040–1053 (2008).55. Nakano, A. Parallel multilevel preconditioned conjugate-gradientapproach to variable-charge molecular dynamics. Comput. Phys.Commun. 104, 59–69 (1997).56. The data accompanying this manuscript is available via: https://zenodo.org/records/16961164 or under https://doi.org/10.5281/zenodo.16961164 (2025).AcknowledgementsThis work is supported by the Austrian Science Fund (FWF) under grantsno. I4323-N36 and Y1298-N, and the TPU Development Program Priority2030. A.M. acknowledges the support from the ERC Starting grantPOL_2D_PHYSICS (101075821). I.S. acknowledges the support of theMinistry of Education, Science and Technological Development of theRepublic of Serbia through the Institute of Physics Belgrade and theEuropean Union through the ULTIMATE-I project partner Senzor Infizdoo, grant ID 101007825. Molecular dynamics calculations were run onthe PARADOX super-computing facility at the Scientific ComputingLaboratory, Center for the Study of Complex Systems of the Institute ofPhysics, Belgrade. K.W. and T.T. acknowledge support from the JSPSKAKENHI (Grant Numbers 20H00354 and 23H02052) and World Pre-mier International Research Center Initiative (WPI), MEXT, Japan. R.G.acknowledges support from the Royal Society, ERC Consolidator grantQTWIST (101001515) and EPSRC grant numbers EP/V007033/1, EP/S030719/1 and EP/V026496/1. A.C. acknowledges support from EPSRCCDTGraphene NOWNANO, grant EP/L01548X/1. RDR acknowledges theAgrochemical Engineering Innovation and Intelligence Base for OasisEcology. Views and opinions expressed are however, those of theauthors only and do not necessarily reflect those of the European Unionor the European Research Council Executive Agency; neither the Eur-opean Union nor the granting authority can be held responsiblefor them.Author contributionsM.A.A. and M.Z.K. prepared the devices, carried out experiments, anddata analysis under the supervision of A.M. I.S. proposed the underlyingmechanism and performed the simulations. G.M. performed the mea-surements for the field remanence. M.A.A., I.S., R.D.R. and A.M. wrotethe manuscript. M.A.A., M.Z.K. and Z.S. carried out etching experimentsunder the supervision of A.L. K.W. and T.T. provided hexagonal boronnitride crystals. R.G. and A.C. provided monolayer graphene on hex-agonal boron nitride. R.D.R.,M.Z.K. andC.T. helped in the internal reviewof the manuscript. A.M. and R.D.R. acquired the main funding for thestudy. All the authors discussed the results and reviewed themanuscript.Competing interestsThe authors declare no competing financial and non-financial interests.Additional informationSupplementary information The online version containssupplementary material available athttps://doi.org/10.1038/s41467-025-65922-6.Correspondence and requests for materials should be addressed toIgor Stanković, Raul D. Rodriguez or Aleksandar Matković.Peer review information Nature Communications thanks Carlo Grazia-netti and the other anonymous reviewer(s) for their contribution to thepeer review of this work. A peer review file is available.Reprints and permissions information is available athttp://www.nature.com/reprintsPublisher’s note Springer Nature remains neutral with regard to jur-isdictional claims in published maps and institutional affiliations.Article https://doi.org/10.1038/s41467-025-65922-6Nature Communications |        (2025) 16:10982 9https://zenodo.org/records/16961164https://zenodo.org/records/16961164https://doi.org/10.5281/zenodo.16961164https://doi.org/10.5281/zenodo.16961164https://doi.org/10.1038/s41467-025-65922-6http://www.nature.com/reprintswww.nature.com/naturecommunicationsOpen Access This article is licensed under a Creative CommonsAttribution 4.0 International License, which permits use, sharing,adaptation, distribution and reproduction in any medium or format, aslong as you give appropriate credit to the original author(s) and thesource, provide a link to the Creative Commons licence, and indicate ifchanges were made. The images or other third party material in thisarticle are included in the article’s Creative Commons licence, unlessindicated otherwise in a credit line to the material. If material is notincluded in the article’s Creative Commons licence and your intendeduse is not permitted by statutory regulation or exceeds the permitteduse, you will need to obtain permission directly from the copyrightholder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.© The Author(s) 2025Article https://doi.org/10.1038/s41467-025-65922-6Nature Communications |        (2025) 16:10982 10http://creativecommons.org/licenses/by/4.0/http://creativecommons.org/licenses/by/4.0/www.nature.com/naturecommunications Ferroelectricity in graphene nanoribbon devices enabled by collective water molecule dynamics Results Observation of the hysteretic response in oxygen-terminated graphene nanoribbon devices Stability of the hysteretic response: excluding charge trap origins Molecular dynamics simulations of the oxygenated and water-terminated graphene nanoribbons Probing the remanence of the water dipole-induced fields in graphene nanoribbons Discussion Methods 2D material heterostructures Organic nanostructures Reactive ion etching Electrical characterization AFM measurements MD simulations Data availability References Acknowledgements Author contributions Competing interests Additional information