Ferroelectric 2D ice under graphene confinement

We here report on the direct observation of ferroelectric properties of water ice in its 2D phase. Upon nanoelectromechanical confinement between two graphene layers, water forms a 2D ice phase at room temperature that exhibits a strong and permanent dipole which depends on the previously applied field, representing clear evidence for ferroelectric ordering. Characterization of this permanent polarization with respect to varying water partial pressure and temperature reveals the importance of forming a monolayer of 2D ice for ferroelectric ordering which agrees with ab-initio and molecular dynamics simulations conducted. The observed robust ferroelectric properties of 2D ice enable novel nanoelectromechanical devices that exhibit memristive properties. A unique bipolar mechanical switching behavior is observed where previous charging history controls the transition voltage between low-resistance and high-resistance state. This advance enables the realization of rugged, non-volatile, mechanical memory exhibiting switching ratios of 106, 4 bit storage capabilities and no degradation after 10,000 switching cycles.

*Corresponding author. Email:yphsieh@gate.sinica.edu.tw Details on finite element simulation Simulation was conducted in a three-dimensional geometry that is based on microscope analysis (Figure 1, Supplementary figure S1). We simulate a structure that contains a 1µm opening in a 40nm dielectric separation layer. The walls of the opening were angled at 140° to account for shadowing effects during evaporation (Supplementary figure S3). Graphene membranes were added to both bottom and top of the separator and were represented by a 1nm thick layer with an elastic modulus of 0.96TPa following previous reports. 1 The graphene's lateral movement was constrained at the edge of the opening. Increasing voltages were applied to both graphene layers until contact was observed.
Upon contact, the strain distribution was extracted from the cutting line through the center of the orifice and convoluted with a gaussian function of FWHM of 1um to provide correspondence with the Raman measurements under a finite laser beam spot 2 .

High frequency resonator measurements
Supplementary figure S4 shows the experimental set-up for radio-frequency read-out and operation. We use three terminals to operate our device. Source and drain electrode where attached to a suspended graphene layer that covers a centimeter-scale array of openings. Highly doped p-type silicon was used as a global back-gate in order to simultaneously operate each device on the chip. A lock-in approach was employed to measure the distance dependent gating effect on graphene (Supplementary figure S4) following previous reports 3 .
Derivation of voltage dependence for PUVD We identify the transport condition between the graphene layers by investigating the currentvoltage behavior. We observe that conduction proceeds by direct tunneling up to 8V before transitioning to Fowler-Nordheim tunneling (Supplementary figure S6 (a)). 4 . This surprisingly high value suggests that a large series resistance exists between the graphene/graphene junction and the electrodes. We estimate the voltage at the junction by considering the tunneling resistance for graphene/hBN/graphene junctions 5 and the directly measured graphene sheet resistance (980Ω/ ). Under these conditions, 90% of the external voltage drops on the graphene leads and the internal electric field is below 3V/nm throughout the experiments. This field strength is within the range, where ferroelectric polarization becomes stable but below the critical fields that lead to dissociation of water and electric breakdown of the dielectric (Supplementary figure  S6(b)). 6 7 . In the direct tunneling regime, the tunneling current in a ferroelectric junction can be described by In the forward bias regime this formula can be approximated by 9 = ( exp(− √ ) − ( + ) exp (− √( + ))) Fits to this formula reveal the tunneling prefactor JT, the tunneling constant A and the barrier height .
The change in barrier height can be related to the ferroelectric polarization of ice through 10 Δ = where c = 2.44 in the units of volts (V) per debye (D), pz is the water molecule's dipole moment, and f is concentration of aligned dipoles. The tunneling prefactor can be furthermore employed to extract the tunneling gap size assuming a unity correction factor 11 = /2 ℎ( ) 2

RF measurements
We employ a vector-network analyzer to extract the complex reflection and transmission coefficients, S11 and S21 for drum-type switch devices (Figure 1(b)). From these measurements we can calculate the dielectric constant following previous reports 12 and the real part of this parameter is shown in Figure 2(e). The frequency dependence of the dielectric constant in a single-oscillator picture is given by the Debye function 13 * ( ) = ′ ( ) − ′′ ( ) = ∞ + − ∞ 1 + 0 and we employ its real part to extract the response time through fitting to

