Competing hydrogen-bond orders drive water’s anomalous surface tension

Domain: physics
Published: 2026-02-20
Editorial outcome: major_revision

Abstract
Water’s surface tension shows a nonlinear temperature dependence, including a reentrant increase in the supercooled regime—a longstanding puzzle in physical chemistry. Using molecular dynamics simulations, we uncover a structural mechanism linking microscopic ordering to macroscopic interfacial behaviour. Surface tension arises from the interplay between ρ-states, characterised by O–H alignment under surface symmetry breaking, and tetrahedral S-states stabilised in the subsurface by negative pressure. Water’s surface tension γ is governed by the interplay of their anisotropies: at intermediate temperatures, ρ-state anisotropy saturates while S-states remain weakly anisotropic, slowing the growth of γ. Upon deeper supercooling, however, S-states acquire orientational order, amplifying anisotropy and producing the reentrant rise. This unified framework explains both inflection points of γ(T) and establishes a structural–mechanical link between local hydrogen-bond motifs and interfacial stress, with implications for nucleation, cryopreservation, and ferroelectric-like ordering, and extending beyond water to other network-forming liquids.

Introduction
Water’s surface tension underpins a wide range of physical, chemical, and biological phenomena1–5. It enables insects like water striders to walk on water6,7, shapes liquid droplets8,9, stabilises liquid bridges between solids10, and drives capillary flow in confined spaces such as plant vessels11–13. By tuning surface tension, it is possible to control bubble formation in gas-evolving reactions14,15 and manipulate fluid motion via the Bénard–Marangoni effect, where temperature gradients create surface tension gradients that drive flow16.

The physical origin of surface tension is the stress anisotropy created near the liquid surface2,17,18. At the interface, cohesive forces on molecules are unbalanced due to the reduced number of neighbouring molecules in the vapour phase compared to the liquid. This imbalance creates a tensile force that pulls molecules inward to minimise interfacial area. The magnitude of surface tension γ depends on the strength of intermolecular interactions. At ambient conditions, water has a relatively high surface tension of γ ≈ 72 mN/m19,20, significantly higher than most organic liquids, such as benzene (γ ≈ 28 mN/m)21. The high surface tension of water arises from the significant energy cost associated with disrupting the hydrogen-bonding network when a surface is formed. Because hydrogen bonds are substantially stronger than typical van der Waals or dispersion interactions, the loss of even a small number of hydrogen bonds at the air–water interface contributes disproportionately to the surface energy density. This explains why water exhibits a much higher surface tension compared with most organic liquids.

While most simple liquids exhibit a nearly linear increase in surface tension with decreasing temperature3, water deviates markedly from this trend. Experiments and simulations show that the rate of increase in γ gradually slows below ~275 K (see, e.g., ref. 22). More notably, in the deeply supercooled regime (T ≲ 250 K), simulations have revealed the emergence of a second inflection point, characterised by an accelerated rise in γ upon further cooling23–25. Experimental evidence for this anomaly remains scarce due to the difficulty of probing water at such low temperatures. However, recent high-precision measurements by Vinš et al.26, using a modified capillary rise technique down to 241.8 K, report subtle but systematic deviations from the smooth temperature dependence observed at higher temperatures, suggesting the onset of interfacial anomalies in this regime. This reentrant, nonlinear behaviour highlights that water exhibits fundamentally distinct interfacial properties from those of simple liquids at low temperatures. Previous studies have proposed diverse explanations, such as stronger and longer-lived hydrogen bonds25, preferential adsorption of high-density liquid (HDL)24, and a minimum in surface excess entropy linked to compact interfacial structuring23. While these interpretations offer important insights into interfacial thermodynamics and kinetics, they do not directly link molecular-scale structural transformations to the local stress anisotropy that defines surface tension, nor do they clarify how these features evolve spatially across the interfacial region. Moreover, previous interpretations have overlooked an intermediate temperature regime where the increase of γ with decreasing T is slower than in the high-T linear regime.

In contrast, our approach is grounded in the two-state model27–29, which describes water as a mixture of locally favoured tetrahedral (S-) and disordered (ρ-) state molecules. We show that the anomalous temperature dependence of surface tension arises from a structural competition between surface-induced OH/dipolar ordering—more easily adopted by the ρ-state—and the formation of tetrahedral S-state—stabilised in the subsurface region by negative pressure due to its larger specific volume. This competition governs interfacial stress anisotropy, the microscopic origin of surface tension.

In this work, using atomistic molecular dynamics simulations, we show that water’s surface tension is governed by the interplay of ρ- and S-state anisotropies. At intermediate temperatures, the ρ-state contribution saturates while interfacial S-states remain largely orientationally disordered and thus weakly anisotropic, moderating the growth of γ. Upon deeper supercooling, however, dipole–dipole interactions drive orientational ordering of S-states near the interface, amplifying anisotropy and producing the reentrant rise in γ. This surface-induced S-ordering also promotes the formation of Ice-0-like motifs with five- and seven-membered rings, thereby facilitating ice nucleation at air–water interfaces30. Our findings thus provide a unified microscopic framework for understanding the nonlinear temperature dependence of water’s surface tension, rooted in two-state structural competition. More broadly, this approach offers general insight into interfacial anomalies and phase transitions in complex liquids exhibiting local structural ordering31,32, with important implications for crystal nucleation and related interfacial phenomena.

Results
Nonlinear T-dependence of water’s surface tension
Our system consists of TIP4P/2005 water molecules33 forming a water film within a three-dimensional periodic box (Fig. 1a; see “Methods” section for additional details). We examine the T-dependence of γ. The TIP4P/2005 water model33 accurately produce γ across a wide temperature range, from T = 220 K to T = 550 K Fig. 1a, aligning well with prior simulation results17,25. Although minor deviations from experimental data remain, the TIP4P/2005 model reproduces the overall trend and magnitude of experimental γ(T)19,20,22, including the inflection behaviour, with satisfactory accuracy. Notably, a clear deviation from the linear trend in γ appears around 275 K, where its temperature dependence becomes nonlinear. Moreover, we reproduce the characteristic more rapid growth in γ for T ≲ 250 K (Fig. 1a), as reported in previous experimental26 and simulation studies23–25. This anomalous, nonlinear behaviour of γ contrasts with that of simple liquids such as Lennard-Jones model fluids34.

Distinct structural states of interfacial water
Hereafter, we analyse the surface tension anomaly based on the two-state model. As T decreases, many properties of water deviate from the normal behaviour of simple liquids39. These deviations, known as water’s anomalies, are widely attributed to the formation of locally favoured tetrahedral structures (LFTS) stabilised by hydrogen bonding27,37,40–47. The competition between two structural forms, LFTS (S-state; see the inset of Fig. 1b) and a disordered structure (ρ-state), has been incorporated into the two-state model27–29, providing physical explanations for the origin of water’s anomalies.

We classify water molecules into S-state and ρ-state using the translational order parameter ζ and orientational order parameter θavg. Here, ζ is defined as the difference between the minimum distance to a non-hydrogen-bonded neighbour and the maximum distance to a hydrogen-bonded neighbour37, while θavg represents the average angle formed between an oxygen atom and its neighbouring oxygens35,36. The parameter θavg captures orientational order within the first coordination shell35,36, while ζ reflects translational order in the second coordination shell37. Water is considered hydrogen-bonded if the donor(D)-acceptor(A) distance is less than 3.5 Å and the H-D-A angle is less than 30∘ (Luzar-Chandler parameters)48.

To demonstrate the presence of two distinct structural states in interfacial water, we analyse the bivariate distribution P(ζ, θavg)38 for water molecules located at z = 14–22 Å [Fig. 1c], where z denotes the distance from the centre of the water film. Consistent with findings in bulk TIP4P/2005 water38, we observe two distinct peaks separated by a pronounced valley [Fig. 1c for T = 300 K; also see Fig. S1 for T = 230–400 K], highlighting the two-state character of interfacial water. The peak near ζ ≈ 0.5 Å and θavg ≈ 109∘ corresponds to the S-state, while the other peak signifies the ρ-state. We remark that the probability P(ζ, θavg) is averaged over all water molecules and time frames. Our earlier work43,49 further analysed the rapid exchange dynamics between the ρ- and S-states, confirming that these two predominant motifs coexist in dynamic equilibrium. Considering the fundamental relation between the free-energy landscape and probability distributions,

F(ζ, θavg) = −kBT ln P(ζ, θavg), (1)

