The purported square ice in bilayer graphene is a nanoscale, monolayer object

The dominant, opposing thermodynamic forces in bulk water are hydrogen bonding enthalpy and entropy,^1–5 favoring rigid ice superstructures and isolated, disordered gas molecules, respectively. At ambient conditions, the liquid phase is the thermodynamic minimum in free energy because it offers the best compromise between these two competing forces, a fact that is critical for the existence of life on our wet planet.⁶ Indeed, experimental calorimetric measurements show that water molecules in the liquid possess H⁰ = −0.34 eV/H₂O of absolute enthalpy and TS⁰ = 0.21 eV/H₂O of absolute entropy at ambient conditions.⁷ In comparison, hexagonal ice is unstable under the same conditions (ΔG = +0.027 eV/H₂O), since for this “superheated” solid, enthalpic stabilization (ΔH = −0.060 eV/H₂O) cannot overcome entropic losses (TΔS = −0.087 eV/H₂O). At the other extreme, gas molecules have significantly enhanced relative entropy (TΔS = +0.354 eV/H₂O), but even less favorable enthalpy (ΔH = +0.44 eV/H₂O), and thus present another unstable phase at ambient pressure and temperature (ΔG = +0.086 eV/H₂O).

At interfaces, the excess surface free energy underlies the physics of various natural processes, such as cloud formation⁸ and protein folding,⁹ as well as useful industrial systems, including emulsions¹⁰ and water repellant surfaces.¹¹ More commonly, the wettability of solid surfaces is known to arise from nonlinear combinations of liquid/vapor (i.e., surface tension), solid/liquid, and solid/vapor surface energies.¹² The so-called “hydrophobic”¹³ surfaces have positive solid/liquid surface energies and are generally characterized by large (>90°) contact angles and diffuse liquidlike interfacial layers.¹⁴ At the nanoscale, the surface area becomes comparable to volume, and interfacial thermodynamics play a defining role in minimizing the free energy, sometimes stabilizing unusual phases.^15,16 In one recent example, both experimental measurements and theoretical calculations show the existence of stable, nanoscale, one-dimensional (1D) ice phases at ambient conditions in carbon nanotubes.^17,18 It is therefore natural to imagine that similar physics might act to stabilize ice phases upon planar (2D) nano-encapsulation.¹ 

Bilayer graphene, a superstructure comprising two layers of carbon atoms, is an attractive host material for such a phase, since, although the individual sheets are strong (Et = 340 N/m)¹⁹ and flexible,²⁰ the structure is held together by weak van der Waals forces. Recently, high resolution transmission electron microscopy (HRTEM) images of water encapsulated in bilayer graphene¹ have been reported, and it was proposed that this is a new water phase: a reduced dimensional analog of ice VII.²¹ This 2D phase is comprised of a square lattice with completely saturated hydrogen bonding in the monolayer, something not guaranteed for other crystallographic arrangements. The interpretation of this work has since been challenged and has not been subsequently demonstrated in the ensuing years.²² Here, we advance a novel structure descriptor and accelerated molecular dynamics (MD) simulations to reversibly form square ice (SqI) structures of various sizes within bilayer graphene, starting from an initial encapsulated water droplet (EWD). We decompose the free energy of the various ensemble structures into the separate enthalpy, entropy, and differential work terms and show that the balance between these forces conspires to restrict the size of square ice to ∼15 nm in diameter. Thus, square ice is a nanoscale object, which should be spectroscopically discernible based on simulated core-level loss spectra, comprising a single, atomic monolayer that can be immobilized by chemical defects in the graphene sheet.

We performed density functional theory (DFT) calculations at the vdw-DF2²³ level of theory to obtain the square ice unit cell in Fig. S1a of the supplementary material. These results are consistent with previous DFT calculations.^1,24 Periodically infinite, square ice sheets were constructed²⁵ (Fig. S1b) under periodic boundary conditions and subjected to extensive equilibrium molecular dynamics (MD) simulations using empirical force fields. Here, we described the water–water interactions using the fluctuating charge, TIP4P-Fq²⁶ water model, while the carbon-carbon interactions employed the QMFF-Cx²⁷ force field for graphene. The critical carbon–water interactions were described in the same fluctuating charge framework,²⁸ based on fitting to DFT calculations of the orientationally dependent, binding energy curves of an isolated water molecule on graphene (Tables S1 and S2). All simulations were carried out using the large-scale atomic/molecular massively parallel simulator (LAMMPS) MD engine.²⁹ 

We found that our infinite square ice MD simulations were only mechanically stable at room temperature through the application of constraints on the separation (center of mass) between the graphene sheets. Unconstrained simulations invariably resulted in the destruction of the square ice motif and the emergence of an encapsulated water droplet (EWD) on nanosecond time scales. This is consistent with other theoretical approaches that were only able to realize square ice through external confining potentials^24,30 or rigid/immobilized graphene sheets.¹ MD simulations of extended structures comprising unit cells with alternative rhombohedral symmetries³¹ led to similar conclusions.

