Introduction

Twisted bilayer materials (TBMs) have become a hot research frontier in material science and condensed matter physics. TBMs with moiré patterns at specific twist magic angles often demonstrate novel and unique electronic properties. For instance, twisted bilayer graphene at a 1.1° magic angle features flat electronic band structures, while its enhanced electron correlation effects can result in novel physical phenomena such as superconductivity and Mott insulator behavior1,2. This ‘magic-angle’ effect underscores the profound impact of a small twist angle on TBM’s electronic properties, highlighting the importance of precise angular control between stacked monolayers. To date, research on TBMs has predominantly focused on two-dimensional (2D) van der Waals (vdW) materials, where weak interlayer interactions appear to favor the formation of twisted structures3,4,5,6,7,8. In contrast, the exploration of TBMs with stronger non-vdW interlayer interactions remains limited9, primarily due to the challenges to achieving requested twisting angles between two monolayers and demanded moiré patterns.

Water and ice are the most common non-vdW matter on earth and hold vital importance from many perspectives10. The intrinsic hydrogen-bonding interactions endow water and ice with remarkable adaptability to environmental conditions, leading to rich phase behaviors in both bulk and low-dimensional forms. Particularly, when water is confined within nanocapillaries with a width of a few angstroms (Å), a variety of 2D ice structures can form above the freezing point at different lateral pressures, exhibiting physical behaviors that differ markedly from their bulk counterparts11,12,13,14,15.

The rich phases of 2D ice highlight the extraordinary adaptability of hydrogen-bonding networks under nanoconfinement and lateral pressure. Notably, some 2D ices have been found to deviate from the Bernal–Fowler ice rules, such as the “coffin” bilayer ice16, zigzag monolayer ice13 and bilayer ice VII13. Indeed, 2D ice exhibited richer phase structures compared to other covalent 2D materials11,17,18,19. To date, some previously predicted 2D ices have been experimentally observed, enabled by advancements in precise control and measurement techniques17,20,21,22,23. While previous studies have explored bilayer ice with either AA or AB stacking patterns24,25,26,27,28,29, the possibility of forming twisted bilayer ice (TBI) with moiré patterns remains unexplored. With advances in nanofabrication techniques, such as vdW layered material assembly30, achieving precise sub-nanometer confinement is increasingly feasible, making the synthesis of a TBI experimentally accessible. Revealing this possibility would significantly advance our understanding of hydrogen bonding in low-dimensional water systems while finding new physical phenomena in TBIs.

In this study, we report the first simulation evidence of TBI in model nanocapillaries through both classical and first-principles molecular dynamics (MD) simulations. The most prevalent twist angles identified for TBIs are 21.8° and 27.8°, with a smaller fraction of TBI having a twist angle of ~15.7°. Our results show that TBI formation is highly sensitive to the confinement condition and can occur within a narrow range of pressures and nanoconfinement widths. These simulation results offer a foundation for future experimental confirmation. Importantly, the evidence of spontaneous TBI formation fundamentally distinguishes it from the formation of twisted bilayer van der Waals materials, which typically requires manual manipulation to achieve desired twist angles. Additionally, detailed analysis shows unique hydrogen-bonding interactions between the two ice monolayers, which can lead to an orderly distribution of flower-like local structures with a mixture of AA and AB stacking patterns of oxygen atoms and a link-like AB-stacking local structure connecting neighboring flower-like local structures. These local stacking patterns appear to be critical for the formation of the TBI with a global moiré pattern, and they are also remarkably robust during the growth of bilayer ice while maintaining the specific twist angle. Finally, the Gibbs free energy analysis is provided to elucidate the structural stability of TBIs, offering further insights into the molecular mechanisms underlying the TBI formation and stability.

Results

We performed a series of classical MD simulations to investigate the spontaneous freezing of TBI from liquid water under nanoconfinement across a range of widths (h) from 7.5 Å to 8.5 Å, where water molecules encapsulated between two parallel graphene sheets (Supplementary Fig. 1). This range of confinement width h was selected based on previous extensive MD simulations12,13 on bilayer water/ice and was also to span configurations from AB-stacked bilayer ice (~ 7.0 Å) to AA-stacked bilayer ice (~9.0 Å). A phase diagram was constructed to map the stability of TBI phases, with h = 8.0 Å being identified as optimal for stabilizing the reported phase. Here, representative results for a fixed confinement width of h = 8.0 Å are presented as follows. The lateral pressure acting on the confined liquid/ice can range from 0 GPa to several GPa, driven by nanocapillary and hydration pressures, as reported in previous studies17,22. To simulate diverse lateral pressures (up to 4.5 GPa), we systematically varied the number of confined water molecules, corresponding to area densities ranging from 22.45 to 30.87 nm⁻². Lateral pressure calculations, validated against established methods31 (Supplementary Fig. 2), confirmed the applied stress regime.