where kB is the Boltzmann constant, the clear bimodality of P(ζ, θavg) observed in Fig. 1 of ref. 38 and Fig. 1c supports the two-state description of both bulk and interfacial water as coarse-grained yet physically consistent representations. This indicates a strong correlation between translational and orientational ordering, supporting a two-order-parameter description31. The free-energy landscape F = −kBT ln P(ζ, θavg) (Fig. S1a) provides a thermodynamic basis for our two-state framework. Importantly, the bimodality of the distribution P(ζ, θavg) and the associated two-state features persist when an instantaneous liquid interface is employed (Willard–Chandler interface; Fig. S2)50. We further emphasise that these two-state features are not artefacts of sampling water molecules exclusively at the air–water interface. Even in the bulk liquid, where there is no net orientational ordering, we consistently observe the coexistence of two distinct structural motifs (Fig. S2), in agreement with our work38. Together, these results demonstrate that the S-state and ρ-state are intrinsic structural features of liquid water, present in both bulk and interfacial environments.

We note that our two-state model serves as a conceptual framework for partitioning fluctuating water structures into two distinct motifs, a framework long used to explain water’s thermodynamic anomalies27–29. Applied to surface tension, this approach allows us to disentangle the partial contributions of these motifs. Importantly, the partitioning is statistical: water molecules are not permanently assigned to a single state. Rather, the respective contributions are obtained by long-time averaging, ensuring thermodynamic consistency and robustness. Although recent simulation studies have questioned the two-state description at room temperature51–54, our work focuses on the supercooled regime, where the distinction between S- and ρ-states becomes clear and the dominant role of the S-state becomes critical.

Roles of ρ- and S-states in the surface tension anomaly
To probe the structural origin of the nonlinear γ(T), we decompose the surface tension into contributions from the S- and ρ-states (Fig. 1a). The S-state contributes substantially less than the ρ-state, even as its fraction increases at lower T, confirming its tension-reducing role. For example, at T = 250 K, the S-state accounts for only 13% of the total γ (Fig. 1a), despite representing approximately 37% of the bulk population (Fig. 1b). This disparity highlights that the symmetric hydrogen-bonded structure of the S-state, being mechanically rigid, produces lower interfacial stress as long as its dipoles remain randomly oriented.

The slowdown of γ between 275–250 K primarily reflects the saturation of the ρ-state contribution (Fig. 1a), as the surface layer becomes fully populated with aligned ρ-states whose anisotropy no longer increases upon cooling (see the inset of panel c in Figs. S3–S7). In this regime, the interfacial s varies only weakly (Fig. 1b), so the S-state contribution also remains nearly constant. Together, these effects lead to the plateau-like behaviour of γ(T).

Moreover, we reproduce the more rapid growth in γ for T ≲ 250 K (Fig. 1a), as reported in previous experimental26 and simulation studies23–25. This behaviour closely matches the inflection point at which the contribution from the S-state increases sharply (red line in Fig. 1a), suggesting that the stress anisotropy associated with the S-state is significantly enhanced in the deeply supercooled regime. We also note that the interfacial region consistently shows a lower S-state fraction (s) compared to the bulk (brown line in Fig. 1b), in agreement with experimental observations of reduced tetrahedral coordination at water interfaces55. Importantly, a sharp increase in interfacial s is observed around T ≈ 250 K (Fig. 1a), which coincides with the onset of the faster rise in γ(T). This reflects the fact that S-states, which generate little stress when randomly oriented, begin to develop orientational order in this regime (as discussed below), thereby enhancing anisotropy and contributing to the reentrant increase of γ.

Spatial distribution of tetrahedral water
To understand the spatial distribution of the S-state near the interface, we decompose the water density profile ρw(z) into contributions from the S- and ρ-states (Fig. 2a). This decomposition allows us to determine the fraction of the S-state, s(z), as a function of the distance z from the centre of the water film (Fig. 2b). Notably, we observe that s peaks at z ≈ 16 Å and decreases significantly at the interface. Furthermore, the variation of γ(z) with respect to z reveals that the stress anisotropy, which is the origin of surface tension2,3, extends down to 10 Å below the interface (Fig. 2c). In the bulk, stress anisotropy is absent due to the isotropic environment. The interfacial thickness estimated from the γ(z) profile (~10 Å) is slightly larger than the 6–8 Å reported by MB-pol simulations and second-order spectroscopic measurements56, and also larger than the 3.5 Å corresponding to two surface layers identified by DFT-based MD simulations57. This discrepancy may arise because our approach quantifies stress anisotropy, a mechanical property that reflects collective, non-local correlations, rather than structural anisotropy probed by optical methods. Stress-based measures may therefore exhibit a longer decay length compared to structurally defined interfacial thickness. Notably, we find that the S-state contributes significantly less to surface tension compared to the ρ-state at the interface. This is not solely because the S-state is less prevalent; we also confirm that the surface tension contribution per water molecule, γw, is consistently smaller for the S-state (Fig. 2d). These results confirm that the stable, symmetric tetrahedral structure of the S-state contributes to a balanced stress environment and thus reduces surface tension, provided their orientations remain random. We perform a similar analysis as in Fig. 2 (at T = 300 K) across a broader temperature range (T = 210–400 K), revealing qualitatively similar roles of the S-state (Figs. S3–S7). Water molecules in the S-state consistently exhibit weaker spatial stress anisotropy (smaller γw) and contribute less to the surface tension compared to those in the ρ-state.

However, we observe that below T ≲ 250 K, the S-state exhibits a sharp increase in γw(z) (Figs. S3–S7). This indicates that the S-state no longer mitigates the rise in surface tension. The trend is further supported by the evolution of the maximum peak γ* with respect to T for both states. Specifically, γ* grows markedly with decreasing T, particularly below T ≲ 250 K (Fig. S8f), consistent with the onset of accelerated growth in γ under deeply supercooled conditions. We attribute this enhancement of S-state stress anisotropy to the emergence of orientational ordering at low temperatures. Details will be discussed in the following sections.

Negative pressure enhancing tetrahedral ordering
We note that the local enhancement in tetrahedral order s in the interfacial region can be expressed as δs = speak − sbulk ≈ 1.3%, where speak and sbulk are the peak and bulk values of the profile Δs(z) (Fig. 3a). This enhancement can be partially attributed to the subsurface negative pressure P, located slightly beneath the surface [see the inset of Fig. 2d for P]. This mechanical negative pressure promotes formation of the S-state relative to the bulk, where P ≈ 0 and sbulk ≈ 24%. However, δs remains nearly constant across temperatures, even as the magnitude of the negative pressure increases at lower T (Fig. 3b). We attribute this constancy of δs to enhanced OH- and dipole-orientational orderings at lower T (see Fig. 3c; details of these orderings will be discussed later), which compete with the negative-pressure induced S-ordering. This scenario also explains the slight inward shift in the peak position of s (Fig. 2b) compared to the position of the negative pressure (the inset of Fig. 2d). Thus, the inward shift of the peak position—while its magnitude remains nearly constant—indicates that the balance between negative pressure and OH ordering is not determined solely by their relative amplitudes but also by their depth-dependent contributions across the interface.

Together, Figs. 2 and 3 reveal a spatial pattern of competing molecular orders at the air-water interface. Near the surface (z ≈ 20–22 Å), water exhibits strong OH-orientational and dipolar alignment, which is enhanced by broken symmetry and reduced local density. This ordering frustrates the formation of tetrahedral hydrogen-bond ordering, thereby suppressing the S-state in this region. Slightly deeper in the subsurface layer (z ≈ 16–20 Å), the local negative pressure stabilises S-state structures, although residual dipolar alignment still imposes orientational constraints. In the bulk-like region, both the OH-orientational ordering and stress anisotropy vanish, and S-states dominate with random orientations.

In the intermediate temperature range (275–250 K), the behaviour of the S-state reflects a balance between opposing effects. The extent of the S-enriched subsurface layer contracts (Fig. 3a), reducing the fraction of molecules that contribute to isotropy. At the same time, the fraction of S-states within the remaining layer rises moderately, mirroring the bulk T-dependence of s (Figs. 1b and S3b–S7b). This competition yields only a weak T-dependence of interfacial s. In contrast, once T falls below 250 K, the sharp increase of bulk s drives a corresponding rise in interfacial s (Fig. 1b).