The existence of such constraints is difficult to rationalize with the original experiments; thus, we adopted an alternative approach, simulating systems comprising various numbers of water molecules embedded inside a 40 × 40 nm² bilayer graphene supercell (corresponding to submonolayer coverage—Table S3). Initial MD equilibration resulted in a single, mixed phase EWD, characterized by a periphery of monolayer square ice (SqI), an interfacial monolayer (in contact with the graphene) and an internal core of disordered water molecules, more clearly evident in larger droplets due to their large aspect ratio. The calculated oxygen–oxygen pair distribution function shows a broad distribution, with a central peak at 2.78 Å, reminiscent of the bulk liquid, though smaller than the 2.85 Å separation in the constrained, extended-phase SqI (Fig. S2). Notably, our equilibrium MD simulations did not result in spontaneous SqI formation on the nanosecond time scale.

In order to affect reversible conversion between the EWD and SqI phases, we performed accelerated MD simulations, employing the Q₂ and Q₅ Steinhardt bond-order parameters,^32,33

which quantify the coordination environment around each water molecule i within a mathematical basis of a vector of m spherical harmonics functions, Y_lm(θ,φ), of neighboring molecules N_b(i). As shown in Fig. S3, this combination of order parameters can unambiguously distinguish square ice from disordered liquid phases. We thus constructed 2D free energy surfaces (FES) by exclusive sampling of the (Q₂,Q₅) subspace^34,35 using the metadynamics scheme.^36,Figure 1 presents the thermodynamics and equilibrated snapshots of the SqI and EWD phases with 1600 water molecules. The (Q₂,Q₅) FES has two deep basins for the EWD and SqI structures, and a significant transition barrier of ∆G^* = 0.567 ± 0.030 eV/H₂O between the two.³⁷ This large barrier would preclude the spontaneous formation of SqI on the nanosecond time scales typical of previous molecular dynamics simulations.

The SqI obtained from accelerated MD sampling of the (Q₂,Q₅) FES contains numerous defects and grain boundaries. Recent work for small systems (<1000 water molecules) has shown that through further equilibration, it is possible to remove many of these defects, through water “flipping” mechanisms.³⁸ Here, we realized defect-free square ice by encapsulating a 20 × 20 supercell inside bilayer graphene and performing constrained dynamics [restricting the motion of the water molecules with a 100 (kcal/mol)/Ų harmonic spring] over 1 ns of MD equilibration. Once the graphene sheet had accommodated the square ice structure, we removed the constraints and evaluated the free energy. We found that this, defect free, square ice structure was lower in free energy by ∼0.014 eV/H₂O than the corresponding, defected motif, similar to the findings of Zhu and co-workers.³⁸ Interestingly, in all cases, the water molecules experience compressive stresses, equivalent to >1 GPa of pressure, due to encapsulation in bilayer graphene, although the stress in the SqI phase is anisotropic and limited to the SqI plane (Fig. S4). Pressures of a similar order of magnitude have been reported experimentally for graphene nanobubbles.³⁹ 

Figures 1(c) and 1(d) and Table S3 show, as expected, that SqI has significantly larger surface area than EWD. While the volume of SqI and EWD is comparable (reflecting the low compressibility of water in condensed phases), the interfacial surface area of SqI necessarily scales linearly with the number of water molecules N, whereas that of the EWD scales as ∼N^2/3 for systems with at least 100 molecules. The EWD size scaling is consistent with continuum mechanical models of blisters⁴⁰ produced by fluids under thin, elastic, adhesive membranes. In this theory, advanced by Hencky and others,^41–44 the blister dimensions (radius and height) are related to the pressure of the confined fluid, the Young’s modulus of the membrane, and the adhesive energy per unit area of the membrane to its substrate. We show in the supplementary material how this continuum theory correctly predicts the size and shape of nanoscale EWDs.

Further insights into the SqI thermodynamics were obtained by means of free energy calculations, employing the Two-Phase Thermodynamics (2PT) method.^45–48 This approach allows us to separately determine the thermodynamics of the water molecules and graphene sheets while also accounting for the separate diffusive, configurational, and vibrational components of the entropy. These calculations reveal, consistent with our previous metadynamics simulations, that for 1600 water molecules encapsulated in bilayer graphene, SqI is more thermodynamically stable than EWD: ∆G = G⁰_SqI − G⁰_EWD = −0.019 ± 0.003 eV/H₂O. Here, gains in the system enthalpy (∆H = −0.37 ± 0.003 eV/H₂O) overcome reduction in the total entropy (T∆S = −0.35 ± 0.005 eV/H₂O). This result stands in contrast to the mechanically unstable infinite SqI sheets and implies a size dependence to the free energy difference, which we elucidate by means of a model Hamiltonian,