At room temperature and confinement width of h = 8.0 Å, a 2D bilayer ice with a twisting angle was observed when the lateral pressure exceeded 3.0 GPa, leading to the formation of several intriguing TBI phases with moiré patterns. Figure 1A, D illustrate the two dominant TBI phases with commensurate twist angles of θ = 21.8° and θ = 27.8°, respectively, along with their corresponding moiré potential distributions. The formation dynamics can be seen in Supplementary Movies 1 and 2. Here, the twist angle θ is defined as the smallest angle between the principal axes of the oxygen lattice in the upper and lower ice layers. Other than the TBI with commensurate twist angle of 21.8° and 27.8° with the area density of 30.8 nm-2 and 29.4 nm-2, respectively, a potentially metastable TBI with θ = 15.7° was observed (Supplementary Fig. 3).

Fig. 1: The spontaneous formation of TBI in graphene nanocapillaries.
figure 1

A Top view of the TBI at θ = 21.8° (upper left panel); dimensionless moiré potential contour plot with unit cell highlighted by black dashed lines. (upper right panel). The dimensionless moiré potentials were generated by superimposing the cosine-modulated interlayer interactions of two rotationally misaligned ice lattices; and top views of the upper (lower left panel) and lower (lower right panel) ice layers. Blue and red spheres represent oxygen atoms in the upper and lower layers, respectively. Hydrogen atoms are omitted for clarity. B Potential energy versus time in the MD simulation (blue curve) and full-width at half-maximum (FWHM) of the oxygen arrangement peak in the 21.8° TBI (red curve). Inset: Snapshot of the system from the MD simulation. C Spatial distribution of oxygen atoms in the upper (red) and lower (blue) layers in the 21.8° TBI. D Top view of the TBI at θ = 27.8° (upper left panel); dimensionless moiré potential contour plot with unit cell highlighted by black dashed lines. (upper right panel); and top views of the upper (lower left panel) and lower (lower right panel) layers. Blue and red spheres represent oxygen atoms in the upper and lower layers, respectively. Hydrogen atoms are omitted for clarity. E Potential energy versus time in the MD simulation (blue curve) and FWHM of the oxygen arrangement peak in the 27.8° TBI (red curve). Inset: Snapshot of the system from the MD simulation. F Spatial distribution of oxygen atoms in the upper (red) and lower (blue) layers in the 27.8° TBI.

Full size image

The monolayer in both TBIs at θ = 21.8° and 27.8° features unique flower-like patterns, characterized by two central AA-stacked oxygen atoms that are surrounded by twelve AB-stacked oxygen atoms, as shown in Fig. 1A, D. This flower-like pattern corresponds closely to the characteristic tiling elements of dodecagonal quasicrystals and resembles the local structure of the AB-II phase observed under high lateral pressure within diamond confinement, as previously reported by Goswami and Singh32. However, in their work, the 2D ice did not exhibit long-range transitional order, and therefore, no moiré pattern was formed. In contrast to the 21.8° TBI which only exhibited regularly distributed flower-like local structure, the TBI at θ = 27.8° displayed a more complex superlattice with each flower-like structure surrounded by 12 link-like local structures. The link-like local structure was formed by two oxygen atoms from the upper and lower layers, and stacked in the AB manner. As will be discussed later, these distinct moiré patterns are closely related to the structural adaptability of hydrogen bonding interactions.