Methods
Water modelling
We perform atomistic simulations of TIP4P/2005 water33, implemented in the Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS; version 29 Aug 2024)85. Our system contains Nw = 2112 water molecules arranged in a slab geometry: a 4 nm-thick water film confined in a three-dimensional periodic box of dimensions Lx = Ly = 4 nm and Lz = 12 nm (Fig. 1a). This setup creates two water-air interfaces along the z-axis, separated by an 8 nm vacuum region. The vacuum spacing is sufficiently large to eliminate spurious interactions between the two surfaces, allowing each interface to be treated as effectively independent. Surface symmetry breaking arises naturally from the abrupt density gradient at the interface and the anisotropic boundary conditions imposed along the z-direction. We use the SHAKE algorithm86 to constrain bond lengths and angles, ensuring the rigidity of water. The system is evolved using the velocity-Verlet algorithm with a time step 1 fs. Electrostatic interactions and dispersion interactions are calculated via the particle–particle particle–mesh algorithm87 with a relative accuracy of 10−5. The temperature T is maintained using the Nosé–Hoover thermostat with damping parameter 102 fs88 in the canonical (NVT) ensemble. We conduct simulations for 1 ns to equilibrate the system, followed by the production run of 200 ns. For temperatures of T ≤ 250 K, we reduce Lx = Ly by a factor of 2 and decrease the number of water molecules to Nw = 1056, while maintaining the film thickness and box size Lz. The system is equilibrated for 20 ns, followed by an extended production run of 200 ns (220 K ≤ T ≤ 250 K) and 400 ns (T = 210 K). Configurations are recorded every 103 fs for sampling the equilibrium properties. The surface tension is calculated using the mechanical definition17,89,

γ = (1/2)Lz [Pzz − (1/2)(Pxx + Pyy)]. (2)