Details on ab-initio simulation
The density functional theory (DFT) calculations were performed using the Vienna ab-initio simulation package (VASP) [14][15][16] . The Perdew-Burke-Ernzerhof (PBE) exchange-correlation functional 17 was used to perform molecular dynamics and structural optimization. For spot checks of energy ordering of different structures the hybrid PBE0 approximation was also used 18 . The D3 dispersion correction of Grimme 19 using the damping function proposed by Becke and co-workers 20 was employed in all the calculations. The plane-wave basis-set cut-off was set to ~400eV and the dipole of the structure in direction perpendicular to graphene was corrected both in energy and potential 21 .
To identify potential ice phases under graphene confinement, we have surveyed a large number of possible geometries in different unit cells, different net dipole orientations (parallel vs perpendicular), and different densities using an automated random structure search at forcefield level. Several low-energy structures for each setting were then selected and reoptimized at the PBE-D3 level. The thus obtained low energy ice phases agree with previous literature reports and are summarized in Suppl. Table 1. We have tested the stability of the zero-dipole "Cairo" and "hexagonal" phases that Chen et al. identified as the most stable for zero or low pressures 22 . The unit cell for Cairo structure was a=29.652 Å, b=9.783 Å, with 112 carbon atoms in each layer and 36 water molecules in the ice layer. For hexagonal ice of Chen at el., the unit cell had parameters a= 16.944Å, b= 19.565 Å with 128 carbon atoms and 32/ 64 water molecules for single/bilayer ice, respectively. The same unit cell was used for the HD-HMI structure. 23 The structures identified by Chen have no dipole. Initial structures with a dipole were obtained by rotating OH bonds that point to -z to +z direction (see Supplementary figure S7). Subsequent optimization was performed using PBE-D3. The ice XI structure and HD-fRMI 24 used a=16.944 Å, b=14.674 Å, =30 Å. We also tested a structure with dipoles aligned and perpendicular to graphene and no hydrogen bonds (last line of T1, Fig S7a). This structure is not stable without an external electric field due to the absence of hydrogen bonds between water molecules, which is consistent with results obtained for single water adsorption on graphene 24 Upon optimization, it transforms to the ferroelectric ice XI structure. This "all-dipole-aligned" structure as well as the ice XI structure and HD-fRMI 23 used a=16.944 Å, b=14.674 Å c=30 Å. We furthermore optimize the dipole-oriented structures under different applied electric fields. In the zero-field case, we observe that the modification of the total energy by a formed dipole, either in-plane or out-of-plane, was negligible (see Suppl. Table 1). This observation agrees with previous simulations on the energetic degeneracy of ice Ih and ice XI phases and originates from the small energy change of rotating hydrogen bonds. Based on these results, the formation of dipole-oriented ice does not occur spontaneously but requires application of an electric field. Our simulations show that ferroelectric ice becomes favorable at fields in excess of 2 V/nm which agrees with simulations by Quian et al. 25 Once formed, the ferroelectric phase is remarkably stable and an energy barrier of 0.2eV per water molecule was calculated for the reorientation of a dipole from the ferroelectric phase. The barrier was obtained using the nudged elastic band method 26 with 5 images between the initial and final structures. As with the other calculations, the PBE-D3 functional was used to obtain the energies. The barrier was obtained in the large simulation cell containing 32 water molecules.
Bilayer ice was built by stacking the obtained net-dipole structure. To form the zero-dipole structure, we again manually rotated the OH bonds to form a structure where no hydrogen interacts with the graphene. All the structures were optimized to find the energy minimum.

Details of Molecular Dynamics Simulation
MD simulations were done in VASP on an NVT ensemble using Andersen thermostat with collision probability of 0.1 per atom and time step. Deuterium mass was used for hydrogens and the time step was 1.0 fs with at least 17 ps run times. The distance between layers is approx. 7.2 Å for the single layer ice and 11.4 Å for the two-layer ice. The simulation cell is cubic with a=16.944 Å and b=14.674 Å, c=30 Å to form a vacuum layer. The graphene unit cell size was chosen to be close to the unit cell of ice and contained 96 carbon atoms per layer and 32 water molecules per layer. Dipole corrections (in energy and potential) along z were used. The planewave basis-set cut-off was 400 eV.

Estimating the numbers of distinct storage levels
In the proposed multilevel storage scheme, the input signal is represented by the priming voltage which encodes a resistance value that can be read out at low voltages. Thus the number of levels (n) is determined by the precision of the measured conductivity Δ = ( − )/ In our case, the precision is limited by the deviation of the measured value from the predicted one. Using a power-dependence to fit the relation between the priming voltage and resistance observed in Figure 3(c) we obtain Δ = | ( ) − 3 × 10 9 −4 | The average Δ in Figure 3(c) was 8.8 Ω and the resistance range for a priming voltage range of 20V was 130 Ω which yields 15 distinct levels of storage. Device optimization could increase this value significantly.

Estimation of stiction force
To calculate the restoring force upon deflection, we employ an analytical formula 27 : where 0 2 , 2 are the 2D pretension and Young's Modulus of the membrane, respectively, is the membrane's Poisson ratio, a the hole diameter, the deflection. Using reported values for graphene 27 we obtain a force of 4nN for a deflection of 20nm. This force can be related to the van-der-Waals force 28 at the graphene/graphene interface at ambient conditions and we find that a circular contact with 8nm radius would provide sufficient stiction to maintain contact between the two membranes.
Evaluation of mechanical switching performance An important parameter to evaluate the power efficiency of transistors and mechanical switches is the on-off ratio between high resistance and low-resistance states 29 . We observe a memristive effect on the on-off ratio which can be tuned through the prior polarization (Supplementary figure S8) and extract a maximum ratio of 10 6 .