To elucidate the growth dynamics of both TBI phases, we analyzed the time-dependent potential energy evolution profiles (Fig. 1B, E) and the spatial distribution of oxygen atoms. In addition, the crystallization process of TBI was quantified using the full width at half maximum (FWHM) of the oxygen atom distribution (Supplementary Fig. 4). A smaller value of the FWHM of the distribution indicates a higher degree of crystallinity. Initially, at 0 ps, a liquid phase with higher energy and disordered oxygen distribution was observed, as shown in Fig. 1B, E. At 25 ps, several flower-like structures were observed showing the nucleation of TBI. At 50 ps, the oxygen atoms in each layer started to exhibit six-fold symmetry, and the growth of TBI was evident through the increasing number of flower-like structures. Lastly, at 10 ns, a nearly perfect TBI formed with lower and converged potential energy, and its FWHM values indicated a perfect twisting at θ = 21.8°. Similarly, for the TBI at θ = 27.8°, the spontaneous liquid-to-solid phase transition occurred on the picosecond timescale. Flower-like structures emerged by 300 ps, and over time, these structures were interconnected by link-like local structures, forming a distinct moiré pattern from TBI at θ = 21.8°. These results highlight the critical role of flower-like structures in the nucleation and growth of TBI, regardless of the twist angle.

To verify the thermodynamic stability of TBIs at first-principles accuracy, we performed machine learning potential (MLP)-based MD simulations (see Methods). The MLP, trained using the vdW-DF2 functional and benchmarked against quantum Monte Carlo simulations33, exhibits the best performance for describing water-water and water-carbon interactions at the GGA-level DFT accuracy. To mitigate the finite-size effects, MLP-based MD simulations used slab areas ranging from 30.5 Å × 26.4 Å to 49.0 Å × 36.4 Å (266–588 water molecules), exceeding the threshold for neglecting the periodicity effects11. The Gibbs free energy of TBIs at arbitrary twist angle was computed (as shown in Fig. 2), where the enthalpy (H) was computed from the average potential energy and pressure-volume work in the constant-temperature-volume (NVT) ensemble, while the entropy (S) was computed using the two-phase thermodynamics method. The TBIs with twist angles of 21.8° and 27.8° exhibited the lowest Gibbs free energies, in excellent agreement with classical MD simulations in which the spontaneous formation of TBIs was observed. However, the small free-energy differences suggest that the formation of structurally perfect TBIs can be highly sensitive to environment conditions such as lateral pressure and confinement width. This sensitivity makes the spontaneous formation of TBI phases challenging if the finite-size effect was not considered and small supercells were used in MD simulations. These findings underscore the need for precise control of experimental conditions to achieve and stabilize TBI phases in the laboratory.

Fig. 2: Computed relative Gibbs free energy of TBI at various twist angles from MLP-based MD simulations.
figure 2

The Gibbs free energy for TBI at θ = 21.8° was subtracted for clarity. Top views of the structures at the corresponding twist angles are displayed in the insets. The error bars represent the standard deviation of energy fluctuations throughout the MD simulations.

Full size image

To further elucidate the origin of thermodynamic stability in TBIs, we conducted a detailed analysis of hydrogen bonds in TBIs. Hydrogen bonds were identified based on the following criteria: (1) the O···H distance (Oi···Hj) is shorter than 2.5 Å, and (2) the O-H···O angle (Oi–Hi···Oj) is greater than 90° (where i and j denote different water molecules)34. We found that TBIs at arbitrary twist angles exhibited stronger interlayer hydrogen bonding, compared to intralayer hydrogen bonding (Supplementary Fig. 5). This indicates that the hydrogen-bonding interactions for maintaining interlayer twist are stronger than intralayer hydrogen-bonding interaction within a monolayer. We therefore further examined the characters of interlayer hydrogen bonds. Take the flowers in the 21.8° TBI as an example, no hydrogen bonds were found between flowers and within each flower, one interlayer hydrogen bond was formed within each pair of water molecules in the opposing layer, resulting in a total of seven intra-flower interlayer hydrogen bonds. Moreover, the interlayer hydrogen bonds in the outer rim of the flower tilt at almost the same angle (bottom panel of Fig. 3A). This tilting feature is essential to the twist stacking, and, once formed, it remains unchanged after growth.

Fig. 3: Interlayer hydrogen bonding networks in TBI at 21.8° and 27.8°.
figure 3

A A top view of the flower-like structure in TBI at θ = 21.8° (top panel) and the corresponding side view of interlayer hydrogen bonding networks (bottom panel). B A top view of the flower-like structure connected via link-like structures in TBI at θ = 27.8° (top panel) and the corresponding side view of interlayer hydrogen bonding networks (bottom panel). C Probability distribution of the interlayer O-H···O angle in the flower-like structure of TBI at θ = 21.8° and 27.8°, as well as in the link structure of TBI at θ = 27.8°. D Probability distribution of the overall interlayer O-H···O angle in TBI at θ = 21.8° and 27.8°.