Here, σαα = (1/V) Σi σiαα, where σiαα are the diagonal components of the per-atom stress tensor of atom i and V = Lx Ly Lz. The microscopic stress tensor was calculated for each atom using the compute stress/atom command in LAMMPS (version 29 Aug 2024)85 (https://docs.lammps.org/compute_stress_atom.html), which implements the Irving–Kirkwood (IK) formulation. The details of how LAMMPS computes stress under periodic boundary conditions are described in refs. 90,91. We acknowledge that the local stress profile is not uniquely defined and depends on the specific method used to localise intermolecular forces92,93. However, our primary conclusions are based on the qualitative features and temperature dependence of the stress profile rather than on the precise quantitative value at any single spatial position. As a result, modest quantitative variations in γ(z) arising from different definitions of the local stress do not affect the overall physical interpretation or the conclusions of this work.

Ring analysis
Ring analysis94 was performed by first constructing a hydrogen-bond network connecting each oxygen atom to its nearest neighbours. For each water molecule i, six triplets are defined, each consisting of i and a pair of its hydrogen-bonded neighbours. A ring is then identified as the shortest closed path that includes all molecules within a triplet. The local topology around a molecule is characterised by the size distribution of these rings94. Ice-0 is distinguished by two characteristic local environments around the central water: (5,5,5,5,6,6) and (5,5,5,7,6,6)30,94. A molecule is classified as “Ice-0–like” if at least five of the six surrounding rings match these signatures.

Data availability
All data that support the findings of this study are included in this paper and its Supplementary Information. Source data are provided with this paper.

Code availability
The LAMMPS input files used for the simulations are provided as Supplementary Software. Additional analysis scripts for calculating the physical quantities defined in the figures are available from the corresponding author upon request.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_3": "In this manuscript Tanaka and co-workers have examined the behavior of the surface tension of the air-water interface in an effort to rationalize/explain its apparent anomalous behavior under supercooling. While the experimental literature observation on this remains inconclusive, there have been several suggestions/indications on the existence of an inflection point in the surface tension and its relationship to thermodynamic anomalies.\nWhile the paper tackles an interesting problem in the physical-chemistry/chemical-physics of liquid water and at aqueous interfaces, I am not entirely convinced about its urgency or broad impact. I believe it sits much better in a more specialized/technical journal like JCP or JPC. In fact, the authors site a previous paper on the topic published in JCP - while the approaches are different, the two works mirror each other quite closely (https://pubs.aip.org/aip/jcp/article/158/5/054503/2871767/High-density-liquid-HDL-adsorption-at-the). Specifically, the paper uses standard models (TIP4P) and standard analysis tools - rings, parameters that have been pioneered by Tanaka. All very interesting, but still standard and would find a better place in specialized journal.\nI now raise my specific issues with the paper which I hope will be useful for the authors for a future publication elsewhere:\n1) The authors use TIP4P/2005 water model. This is fine and perfectly okay. How much do their results translate to another more accurate model such as MB-pol. I am not suggesting they repeat all their analysis with MB-pol but some discussion on the limitations are needed. Furthermore, some benchmarks with a polarizable model would also make a much more convincing case.\n2) More conceptually, I am really not convinced by the s-vs-p state water partitioning that is done here. It seems like a real desire to fit into previous two-state models of water invoked to try to rationalize the putative second critical point of water/and anomalies. In particular, have the authors looked at the dynamics (correlation function) associated with the theta and eta parameters that they claim yield a two-state underlying thermodynamics. Related to this point, the authors should prove using a committor analysis that they have the correct free energy landscape - Figure 1b. Until these dynamical properties are examined which can unequivocally show the two-state behavior the arguments cannot be used to rationalize the fluctuations at the air-water interface.\n3) Related to point 2), there have been several studies challenging the notion of 2-state water at least at room temperature and in particular the sensitivity of the results to the choice of collective variable/order parameter - the authors should look and comment on the following papers - https://pubs.acs.org/doi/10.1021/acs.jctc.1c01292, https://pubs.acs.org/doi/10.1021/acs.jpcb.2c07169, https://pubs.acs.org/doi/10.1021/acs.jctc.2c01205, https://pubs.acs.org/doi/10.1021/acs.jpclett.4c00383).\n4) Related to my previous point - what is known from SFG spectroscopy on the water orientations as a function of supercooling? I know the authors look at the dipole orientations and make some connections with SFG data - however, wouldn't a simpler argument for all this analysis simply be rooted in fluctuations in one free energy minimum that shifts as a function of temperature?\n5) I have not done a thorough literature review but I think there are papers that should be examined more carefully on the topic of surface tension and how it changes under supercooling - on such paper that is missing from the current paper are the following - https://pubs.acs.org/doi/full/10.1021/acs.jpclett.0c01163, https://pubs.aip.org/aip/jcp/article-abstract/150/23/234507/197649/Surface-tension-of-supercooled-water-nanodroplets?redirectedFrom=fulltext, https://pubs.acs.org/doi/abs/10.1021/acs.jced.7b00519)\n6) One other concern/issue that I have related to the first point is the role of NQEs (nuclear quantum effects) on the apparent anomalies which may become more enhance under supercooling.", "reviewer_4": "The authors study the temperature-dependence of the water-vapor interfacial tension using molecular dynamics simulations deep into the supercooled regime. The observed non-linear behavior including two inflection points is rationalized in terms of the competition of two distinct water populations (tetrahedral vs. random), of which the composition at the interface (and in bulk) changes with temperature. The surface tension profile agrees well with the differences in water ordering, thus the authors conclude that the competition between dipolar ordering and tetrahedral structure formation governs the resulting surface tension.\nThis conclusion is appealing since it rationalizes important aspects of water’s anomalies merely in terms of structural properties, while also not contradicting thermodynamic arguments involving a liquid-liquid critical point. The topic is of high interest and well written, however, especially the second half of the discussion (after the discussion of Figure 3) is very technical. While personally I find all points generally interesting and relevant, I am not sure if this is too detailed for the general audience of Nature Communications and thus the manuscript in its present form might be better suited for a more specialized journal.\nDetailed comments that should be addressed before publication:\n- page 3 „such as plant vessels such as plant vessels“: repetition\n- the inflection points in Fig1a are barely visible, maybe they could be indicated with arrows?\n- the abbreviation LFTS is never introduced and cannot be anticipated with regards to a broad readership\n- The definition of a hydrogen bond in the present context is missing (I assume the Luzar-Chandler parameters are used?)\n- The definition of the S-state in terms of theta, xi is not given\n- Details on the calculation of the microscopic stress tensor are missing, this is highly relevant since different routes (eg. IK vs. Harasima) might lead to different microscopic details (see 10.1021/acs.jctc.6b00576, 10.1063/5.0132487)\n- Fig.3 d,e are not clear, what message should be taken here?\n- detection of ice0-like structures: the manuscript should stand on its own and therefore at least briefly outline how these structures are detected. The same applies to the m-membered rings.\n- A slab system is simulated - has any correction for the 3d-periodic PPPM like the one proposed by Yeh and Berkowitz been employed?\n- The surface tension contribution per molecule in Fig. 2d should have units", "reviewer_1": "In this paper, the authors use molecular dynamics simulations of the air-water interface to investigate the nonlinear behavior of its surface tension as a function of temperature. Using the two-state model of water, they carefully analyze the depth-dependent microscopic structural motifs, aiming to link them to macroscopic thermodynamic properties of the interface. From this analysis, they reveal a competition between surface-induced OH ordering and the preferred tetrahedral motifs that minimize surface tension, with the latter becoming increasingly populated at lower temperatures. However, at lower temperatures, they also identify a mechanism by which these tetrahedral structures no longer mitigate the surface tension through orientational ordering. With this insight, the authors clearly identify a correlation between the change in anisotropic structure with temperature and the corresponding inflection points in the surface tension behavior, thereby suggesting a mechanism connecting the microscopic structure to the macroscopic properties. Overall, I think the manuscript offers clear insight into this important topic with interesting implications for the structural precursors for controlling ferroelectric ordering in the super-cooled regime. It is thus a welcome contribution to the field. Nevertheless, I have several points that I think require attention and revisions to the manuscript before it is publishable in Nature Communications.\n\n• In their paper the authors use a 2-state model and extract partial contributions to the surface tension from 2 different water structures. To me it is not entirely clear how this decomposition works for water as it assumes the presence of 2 distinct and static structural motifs. In liquid water, however, these structures are at best limiting cases of a broad distribution which is additionally highly dynamic. Furthermore, the description of water with a two-state model has been discussed in the literature for a while and appears very controversial. It would be nice if the authors could include a paragraph that addresses this controversy and justify the applicability of the two-state model for their study.\n\n• Figure 1 – Regarding part (a), the authors show the temperature dependence of the surface tension, comparing their simulated results to those from previous simulations as well as experimental measurements. This plot is an essential part of the paper as it should demonstrate the deviations of the T-dependent surface tension from a straight line. Unfortunately, the way the curves are presented this deviation is very difficult to make out for the reader. Furthermore, with the presented scale for the surface tension, it is even more difficult to identify the described inflection points. While I acknowledge that the authors want to show the breakdown into the S- and ρ- contributions, it would be very helpful for the reader if they make use of the full y-axis scale to show the temperature dependence. One suggestion would be to have a break in the axis from 20-45 mNm-1 so that the s-state contribution can remain at the bottom, but then use the remaining space to show the 45-90 mNm-1 range.\n\n• Figure1 – For the comparison between the experimental and simulated data the authors make two claims: (1) the error bars are smaller than the symbol size, and (2) the calculated values align well with the experimental values. Without further clarification of this point, it might seem that these two statements are contradictory given that there are clearly observable deviations between the two datasets which must therefore fall significantly outside several standard deviations defined through the error margins. The authors should clarify this point more carefully to clearly describe what they define as ‘aligning well’.\n\n• Figure1 – It is difficult to correctly interpret the evolution of S- and p- contributions with temperature as they include both, the temperature dependence of the molecular contributions to the surface tension (given in the SI) and the temperature dependent population ratios. The latter are given in figure 1b, however, only for bulk water. There is as far as I could see no graph presented that shows the temperature dependent population ratio between S- and p-state water for the relevant interfacial region. It would be very helpful to the reader if the authors could include such a figure.\n\n• Page 3 – in the first paragraph, the authors refer to the “thermal Marangoni effect”. It would be more appropriate to refer to this as the “Bénard–Marangoni effect”.\n\n• Page 3 – The authors discuss the origin of the relatively high surface tension of air-water, attributing it to the “strong electrostatic interactions and the directional nature of hydrogen bonds”. I think it would be worth clarifying that the bond energy of the hydrogen bonds, in relation to other intermolecular interactions, is also a substantial factor here given that the surface tension represents a surface energy density that is contributed to by the loss in H-bonding.\n\n• Figure 2 – In this figure, as well as in the text, the density is given the symbol ρ. As this is also used for one of the two states of water (ρ-state), I would recommend the authors change one or the other for better clarity.\n\n• Figure 2 – I think the comparison between the plots in Figure 2 with those in Figures S2-S6 could perhaps benefit from consistent axis scaling.\n\n• Figure 2 – For clarity, I think it would perhaps be clearer to not use black and blue for the pressure contributions in the inset to part (d) given that these are used to represent the total and s-state contributions elsewhere in the figure. The same obviously applies to figures S2-S6.\n\n• Page 8 – the authors claim that the anisotropy extends “10 Å below the interface”. I think two points of clarification are warranted here. Firstly, describing it as extending “below” the interface is perhaps somewhat misleading given that this 3D region somehow defines the interface. This length-scale is often described relative to the Gibbs-dividing-surface (GDS), but this exists part-way through this anisotropic region. Secondly, this value of 10 Å is higher than values quoted elsewhere, particularly in relation to the structural anisotropy probed by second-order optical methods. For example, the DFT-based MD simulations in ref. 57 show only ~two contributing ‘layers’ or around 3.5 Å. Alternatively, other recent MD simulations using MB-pol along with experimental measurements using second-order optics concluded around 6-8 Å [Langmuir 2024, 40, 35, 18760–18772]. It would be helpful to the reader if the authors put their obtained value into the context of these other measurements, particularly highlighting why their value might differ with it being not based on the structural anisotropy probed by second-order methods.\n\n• Page 8 – I find the last paragraph generally lacking in clarity and justification. For example, it is not clear to me why the stress anisotropy should be quantified by the maximum peak in the γw profile – is this not essentially only considering the molecules with the largest contribution to the surface tension at the expense of the (like more numerous) other contributors? It is also not clear to me what the justification is for attributing the sharp rise in stress anisotropy for the S-state to the orientational ordering at low T. I think general readers would benefit from additional clarification in this section.\n\n• Page 9 – I think the explanation for why the competition between enhancement of Δs due to the negative pressure and the enhancement of OH ordering leads to the essentially constant, but inward shifting peak position needs more clarity. Unless I am mistaken, the increase in OH ordering essentially contributes negatively to Δs, while the negative pressure contributes positively to it. Therefore, if both increase equally, no change in the profile should occur. The fact that the peak position is shifted to lower z-positions whilst maintaining the overall peak height then suggests not only a change in amplitude of these two factors, but also in their depth-dependence.\n\n• Page 9 and Figure 3 – The authors should clarify the definition of Δs. On page 9 it is defined as the difference between the peak and the bulk value, therefore being independent on z, while in figure 3 it is defined differently as s(z) – s bulk. Based on the discussion on page 9 about having a constant Δs for different temperatures, the former definition seems to be more appropriate. Therefore, the authors should change figure 3(a) and its caption accordingly.\n\n• Figure 3 – I find parts (d) and (e) here not at all clear in terms of what they are trying to display / show. I suggest the authors make changes to this figure to make it clearer the purpose of these sub-figures.\n\n• Figure 4 – In addition to figures 5, 6, S8, and S9, I think the way the structural motifs of different molecules are presented by showing the 2D plots of one O-H angle vs. the other are not at all intuitive. This is largely due to the non-independent relationship between the two angles. A good example of this is the isotropic distribution shown in Figure 4(f). When one thinks of an isotropic distribution in terms of molecular orientations, one naturally imagines a perfectly smooth plane where all orientations are equal. With this cos(θ 1) vs cos(θ 2) way of showing the orientations, this isotropic picture is far from what is obtained. It seems to me that the desired structural information could equally be obtained by describing the structures in terms of their two independent Euler angles: tilt of the molecular dipole, and twist around the dipolar axis. This description can be found in recent simulation work on the air-water interface by the Netz group [arXiv:2505.12962 and arXiv:2505.19856]. With such a description, the isotropic distribution yields a perfectly smooth image that seems much more intuitive to the reader. Therefore, while I realize that a reasonably large effort is required in changing the figures from the current description to the tilt-twist description, I think it would be highly valuable for the understanding and general readability of the paper. It also has the additional benefit that the tilt angle then describes perfectly the distribution of dipole orientations.\n\n• Page 11 – I think the authors should refine their statement on the water orientation at the interface. Currently, they state that “the OH bonds in the top layer orient upward towards air”. While I agree that similar statements are often used to describe the air-water interface, this is not an accurate representation. Instead, the top layer has an average dipole orientation pointing towards air, but nevertheless has a broad distribution. Therefore, simply stating that the OH bonds point towards air is overstating the population of these moieties. Furthermore, as shown in the orientational profiles, the most populated structural motif in the ‘top layer’ is something akin to #2, which only has a single of its OH groups directed towards air. The statement as it stands, however, makes it sound like both OH groups are pointing towards air, more akin to structure #1. While these may seem like subtle points, it completely changes the ‘picture’ of the interfacial structure’\n\n• Page 12 – The authors discuss the non-negligible proportion of in-plane water molecules. Firstly, I find it very difficult with the current way this is presented (see above comment) as well as the chosen color-scale to identify this “non-negligible” population. As above, I recommend that the authors present this differently. Secondly, the authors discuss the energetic favorability of the in-plane structural motif, but make no contrast to other thermodynamic factors. After all, the in-plane structures seem not to be the dominant orientation, so if they are energetically favorable compared to pointing up/down dipoles, the authors need to clearly qualify why they are not the predominant structure.\n\n• Page 13 – The authors state that they “detect a significant number of water molecules displaying structures indicative of dipole alignment”. I would again suggest that this can be more clearly shown with the tilt-twist description given that the distribution of tilt vs z perfectly describes the extend of dipolar alignment.\n\n• Page 13 – The authors state that the “S-state consistently exhibits less pronounced OH-orientational order, characterised by a more dispersed P(cos(θ 1), cos(θ 2))”. Firstly, I don’t think this is immediately obvious from the figures in the way they are currently presented. Take for example Figure 5(a) and (f). It is hard to qualitatively compare the orientational order simply by looking at these two plots. There should be a way to better quantify this and provide a more robust statement on this point. Secondly, the disperse nature of P(cos(θ 1), cos(θ 2)) is not automatically connected to orientational order. An example of this is a distribution with two preferred, but different, orientations. This can make it seem more disperse while actually the extent of the orientational order could be very high.\n\n• Figure 7 – In parts (e) and (f), don’t use pink as the color for 16-17Å since this color represents 21-22Å in (c) and (d), with all the other colors being consistent.\n\n• Figure 9 – The authors state that “it is evident that the S-state is particularly favoured in forming 5-, 6-, and 7-membered rings,…”. Given that f S < f ρ I would suggest rewording this as the rings are actually strictly favored by the ρ-state. Unless I am mistaken, what the authors mean here is that the 4-membered ring state is the only state showing a lower population than the average, indicating it is less energetically favorable than the average s-environment.\n\n• Page 16 – when discussing the orientational distributions of the H-bonded O-Hs, the authors show that it is a broad, symmetric distribution around 0. Firstly, I am surprised that the distribution is symmetric. Given the asymmetry of the z-dependent anisotropy at the interface, I would naturally expect the distribution to reflect this. Perhaps the authors could add a comment on this. Secondly, the authors state that these distributions indicate a “low probability of being dangling”. I find this somewhat of an irrelevant statement since they are discussing the orientational distribution of H-bonded O-H groups, i.e. one should expect that they are not ‘dangling’ simply by definition.\n\n• Page 17 – The authors claim that they “find that the S-state is more favourable for the formation of these nucleation precursors”. I would argue on the contrary as the ρ-state has far higher f-fractions for the ring structures. I therefore think that the authors should expand on this point.\n\n• Page 20 – In the conclusions, the authors propose the possibility of validating these simulation results using second-order spectroscopy methods. I completely agree that this would be a valuable approach to test the proposed structural mechanism. As such, it would also be interesting to calculate the depth-dependent χ(2) in such low temperature conditions, particularly splitting the contributions for the s-state and ρ-state water. This would give clear indicators for the potential of future experimental measurements. While I guess this is beyond the scope of this work, it might be useful to comment on this as future potential.", "reviewer_2": "I co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is part of the Nature\nCommunications initiative to facilitate training in peer review and to provide appropriate recognition for Early Career\nResearchers who co-review manuscripts."}, "author_response": ""}, {"round": 2, "reviewer_comments": {"reviewer_3": "I have re-examined the paper by Tanaka and co-workers. It has certainly improved and the authors have made significant efforts to improve the clarity of the discussions and make a stronger connection and distinction with previous works in the field. Overall I am very happy with the changes made. Having looked at the paper again now with the changes I can definitely see that it could land a place in Nature Communications. There are still some issues that I have that the authors must address before it is considered for publication:\n1) While I understand the need to use larger boxes and longer simulation times, there are now off-the-bench models trained with neural networks that allow for doing quantum-quality with polarizable models, simulations like the ones done here. The authors should take a look for example at the following work that studies the surface tension with an NN potential to see if they can at least make a comment to connect with that result - https://pubs.aip.org/aip/jcp/article-abstract/162/24/244705/3350869/Charge-transfer-at-air-water-interfaces-A-machine?redirectedFrom=fulltext\n2) Regarding the ferroelectricity issue of supercooled water, this is another interesting point that deserves much more background discussion. In particular, there are suggestions by previous Russian groups suggesting that the putative LLCP in supercooled water involves a ferroelectric transition (https://link.springer.com/article/10.1134/S0036024411050220). Then more recently, there have also been a series of papers investigating the possible Ferroelectric origins in bulk supercooled water - https://www.pnas.org/doi/10.1073/pnas.2407295121, https://www.pnas.org/doi/10.1073/pnas.2412456121).\n3) Regarding Figure R1 in the responses, how was it constructed exactly? Was it made using \"ALL\" the water molecules in the box or only those at the interface? How does that same plot look like for bulk water (assuming one constructs the same thing looking at the orientational polarization along a specific direction even in the bulk?). As already pointed out in previous papers and in my critique in the first round, the free energy distributions along specific order parameters/variables in reduced dimension can lead to all sorts of artefacts (due to dimensionality reduction). I am basically not convinced that the 2 minima in the FES shown are really proof of \"two-state\" or just because when you get to the interface their is net orientation/polarization. If you had any surface/interface like a hard wall, you would also see the same thing. There have been many studies examining the water topology at the air-water interface - like waters lying parallel or perpendicular. Are the features the authors present here some how related to those? See work by Adam Willard on 3-body defects at the air-water interface.\n4) Related to the point 3, did the authors look at the FES as a function of interfacial distance for example from the Willard-Chandler interface? How does the \"bi-modality\" in the FES look like as a function of proximity to the surface?\nOnce the authors clarify/address these points, I would be happy to consider for publication.", "reviewer_4": "The authors have well addressed my previous comments and together with the incorporation of comments provided by the other reviewers the quality and readability have considerably improved.\nOnly a small detail already pointed out in my previous comments and which should be stressed for clarity remains regarding the computation of a local stress profile - this is by no means unique, see e.g. https://doi.org/10.1063/5.0132487 and https://doi.org/10.1021/acs.jctc.6b00576: I must admit, however, that the slight quantitative changes one would expect in the profile would not change the qualitative conclusions made in the manuscript.", "reviewer_5": "The authors have carefully considered the comments of all reviewers, including Reviewers #1 and #2. In this way, they have improved and streamlined their manuscript, with a clear focus on providing a microscopic explanation for the nonlinear temperature dependence of water’s surface tension. With their intention to stress the novelty and relevance of their microscopic framework – statements with which I agree – some text on these aspects may even have become too redundant.\nA suggestion by Reviewers #1 and #2 that the authors did not take up concerns using Euler’s tilt and twist angles to characterize molecular orientations rather than employing the theta1 and theta2 angles to describe the orientations of the two O-H bonds of a water molecule (see Fig. 4). Both representations are valid alternatives. I can understand the authors’ arguments for sticking with the theta1-theta2 representation because it has some advantages in the present context.\nAnother suggestion that the authors did not implement relates to using identical axis scaling in Figs. S2-S6. I agree with the authors that such uniform scaling would infer with clear visibility of the data. However, contrary to what was promised in their reply and what is useful, the authors did not explicitly specify the respective axis scaling in the captions (“To aid cross-comparison, however, we have explicitly noted in the revised captions the numerical ranges of each axis.”)\nAdditional comments:\n- Fig. S4 displays an incorrect graphic.\n- Fig. S12: (i) Unlike stated, there is no inset in panel (e). (ii) What atoms define the vector connecting two hydrogen-bonded water molecules? (iii) “… the peak in P(cos theta_H) occurs at a nonzero value”: there are actually two peaks."}, "author_response": "We thank all reviewers for their careful reading of our manuscript and for their constructive\nand insightful comments. We have carefully considered all suggestions and revised both the\nmanuscript and Supplementary Information accordingly.\nFor the reviewers’ convenience, all modifications in the revised manuscript and\nSupplementary Information are highlighted in blue. Detailed point-by-point responses are\nprovided below.\nI have re-examined the paper by Tanaka and co-workers. It has certainly improved and the\nauthors have made significant efforts to improve the clarity of the discussions and make a\nstronger connection and distinction with previous works in the field. Overall I am very happy\nwith the changes made. Having looked at the paper again now with the changes I can\ndefinitely see that it could land a place in Nature Communications. There are still some issues\nthat I have that the authors must address before it is considered for publication:\nWe sincerely thank Reviewer 3 for his/her positive reassessment of our revised manuscript.\nWe are pleased that Reviewer 3 recognises our efforts to improve the clarity of the\ndiscussion and agrees that the manuscript is potentially suitable for publication in Nature\nCommunications. We appreciate Reviewer 3 taking the time to provide the additional\ninsightful comments.\nIn the following, we address Reviewer 3’s comments one by one and highlight the\ncorresponding revisions implemented in the manuscript.\n1) While I understand the need to use larger boxes and longer simulation times, there are\nnow off-the-bench models trained with neural networks that allow for doing quantum-quality\nwith polarizable models, simulations like the ones done here. The authors should take a look\nfor example at the following work that studies the surface tension with an NN potential to\nsee if they can at least make a comment to connect with that result -\nhttps://pubs.aip.org/aip/jcp/article-abstract/162/24/244705/3350869/Charge-transfer-at-air-\nwater-interfaces-A-machine?redirectedFrom=fulltext.\nWe thank the reviewer for drawing our attention to recent advances in machine-learning\npotentials (MLPs) for water, which indeed represent an important step toward quantum-\nmechanical accuracy at reduced computational cost. We have now cited the suggested work\nand added a discussion in the revised manuscript to clarify how our results relate to these\nsimulations.\nIn particular, the MLP study focuses on the high-temperature regime (300 K and above) and\ndoes not report surface-tension calculations. In contrast, our simulations span a much wider\ntemperature range, extending into the deeply supercooled regime down to 210 K, where the\nanomalous and non-linear behaviour of surface tension is most pronounced. Simulations\nwith quantum-quality MLPs at such low temperatures, system sizes, and time scales remain\ncomputationally prohibitive at present.\nBelow we reproduce the relevant revision made in the manuscript.\nFor example, recent simulations using machine-learning potentials (MLPs) for water with\nquantum-mechanical accuracy have reported that water molecules in the outermost layer in\ncontact with air tend to be positively charged, while those in a sub-interfacial layer are\n2\nnegatively charged, reflecting interfacial charge transfer effects [81]. These simulations are\nparticularly valuable in that they provide microscopic structural information at the\nquantum-mechanical level while being significantly faster than conventional density\nfunctional theory calculations.\nAlthough our results are obtained using a classical water model, we observe qualitatively\nconsistent interfacial behaviour: the net dipole moment near the interface points toward the\nvapour phase (inset of Fig. 2c), and the outermost layer is enriched in dangling OH bonds\n(Fig. S9, Supplemental Information, at 300 K). The presence of these dangling OH bonds is\nstructurally consistent with the positive charging tendency reported in Ref. [81], as the\nexposed hydrogen atoms carry partial positive charge. Together with the close quantitative\nagreement between our calculated surface tension γ and experimental values (Fig. 1a), these\nresults support the suitability of the TIP4P/2005 model [33] for investigating the structural\nand thermodynamic properties of interfacial water.\nWe note that the MLP study [81] focused on interfacial structure at temperatures of 300 K\nand above and did not report surface-tension calculations. In contrast, the present work\naddresses the temperature dependence of surface tension over a broader range, including\nthe deeply supercooled regime down to 210 K.\n2) Regarding the ferroelectricity issue of supercooled water, this is another interesting point\nthat deserves much more background discussion. In particular, there are suggestions by\nprevious Russian groups suggesting that the putative LLCP in supercooled water involves a\nferroelectric transition (https://link.springer.com/article/10.1134/S0036024411050220).\nThen more recently, there have also been a series of papers investigating the possible\nFerroelectric origins in bulk supercooled water -\nhttps://www.pnas.org/doi/10.1073/pnas.2407295121,\nhttps://www.pnas.org/doi/10.1073/pnas.2412456121).\nWe agree with the reviewer that the possible ferroelectric ordering of supercooled water is\nan important and timely topic that merits broader discussion, particularly in the context of\nthe deeply supercooled regime addressed in this work. We have therefore expanded the\nbackground discussion to place our results in the context of earlier theoretical proposals and\nrecent simulation studies.\nSpecifically, we now explicitly discuss early suggestions by Russian groups that the putative\nliquid–liquid critical point (LLCP) in water may be accompanied by, or closely related to, a\nferroelectric transition, as well as more recent molecular dynamics studies that have\ninvestigated ferroelectric correlations and polarisation instabilities in bulk supercooled water\nusing both classical and machine-learning-based potentials. In our revised manuscript, we\nclarify how these ideas connect to our two-state framework and emphasise that our results\nidentify the tetrahedral S-state as a structurally well-defined motif that naturally supports\nshort-range dipolar correlations.\ninterface. Rather, we propose that the S-state provides a microscopic structural basis for local\ndipolar correlations, which may act as precursors to ferroelectric ordering under conditions\nof deep supercooling, confinement, or reduced orientational frustration. The revised text\nsituates our findings as complementary to LLCP- and ferroelectric-transition scenarios,\noffering a real-space structural perspective that links tetrahedral ordering, dipolar\ncorrelations, and interfacial stress anisotropy.\nThere is a growing body of evidence suggesting the possible existence of a ferroelectric\nphase transition in supercooled water [71-74]. Early theoretical work proposed that\nsupercooled water may support a ferroelectric state, driven by collective dipolar ordering\nwithin the hydrogen-bond network [71]. More recently, molecular dynamics studies based\non the classical TIP4P/Ice model have investigated the ferroelectric origins of bulk\nsupercooled water, proposing that the ferroelectric transition and the liquid–liquid phase\ntransition (LLPT) can be viewed as two facets of the same underlying phenomenon governed\nby coupling between density and polarisation degrees of freedom [74]. Consistent with this\npicture, simulations employing deep neural-network potentials have shown that the low-\ndensity liquid exhibits a strong propensity toward spontaneous polarisation, providing\nfurther support for the connection between tetrahedral ordering and dipolar correlations\n[73].\nOur results offer a complementary, structurally resolved perspective on this emerging\npicture. We find that a robust short-range dipolar correlation consistently exists for the\ntetrahedral S-state (the low-density component) in both liquid water and ice (Fig. S18). The\nmagnitude of this correlation is comparable to that observed in crystalline ice, indicating\nthat S-state water molecules form locally polar structural motifs even in the absence of long-\nrange ferroelectric order. Importantly, these dipolar correlations between S-structures can\npropagate in space at low temperatures, where the S-state becomes the dominant species,\nthereby facilitating the emergence of longer-ranged ferroelectric ordering.\nThus, our results provide a plausible structure-based explanation for the reported\nferroelectric ordering in supercooled water [72–74], without interfering with the existence\nof a liquid–liquid critical point (LLCP). Within this framework, ferroelectric correlations\nand liquid–liquid criticality can naturally coexist, with the tetrahedral S-state acting as the\nmicroscopic structural motif that links tetrahedral ordering, dipolar correlations, and\nmacroscopic polarisation tendencies. Moreover, ferroelectric ordering can be viewed not\nmerely as a byproduct of the liquid–liquid transition, but as an active driving mechanism\nthat promotes it through the cooperative growth of dipolar correlations, pointing to a unified\nphysical picture in which polarisation and density ordering are intrinsically coupled and\nmutually reinforcing.\n3) Regarding Figure R1 in the responses, how was it constructed exactly? Was it made using\n4\n\"ALL\" the water molecules in the box or only those at the interface? How does that same\nplot look like for bulk water (assuming one constructs the same thing looking at the\norientational polarization along a specific direction even in the bulk?). As already pointed\nout in previous papers and in my critique in the first round, the free energy distributions\nalong specific order parameters/variables in reduced dimension can lead to all sorts of\nartefacts (due to dimensionality reduction). I am basically not convinced that the 2 minima\nin the FES shown are really proof of \"two-state\" or just because when you get to the interface\ntheir is net orientation/polarization. If you had any surface/interface like a hard wall, you\nwould also see the same thing. There have been many studies examining the water topology\nat the air-water interface - like waters lying parallel or perpendicular. Are the features the\nauthors present here some how related to those? See work by Adam Willard on 3-body\ndefects at the air-water interface.\n4) Related to the point 3, did the authors look at the FES as a function of interfacial distance\nfor example from the Willard-Chandler interface? How does the \"bi-modality\" in the FES\nlook like as a function of proximity to the surface?\nWe thank the reviewer for these careful and technically important questions regarding the\nconstruction and interpretation of the free-energy surface (FES) analysis. We address points\n(3) and (4) together below.\nConstruction of Fig. R1 and choice of molecules.\nFigure R1 (now included as Fig. S1 in the Supplemental Information) was constructed by\nsampling all water molecules located within the interfacial region, defined as the slab\nbetween Aand A(as specified in the caption of Fig. 1). The free-energy\nsurface is defined through the standard statistical-mechanical relation,\n𝑧𝑧 = 14 ˚ 𝑧𝑧F=(ζ2 ,2 θ ˚ ) = − k T ln P (ζ , θ ),\navg B avg\nwhere and quantify local translational (density-related) and angular (tetrahedral) order,\nrespectively, as explicitly defined in the main text. The bimodality observed in both the\n𝜁𝜁 𝜃𝜃avg\nprobability distribution (Fig. 1c) and the corresponding FES (Fig. S1a) reflects the\npresence of two well-defined basins separated by a free-energy barrier under the studied\n𝑃𝑃(𝜁𝜁,𝜃𝜃avg)\nconditions, indicating the coexistence of two distinct structural motifs within the interfacial\nregion. Furthermore, we confirm the two-state features persist even when employing an\ninstantaneous liquid interface (Willard-Chandler interface; incorporated into the Fig. S2 in\nthe supplemental information).\nRole of dimensionality reduction and interface-induced artifacts.\nWe fully agree that free-energy projections onto reduced sets of order parameters must be\ninterpreted with care, as dimensionality reduction can in principle generate apparent features.\nFor this reason, we have performed several consistency checks. Crucially, the bimodal\nfeatures do not originate from orientational polarisation induced solely by the presence of\nthe interface. Even in the bulk liquid—where there is no net orientational ordering—we\nobserve the same two coexisting structural basins in (Fig. S2d), consistent with\n5\n𝑃𝑃(𝜁𝜁,𝜃𝜃avg)\nour previous work [Ref. 38, J. Phys. Chem. B 127, 3452–3462 (2023)]. This demonstrates\nthat the S-state and -state are intrinsic structural features of liquid water rather than artifacts\nof interfacial polarisation.\n𝜌𝜌\nDependence on distance from the interface.\nTo further address the reviewer’s question regarding spatial localisation, we examined the\nFES as a function of proximity to the interface using both fixed- slabs and an instantaneous\nWillard–Chandler interface. The two-state character persists when the latter definition is\nemployed (Fig. S2), confirming that the observed bimodality is n𝑧𝑧ot an artifact of a particular\ninterfacial definition. As one approaches the interface more closely (e.g., – A; Fig.\nS1b), the two-dimensional distribution becomes more smeared compared to the bulk case,\nreflecting a gradual reduction in the population of the S-state due to the r𝑧𝑧ed=uc1e8d 2hy2dr ˚ ogen-\nbond coordination near the interface. Nevertheless, the two distinct basins remain clearly\nidentifiable, supporting the robustness of the two-state framework. This continuous\nevolution of populations is also consistent with the depth-dependent behaviour shown in Fig.\n2b.\nRelation to orientational ordering and prior interfacial studies.\nWe thank the reviewer for pointing out prior work on water topology at the air–water\ninterface, particularly the study by Shin and Willard (J. Phys. Chem. Lett. 9, 3181–3186\n(2018)), which demonstrated that three-body hydrogen-bond defects play a key role in\ndetermining interfacial polarisation and dielectric response. We view our results as\ncomplementary to this body of work. While previous studies have primarily characterised\norientational features (e.g., parallel versus perpendicular alignment relative to the surface\nnormal), our two-state framework instead classifies local environments based on hydrogen-\nbond connectivity and tetrahedral order.\nIn this context, the orientational anisotropy reported in earlier studies can be understood as\na geometric consequence of the interface selectively stabilising different intrinsic structural\nspecies. The highly ordered S-state tends to preserve bulk-like hydrogen-bond connectivity,\nwhereas the more disordered -state naturally accommodates broken bonds and dangling\nOH groups, which are known to contribute strongly to interfacial polarisation. Thus, our\nframework does not replace ori𝜌𝜌entational descriptions but provides a structural underpinning\nthat links orientational polarisation to local tetrahedral order.\nOverall interpretation.\nTaken together, these analyses demonstrate that the two minima observed in and\nreflect intrinsic structural coexistence in liquid water rather than artifacts arising\n𝑃𝑃(𝜁𝜁,𝜃𝜃avg)\nfrom dimensionality reduction or interfacial sampling. By systematically tracking how the\n𝐹𝐹(𝜁𝜁,𝜃𝜃avg)\npopulations of these states evolve with temperature and proximity to the interface, our\napproach provides a thermodynamic and structural basis for understanding how hydrogen-\nbond connectivity and dipolar ordering jointly influence the anomalous temperature\ndependence of surface tension in deeply supercooled water.\nbond defects. The following text has been added to the main text to clarify our viewpoint:\nThe free-energy landscape (Fig. S1a) provides a thermodynamic\nbasis for our two-state framework. Importantly, the bimodality of the distribution\n𝐹𝐹 = −𝑘𝑘𝐵𝐵𝑇𝑇𝑇𝑇𝑇𝑇 𝑃𝑃(𝜁𝜁,𝜃𝜃𝑎𝑎𝑎𝑎𝑎𝑎)\nand the associated two-state features persist when an instantaneous liquid\ninterface is employed (Willard–Chandler interface; Fig. S2) [50]. We further emphasise that\n𝑃𝑃(𝜁𝜁,𝜃𝜃𝑎𝑎𝑎𝑎𝑎𝑎)\nthese two-state features are not artefacts of sampling water molecules exclusively at the air–\nwater interface. Even in the bulk liquid, where there is no net orientational ordering, we\nconsistently observe the coexistence of two distinct structural motifs (Fig. S2), in agreement\nwith our previous work [38]. Together, these results demonstrate that the S-state and ρ-state\nare intrinsic structural features of liquid water, present in both bulk and interfacial\nenvironments.\nPrevious work [66] has shown that non-ideal hydrogen-bond geometries are required to\naccount for the depth-dependent variation of orientational polarisation at the air–water\ninterface at room temperature and identified triangular three-body hydrogen-bond defects\nas key structural elements. We propose that the two states identified in our work provide a\nnatural structural basis for these observations. The S-state, representing the highly ordered,\nlow-density component, acts as a structural precursor that largely preserves ideal\ntetrahedral hydrogen bonding. In contrast, the ρ-state (the high-density component)\ncorresponds to disordered local environments characterised by non-ideal hydrogen-bond\ngeometries and associated defects, which contribute significantly to the net orientational\npolarisation at the interface. In this sense, our two-state framework extends previous\nstructural interpretations by systematically exploring their temperature dependence down to\nthe deeply supercooled regime and by providing a structural foundation for understanding\nthe anomalous behaviour of water’s surface tension.\nPrevious studies have also characterised the air–water interface through analyses of\nmolecular orientation [60, 61, 66]. Our results are fully consistent with these findings while\noffering a complementary, topological perspective. Whereas orientational metrics quantify\nhow individual molecules align relative to the interface normal, the S-state and ρ-state\nframework captures the connectivity and local hydrogen-bonding environment. For example,\nthe broken hydrogen bonds associated with perpendicular or “dangling” water molecules\nare characteristic of the disordered ρ-state. Thus, our two-state description deepens earlier\norientational analyses [60, 61, 66] by providing a thermodynamic and structural link\nbetween orientational ordering and local tetrahedral ordering at the air–water interface.\nOnce the authors clarify/address these points, I would be happy to consider for publication.\nWe hope that the reviewer finds the revised manuscript suitable for publication in Nature\nCommunications.\nThe authors have well addressed my previous comments and together with the incorporation\nof comments provided by the other reviewers the quality and readability have considerably\nimproved.\nOnly a small detail already pointed out in my previous comments and which should be\nstressed for clarity remains regarding the computation of a local stress profile - this is by no\nmeans unique, see e.g. https://doi.org/10.1063/5.0132487 and\nhttps://doi.org/10.1021/acs.jctc.6b00576: I must admit, however, that the slight quantitative\nchanges one would expect in the profile would not change the qualitative conclusions made\nin the manuscript.\nWe sincerely thank Reviewer 4 for their careful reading of our manuscript and for their\npositive assessment. We are pleased to hear that, following revision, the quality and\nreadability of the manuscript have been considerably improved.\nWe also thank the reviewer for the final technical comment regarding the computation\nof the local stress profile. We fully agree that the local stress tensor is not uniquely\ndefined and that its form depends on the specific procedure used to localise\nintermolecular forces. We have now explicitly acknowledged this non-uniqueness in the\nmanuscript and cited the references suggested by the reviewer [91, 92], which discuss\nthis issue in detail.\nWe acknowledge that the local stress profile is not uniquely defined and depends on the\nspecific method used to localise intermolecular forces, as discussed in detail in previous\nstudies [92, 93]. However, our primary conclusions are based on the qualitative features\nand temperature dependence of the stress profile rather than on the precise quantitative\nvalue at any single spatial position. As a result, modest quantitative variations in ( ) arising\nfrom different definitions of the local stress do not affect the overall physical interpretation\nor the conclusions of this work. 𝛾𝛾𝑧𝑧\nWe hope that the reviewer finds the revised manuscript suitable for publication in Nature\nCommunications.\nThe authors have carefully considered the comments of all reviewers, including Reviewers\n#1 and #2. In this way, they have improved and streamlined their manuscript, with a clear\nfocus on providing a microscopic explanation for the nonlinear temperature dependence of\nwater’s surface tension. With their intention to stress the novelty and relevance of their\nmicroscopic framework – statements with which I agree – some text on these aspects may\neven have become too redundant.\nWe thank Reviewer 5 for their careful evaluation of our work and for recognising the\nimprovements made in response to the comments of all reviewers. We are pleased that the\nmanuscript is now clearer and more streamlined, with a focused presentation of our\nmicroscopic framework for explaining the non-linear temperature dependence of water’s\nsurface tension. We also appreciate the reviewer’s acknowledgement of the novelty and\nrelevance of this framework.\nA suggestion by Reviewers #1 and #2 that the authors did not take up concerns using Euler’s\ntilt and twist angles to characterize molecular orientations rather than employing the theta1\nand theta2 angles to describe the orientations of the two O-H bonds of a water molecule\n(see Fig. 4). Both representations are valid alternatives. I can understand the authors’\narguments for sticking with the theta1-theta2 representation because it has some advantages in the present context.\nWe thank Reviewer 5 for their thoughtful comment regarding the representation of molecular\norientations. We appreciate the recognition that both Euler tilt–twist angles and the –\nangles describing the orientations of the two O–H bonds are valid alternatives. As noted by\nthe reviewer, we chose to retain the – representation because it provides a more 𝜃𝜃d 1 ire𝜃𝜃c 2 t\nand physically transparent connection to hydrogen-bond geometry and orientational ordering,\nwhich is central to our discussion of 𝜃𝜃th 1 e 𝜃𝜃d 2 istinction between the two structural motifs.\nAnother suggestion that the authors did not implement relates to using identical axis scaling\nin Figs. S2-S6. I agree with the authors that such uniform scaling would infer with clear\nvisibility of the data. However, contrary to what was promised in their reply and what is\nuseful, the authors did not explicitly specify the respective axis scaling in the captions (“To\naid cross-comparison, however, we have explicitly noted in the revised captions the\nnumerical ranges of each axis.”)\nThe reviewer is entirely correct in noting that, although we stated in our response that the\nnumerical ranges of the axes would be specified in the figure captions, this detail was\ninadvertently omitted in the revised manuscript. We apologise for this oversight. We have\nnow revised the captions of Figs. S3–S7 in the Supplementary Information to explicitly state\nthe numerical ranges of each axis.\n- Fig. S4 displays an incorrect graphic.\nWe sincerely apologise for this oversight. The correct version of Fig. S4 (now Fig. S5 since\nwe add a new Fig. S2 in response to the comment from Reviewer #3) has now been included\nin the revised Supplementary Information.\n- Fig. S12: (i) Unlike stated, there is no inset in panel (e). (ii) What atoms define the vector\nconnecting two hydrogen-bonded water molecules? (iii) “… the peak in P(cos theta_H)\noccurs at a nonzero value”: there are actually two peaks.\nWe thank the reviewer for identifying these issues in Fig. S12 (now Fig. S13 since we add a\nnew Fig. S2 in response to the comment from Reviewer #3). (i) The reviewer is correct that\nno inset is present in panel (e). We have therefore removed the reference to the inset from\nthe caption. (ii) In Figs. S13e–f, the vector connecting two hydrogen-bonded water\nmolecules is defined using the positions of the oxygen atoms. (iii) We also acknowledge that\nthe distribution exhibits two non-zero peaks rather than a single one. We have\ncorrected the text accordingly.\n𝑃𝑃(cos𝜃𝜃H)\nThe presence of two non-zero peaks indicates that water molecules in the topmost interfacial\nlayer form hydrogen bonds both upward (towards the vapour phase) and downward (towards\nthe sub-surface liquid). This bidirectional bonding geometry is consistent with the distorted\nhydrogen-bond network characteristic of the air–water interface. The caption of Fig. S13 has\nbeen revised to clarify these points.\n…pointing from the oxygen atom of the central water molecule to the oxygen atom of its\nhydrogen-bonded water neighbour…\n…Notably, the two non-zero peaks in…\nWe hope that the reviewer finds the revised manuscript suitable for publication in Nature\nCommunications.\nThe authors have addressed all my comments and concerns. I believe the manuscript is now\nready to be accepted in Nature Communications and I am happy to recommend for\npublication.\nWe are pleased that Reviewer 3 finds the revised manuscript suitable for publication in\nNature Communications, and we thank the reviewer for their positive assessment.\nResponses to the comments of Reviewer 4\nThe authors have well addressed all my points raised in the previous report. In line with the\nincorporation of all the comments provided by the other reviewers the quality and readability\nhave considerably improved and I believe the manuscript is now suitable for publication in\nNature Communications.\nWe are pleased that Reviewer 4 finds the revised manuscript suitable for publication in\nNature Communications, and we thank the reviewer for their constructive feedback.\nResponses to the comments of Reviewer 5\nThe authors have carefully addressed my (and other) concerns and, in this way, further\nimproved their manuscript, rendering it suitable for Nature Communications.\nThere is only a minor point relating to the added Fig. S2, which should be clarified before\npossible publication: The probability distributions shown in Fig. S2 (b)-(d) imply that the S-\nstate (at theta=109°) has a higher probability at the interface (L1) than in the interior (L3),\nwhile the results of the main text, in particular, Fig. 1(b) show the opposite.\nWe are pleased that Reviewer 5 finds our manuscript suitable for publication in Nature\nCommunications, and we thank the reviewer for this careful and helpful comment.\nTo ensure a direct and consistent comparison with Fig. 1b of the main text, we explicitly\nidentify the outermost sublayer of the instantaneous interface, denoted L , based on\n1A\ncomparison with the ideal ice–air interfacial structure (Ref. 63). When this definition is used,\nwe find that the fraction of the S-state in the outermost interfacial layer is lower than in the\nbulk, in full agreement with Fig. 1b of the main text.\nFigure S2 in the Supplementary Information has been updated accordingly to clarify this\npoint and to avoid any potential ambiguity in layer identification.\n2"}, {"round": 3, "reviewer_comments": {"reviewer_3": "The authors have addressed all my comments and concerns. I believe the manuscript is now ready to be accepted in Nature Communications and I am happy to recommend for publication.", "reviewer_4": "The authors have well adressed all my points raised in the previous report.\nIn line with the incorporation of all the comments provided by the other reviewers the quality and readability have considerably improved and I believe the manuscript is now suitable for publication in Nature Communications.", "reviewer_5": "The authors have carefully addressed my (and other) concerns and, in this way, further improved their manuscript, rendering it suitable for Nature Communications.\nThere is only a minor point relating to the added Fig. S2, which should be clarified before possible publication:\nThe probability distributions shown in Fig. S2 (b)-(d) imply that the S-state (at theta=109°) has a higher probability at the interface (L1) than in the interior (L3), while the results of the main text, in particular, Fig. 1(b) show the opposite."}, "author_response": ""}]

Decision letter