Here, the relative energy of the SqI and EWD phases results from (i) reduced graphene stretching to accommodate a flat 2D flake rather than a 3D droplet (Fig. S5), (ii) loss of favorable graphene interlayer binding due to the increased projected area (Fig. S6), and (iii) an increase in favorable graphene-water interfacial binding (Fig. S7), for the same reason. The last two terms in Eq. (2) reflect differences in (iv) the water enthalpy (Fig. S8a) and (v) the water entropy (Fig. S8b) due to obvious changes in the degree of hydrogen bonding and disorder between the solid and liquid, respectively.

Analysis of MD trajectories of the various motifs constrained to their respective region of the (Q₂,Q₅) FES confirms our hypothesis, revealing that terms (i), (iii), and (iv) are negative, favoring SqI formation, while terms (ii) and (v) are destabilizing forces (Fig. S9). The surface area scaling, noted previously, coupled with the coefficients from the various energy terms in Eq. (2), thus leading to an equation for the relative energies as a function of N and temperature T (Table S4),

Equation (3) is plotted in Fig. 2(a) and found to predict the relative stability of SqI with increasing N at 298 K consistent with the free energy MD simulations. Clearly, the functional form of Eq. (3) implies correctly that the extended 2D SqI phase is ultimately unstable with respect to the 3D EWD phase, underscoring the many necessary constraints that have been applied in previous studies to stabilize SqI. A natural consequence of Eq. (3) is to place restrictions on the maximum size of SqI flakes that are thermodynamically stable compared to EWD (i.e., ΔG ≤ 0). This we find has two main components: considering only the relative energies of the water molecules, enthalpic stabilization dominates at smaller sizes but eventually is overcome by entropic losses so that SqI would be thermodynamically stable at 298 K only for N < 8000 (Fig. S8d). However, the destabilizing loss of graphene–graphene interlayer binding in SqI reduces the overall free energy of the system, dramatically, limiting the maximum size to N ∼ 2300 molecules (∼15 nm diameter). We note that SqI would be entirely unstable without a significantly stabilizing graphene-water interaction, assisted to some degree by fractional electronic charge transfer. Similarly, stiffer materials, or materials with increased sheet–sheet binding would disfavor SqI formation. Therefore, in many ways, bilayer graphene possesses the optimal mechanical properties to facilitate large SqI flakes at ambient temperature.

Using Eq. (3), we present the T,N phase diagram in Fig. 2(b) to establish the zones of stability and theoretical maximum size of SqI at various temperatures T. Figure 2(c) predicts the melting temperatures of SqI flakes, which we find can be as high as 600 K for an 800 molecule system. Of course, experimentally realized square ice may coexist with various atomic site defects in graphene, which may further stabilize SqI flakes, or they may be trapped out of equilibrium for kinetic reasons. These favorable factors could lead to flakes larger than predicted here. Finally, we acknowledge that the properties of 2D water/ice in MD simulations are sensitive to the force field used. However, we verified that our results are robust by comparing the thermodynamics of 1600 water molecule SqI and EWD phases with various popular water force fields. In all cases, the SqI phase was found to have a lower free energy (Table S5).

The overall nanoscale morphology of SqI was determined by simulated oxygen K-edge EELS spectroscopy, using constrained orbital-occupancy DFT^49,50 within the excited electron and core-hole (XCH) formalism.⁵¹ As shown in Fig. 3, monolayer SqI flakes have a spectroscopic signature that differs significantly from disordered liquid phases, such as EWD. EELS measurements based on high kinetic energy electrons are naturally dominated by the polarization in the direction of propagation, i.e., perpendicular to the plane of SqI (or equivalently of the graphene bilayer, presumably resting on a TEM grid or sample stage). In the isotropic, bulk-phase liquid, there is no polarization dependence in the EELS spectrum, which exhibits characteristic peaks at 535 eV (due to a population of broken H-bonds; the pre-edge) and 538 eV (the main-edge).⁵² Similarly, the simulated EELS spectrum of EWD does not present strong polarization dependence [Fig. 3(a)], but instead is characterized by a single, red-shifted main-edge at 536 eV, due to increased water–water distances and a suppressed pre-edge. In contrast, the simulated EELS of SqI, with polarization normal to the plane, presents as a single main-edge feature near 541 eV [Fig. 3(b)], in excellent agreement with previous experimental measurement,¹ while in-plane polarization presents a prominent pre-edge peak near 532 eV (Fig. S10). The pre-edge transitions are to antibonding electronic excited states aligned along the OH bonds, perpendicular in SqI to the polarization of the incoming electron beam and hence absent in the measured EELS. Experimentally, we propose that this feature may be seen by tilting the sample stage significantly.