Full size image

For the 27.8° TBI, each flower-like local structure—comprising the seven flowers shown in Fig. 3B—is linked to six adjacent flower-like structures through AB-stacked oxygen atoms, herein termed the ‘link’. Owing to the slightly larger twist angle, the tilting angles of the interlayer hydrogen bonding within each flower were larger, thereby in some of the flowers, locally reorganized hydrogen bonds were observed, such as the one shown in the left panel of Fig. 3B. It is seen one water molecule (marked in yellow) in the lower layer forms two hydrogen bonds with two water molecules in the upper layer, leaving another water molecule (marked in purple) in the lower layer not participating in the interlayer hydrogen bonding. The paired bonding mechanism in the link is similar to that of the flower region in the 21.8° TBI, with a total of 24 hydrogen bonds counted.

Based on the probability distribution of the O-H···O angle shown in Fig. 3C, the strength of interlayer hydrogen bonds follows this order: 27.8° link >21.8° flower >27.8° flower. Consequently, as illustrated in Fig. 3D, the interlayer hydrogen bond strength in the 21.8° TBI, which consists solely of flower-like structures, is approximately equal to that in the 27.8° TBI, including both flower-like and link structures. We anticipate that the similar and strong interlayer hydrogen bond interactions are the dominant factors that make 21.8° and 27.8° twist angles stand out from the others.

Interestingly, the identified twist angles of 21.8° and 27.8° have also been viewed as large commensurate angles in vdW twisted bilayer systems6. To investigate the origin of the thermodynamic stability of TBIs at θ = 21.8° and 27.8°, compared to smaller twist angles, we modeled artificially constructed TBIs with twist angles ranging from 15° to 2° and conducted a series of MLP-based MD simulations at first-principles accuracy. As shown in Supplementary Fig. 6, for the TBI at θ = 2°, we observed a polycrystalline structure with two distinct regions: one exhibiting the AB phase and the other AA phase. Such a configuration is difficult to achieve during the ripening process, where one phase could transform into the other.

Furthermore, we computed the standard deviation of the oxygen distribution along the z direction in the TBI (Supplementary Fig. 6), revealing an increase in roughness of the water-wall interactions as the twist angle decreases, while the long-range order weakens. This roughness renders the energy unfavorable for water-wall interactions. However, a TBI with a small twist angle could potentially form on a templated rough surface, presenting an additional challenge for controlled growth of TBIs with small twist angles in nanocapillaries. In addition to surface roughness, we also examined the surface-water interactions by adjusting the carbon-oxygen interaction strength with ±20% variation. We found that TBIs can be stable for εC-O value being within the range from 0.07 kcal/mol to 0.11 kcal/mol (Supplementary Fig. 7). In addition, local defeats also show negligible effects on the overall twist angle with the hydrophobic confinement (Supplementary Fig. 8).

To further guide future experimental efforts, we constructed a semi-quantitative phase diagram (Fig. 4) based on the percentage of AA-stacking oxygen atoms (P), where AA-stacking was defined as oxygen atoms with a projected xy distance shorter than 0.5 Å (Supplementary Fig. 9). The phase diagram exhibits three distinct regions: 1) AA-stacking bilayer ice: Observed at confinement width h > 8.3 Å, with more than 40% AA-stacking oxygen atoms. 2) AB-stacking bilayer ice: Observed at narrower confinements (h < 8.0 Å at 26 nm⁻² and h < 7.8 Å at 30 nm⁻²) with less than 5% AA-stacking oxygen atoms. 3) TBI: Restricted to a narrow range of confinement widths (8.0 Å <h < 8.3 Å), highlighting the strong sensitivity of TBI formation to confinement conditions. These findings underscore the importance of precise control over nanocapillary width and lateral pressure for the experimental realization of TBIs. Additionally, to facilitate experimental validation through spectroscopy, we computed the vibrational density of states (VDOS) for TBI at twist angles of 21.8° and 27.8°, respectively, providing predicted spectral signatures for detection (Supplementary Fig. 10). These methods and data offer promising pathways toward synthesizing and probing TBI structures under controlled conditions.

Fig. 4: A semi-quantitative phase diagram of confined bilayer ice at room temperature.
figure 4

The phase diagram is plotted for confinement width versus area density. The AA phase (with probability P ≥ 0.4) is represented by the green region, the TBI phase by the blue region, and the AB phase (P ≤ 0.05) by the orange region.

Full size image

In conclusion, we conducted a series of classical and MLP-based MD simulations to investigate the formation of TBIs. We observed that TBI can form spontaneously under nanoconfinement at room temperature, exhibiting commensurate twist angles of 21.8° and 27.8°. This spontaneous formation of TBIs contrasts with the manual twisting required for the production of vdW materials, highlighting a unique property of hydrogen-bonded bilayer systems. The stability of TBI at the two twist angles is further supported by Gibbs free energy calculations, showing that 21.8° and 27.8° correspond to the lowest free-energy configurations. The thermodynamic stability of TBIs at both twist angles can also be attributed to their long-range order and the adaptability of hydrogen bonding. Notably, the interlayer hydrogen bonding in the link-like local structures is significantly strengthened, contributing to the overall stability. However, both TBI structures at θ = 21.8° and 27.8° exhibit an average of 3.5 hydrogen bonds per molecule (Supplementary Table 1), indicating that both TBIs do not fully satisfy the Bernal-Fowler ice rule. Overall, this work deepens our understanding of 2D ices and 2D non-vdW materials by demonstrating that TBIs can form spontaneously and exhibit marked stability at specific twist angles. These findings have important implications for the design and application of twisted 2D materials, particularly those involving hydrogen-bonding systems. We expect that the first simulation evidence of spontaneous formation of TBIs will open new research opportunities in the field of twisted non-vdW materials, with potential applications in nanofluidics and nanomaterial design.

Methods

Classical molecular dynamics

The simulation box dimensions were set to 83.64 Å ×85.2 Å ×18.0 Å. Two graphene sheets were positioned symmetrically at z = 0 and remained rigid throughout the simulations. The center-to-center distance between the nanocapillaries is varied from 7.5 Å to 8.5 Å, which is sufficiently wide to allow for the formation of bilayer ice. The number of water molecules N varies from 1600 to 2200, corresponding to an area density ρs from 22.45 to 30.87 nm-2. All simulations were performed using the LAMMPS package35 in the NVT ensemble at 300 K, with temperature control achieved via a Nosé-Hoover thermostat. The lateral pressure was maintained within the range of 2.0 to 4.5 GPa by varying the number of water molecules (Supplementary Fig. 2). A time step of 1.0 fs was used with the velocity-Verlet integrator and periodic boundary conditions were applied in the x and y directions. Fixed boundaries were maintained in the z-direction.

Water molecules were modeled using the TIP4P/2005 potential17,36, which accurately describes the phase behavior of bilayer ice. The interatomic interactions include a short-range Lennard-Jones (LJ) 12-6 potential and a long-range Coulomb potential. Long-range Coulomb interactions were computed using a particle-particle particle-mesh (PPPM) algorithm with an accuracy of 10−4. The cutoff distances were set at 12.0 Å for LJ interactions and 8.5 Å for Coulomb interactions. Interactions between water molecules and graphene were described using the LJ potential. The parameters for the interaction between the oxygen atoms of water and the carbon atoms of graphene were ϵco = 0.093 kcal/mol and σco = 3.19 Å37. The LJ parameters for carbon-carbon interactions were ϵcc = 0.047 kcal/mol and σcc = 3.22 Å. Simulation systems were initially relaxed for 5 ns to ensure system equilibration, followed by an additional 10 ns MD simulation for data collection.

MLP-based molecular dynamics

A machine learning potential (MLP) was developed using the n2p2 package with the Behler-Parrinello neural network framework. The Lennard-Jones (LJ) 9-3 potential was used to model the confinement interaction between graphene and water molecules, with parameters fitted against DFT calculations38,39. The training data, previously used to produce accurate phase diagrams of 2D ice, demonstrated good accuracy in energy prediction and force prediction, with root-mean-square errors (RMSE) of 1.2 meV and 72.2 meV/Å38, respectively. The Gibbs free energy was calculated as follows:

$$G=H-{TS}$$
(1)

where H is the enthalpy of the system, and the entropy S is estimated using the two-phase thermodynamics method40.