One interesting aspect of the MD simulations of SqI nanostructures is that they are randomly self-diffusive as a structural unit at a rate of ∼2.5 × 10⁻⁶ cm/s² (∼1/10 that of bulk water) due to the center of mass displacement at finite temperature. This motion would be inconsistent with the published HRTEM images,¹ which require no net diffusion within the seconds of data collection required for clear contrast. Thus, we considered “perfect” SqI flakes (i.e., without grain boundaries), which we immobilized computationally through the inclusion of oxygen-containing chemical defects on the interior graphene surfaces to act as hydrogen bonding anchor points. We found that one defect was sufficient to pin the nanoflake, while two constrained the azimuthal angle. There is, in fact, evidence for the presence of oxygen defects in previous experiments. Oxygen K-edge EELS shows the clear signature of graphene oxide motifs⁵ comprising hydroxyls, epoxides, carboxylates, etc., even in regions that did not include visible square ice [Fig. 3(c)]. The chemical specificity of EELS measurements affords definitive peak assignments to graphene oxide, as opposed to other inorganic oxides such as SiO₂.²² 

All of our MD simulations resulted in monolayer SqI, not multilayer AA-stacked structures as have been proposed based on the interpretation of contrast variations in HRTEM images.¹ Ab initio potential energy scans (Fig. S11) on a two layer square ice unit cell revealed the AB-stacked structure as the most energetically stable by 8.5 meV/H₂O, consistent with previous results.³⁰ While various multilayer, nanoscale encapsulated, ice motifs have been realized, such structures do not appear to be thermodynamically accessible when encapsulated in bilayer graphene at room temperature. Moreover, as noted previously, these multilayer ice structures would present out-of-plane H-bonds, which would result in a prominent pre-edge feature in the EELS spectrum, which would be inconsistent with existing measurements. As an alternative explanation of the image contrast in the HRTEM experiments, we performed time-averaged variance analysis on the simulated oxygen atom positions of an immobilized SqI flake. This revealed that interior water molecules have a significantly reduced variance (i.e., they move less) than perimeter waters and we present the sampled distribution as a possible source of image contrast in Fig. 4.

In conclusion, we elucidate the thermodynamic origins of nanoscale (up to 15 nm diameter) 2D square ice flakes as competitive with encapsulated liquid water, driven predominantly by the favorable graphene-water interaction but ultimately destabilized at larger sizes due to entropic losses in the water molecules. We propose that additional, orientation dependent core-level spectroscopy measurements should be conducted to verify the presence of SqI within bilayer graphene. One potential technological application could be in molecular electronics as a nanoscale switch: the rigid H-bonding in SqI should make for an excellent 2D proton conductor, as suggested in the 1D case.^53,54 Based on the stability analysis presented here, phase transitions to the relatively nonconducting EWD phase could be affected either by increased humidity and temperature or by destabilizing the water-carbon interaction through electron doping (charging) of the semimetallic graphene sheets. It would also be interesting to explore the nature of tethering flakes to chemical defects and the possibility of patterning graphene substrates with water, making wet and dry nanoscale regions at will. Future work will utilize kinetic Monte Carlo simulations to explore the role of various chemical defects in possibly stabilizing larger SqI flakes than found here.

See supplementary material for detailed computational methods, graphene blister analysis, Tables S1–S5 and Figs. S1–S11, and unit cells (CIF format) of extended SqI and nanoscale SqI and EWD (1600 water molecules).

Work by T.A.P., C.P.S., and D.P. at The Molecular Foundry was supported by the Office of Science, Office of Basic Energy Sciences, of the U.S. Department of Energy under Contract No. DE-AC02-05CH11231. T.A.P. was supported by the Assistant Secretary for Energy Efficiency and Renewable Energy, Office of Vehicle Technologies of the U.S. Department of Energy under Contract No. DE-AC02-05CH11231 under the Batteries for Advanced Transportation Technologies program. C.P.S. was supported by EERE Bridge Project No. 25860. K.V.L. acknowledges support through the National Nuclear Security Administration under the Stewardship Science Academic Alliances program through DOE Cooperative Agreement No. DE-NA0001982. K.V.L. would also like to thank the National Supercomputing Center for Energy and the Environment (NSCEE) for computational resources used to explore the DFT energy landscape for stacking square ice. The EELS calculations were performed at the National Energy Research Scientific Computing Center (NERSC), a DOE Office of Science User Facility, supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC02-05CH11231. All remaining molecular dynamics and statistical mechanics calculations were performed on Molecular Foundry computer clusters, Nano and Vulcan, managed by the High Performance Computing Services Group at the Lawrence Berkeley National Laboratory.