Liquid structure adjacent to solid surfaces follows the superposition principle

1 Liquid structure adjacent to solid surfa ces follows the superposition principle Qian Ai 1,2 , Haiyi Wu 3 , Lalith Krishna Samanth Bonagiri 2,4 , Kaustubh S. Panse 1,2 , Shan Zhou 1,2 , Fujia Zhao 1,2 , Yitong Li 1,2 , Kenneth S. Schweizer 1,2,5,6 , Narayana R. Aluru 3 & Yingjie Zhang 1,2,7 * 1. Department of Materials Scienc e and Engineering, University of Illinois, Urbana, Illinois 61801, United States 2. Materials Research Laboratory, University of Illinois, Urba na, Illinois 61801, United States 3. Walker Department of Mechanical Engineer ing and Oden Institute for C omputational Engineering & Sciences, The University of Te xas at Austin, Austin, Texas 78712, United States 4. Department of Mechanical Science and Engin eering, University of Il linois, Urbana, Illinois 61801, United States 5. Department of Chemistry, University of Illinois, Urbana, Illinois 61801, United States 6. Department of Chemical and Biomolecular Engin eering, University of Illinois, Urbana, Illinois 61801, United States 7. Beckman Institute for Advan ced Science and Technology, Univ ersity of Illinois, Urbana, Illinois 61801, United States *Correspondence to: yjz@illinois.edu Summary paragraph: Liquid structure at solid–liquid interf aces is critical for many natural and engineered processes ranging from biological signal transductio n to electrochemical energy conversion 1,2 . Advanced experimental and computational methods have provided insights into the structure of liquids adja cent to planar substr ates at the nanoscale 3,4 . However, realistic solid–liquid interfaces are inevitably inhomogene ous across multiple length sc ales, presenting a complexity that surpasses the capa bilities of existing approaches 5 . Here we bridge the complexity gap by discovering and utilizing a hitherto hidden principle of interfacial liquid—superposition. Experimentally, we use 3D atomic force microsco py (3D-AFM) to image the interfacial structure of a wide range of organic an d aqueous solvents and electrolytes, uncovering universal liquid density oscillations and emergent liquid layer reconfigurations at heterogeneous substrate sites. We further develop an analytical mod el, coined so lid–liquid superposition (S LS), which solves the interfacial liquid density distribution based on a key descript or: the effective total corre lation function (ETCF) between a liquid molecule and nearby solid atoms. SLS not only explains all the experimentally observed interfacial liquid distribution profiles fro m the angstrom to near- micron scale, but also predicts more precise atomic -scale interference patterns which are fu rther corroborated by molecula r dynamics (MD) simulations. This study unveils a key structural descriptor of interfacial liquids, and establishes a theoretical framework fo r rapidly and accurately predicting liquid structures adjacent to solid surfaces with arbitrary morpho logy and size scale. A major challenge in our understa nding of solid–liquid interfaces is the complex ity gap. Real- world interfaces in fields such as biology, renewable energy, and water purification are 2 unavoidably heterogeneous across ato mic, nano, micron, a nd larger scales 6,7 . Existing computational and experimental methods have enabled the accurate determination of th e multiscale structure of s olid surfaces 6,8 . However, our knowledge of interfacial liquids has been limited to those at ato mically flat or nanoscale substrates 9,10 . Computational approaches face an accuracy–scalability trade-off: the most accurate simulations (e .g., ab initio molecular dynamics) are restricted to nanoscale systems 11 , while scalable models (e.g., m ean-field approaches such as the Gouy-Chapman-Stern theory 12 ) are highly inaccurate within the crucial first 1–2 nm from the solid surface 11,13,14 . Experimentally, existing too ls are constrained by the resolu tion–complexity trade-off: high-resolution imaging methods (e .g., 3D-AFM) face ch allenges when measuring complex interfaces w ith large roughness 2,9,10 , while interfac e-sensitive spectroscopie s (Raman, infrared, X-ray, etc.) can probe la rge, complex solid–liquid in terfa ces but la ck spatial resolution 15– 17 . To bridge the complexity gap, here we take a descriptor-based appr oach, with the goal of understanding and predicting the interfacial liquid structure th rough analytical formulas and parameters that are applicable regardless of th e morphology and size scale of th e adjacent solid. We first experimentally image the liquid struct ure at mo del solid surfaces with contro lled heterogeneity. Based on the observe d key features, we propose struct ural descriptor s and integrate them into an analytic al model, SLS, to explain the in terfacial liquid struct ure. We then verify the predictive capabilities of SLS for arb itrary, complex solid–liquid interfaces. Experimental interfacial liquid structure at crystalline and atomic step sites We start by conducting 3D-AFM imaging (AC m ode, amplitude modulat ion) of a model system, pure diethyl carbonate (DEC ) solvent at highly oriented pyrol ytic graphite (HO PG) surface. This system both mimics the nega tive electrode–electrolyte interfa ces in lithium-ion batteries and hosts a series of well-defin ed substrate heterogeneities 18 . During 3D imaging, the AFM cantilever’s amplitude and phase ar e recorded as a function of x, y, and z, from which we construct the background-removed conservative force ( Δ F) 2,19–21 (see Methods). At flat, crystalline basal plane areas, we observe damped oscillations in the Δ F vs z (vertica l distance from the HOPG surface) profile (Fig. 1a,b). The overall damped os cillation patte rn is consistent with previously reported interfacial liquid structur e of pure water , aqueous solutions , or ganic electrol ytes, and ionic liquids on flat substrates, as va lidated by 3D-AFM, X-ray scatteri ng, and/or atomistic simulations 19 . Liquid density oscillation near solid surfaces is a uni versal phenomenon, aris ing from the finite size of the liquid molecules together with th e entropic constraints imposed by the substrate 2,19,22,23 . The oscillation period of DEC/HOPG is found to be nearly constant at ~4.0 Å (Fig. 1b), which corresponds to the center-to-center distance of adjacent DEC molecules 24 . To investigate the impact of surface heterog ene ity on the nearby liquid structure, we scan the HOPG surface and identify a series of clean atomic step sites, including mono-step and multi-step edges (Extended Data Fig. 1). 3D imaging is then p erformed at th ese sites. In the 3D-AFM community, both the local phase minima and for ce maxima have been conventionally used as indicators of liquid layer positions, with phase ha v ing higher signal-to-noise ratio and f orce being 3 more quantitative 2,9,19–21 . From the x-z force/phase maps (Fig . 1c–e, Extended Data Fig. 2, and Supplementary Figs. 1–3), we observe remarkable laye r tran sitions across the s tep edge areas, with two key features: 1) The liquid layers above the u pper HOPG surface (left of the step edge in Fig. 1c–e) “slip” at the step edge, connecting to a nearby liquid laye r above the lower HOPG surface; we name this effect as “layer crossover”, as the liquid layer number chan ges across the step edge. 2) Above the lower surface, the liqu id layers at z positions lower than that of the upper HOPG surface terminate near the step edge, although the exact termination spot cannot be determined in the bi-step and tri-step sy stems due to the finite size of the AFM tip (Extended Data Fig. 3 and Supplementary Note 1). Overall, the liquid layers exhibit non-conformal coverage near the step edge site—they do not follow the topogr aphic contour of the substrate. Fig. 1 | 3D-AFM force maps at different sites of the DEC/HOPG interface. a , Flat area x-z force ( Δ F) map. HOPG surface is marked by the gray line at the bottom, while DEC liquid layers are labeled by the corr esponding layer numbers. b , Force vs z profile by averaging along the x direction in a . c – e , 3D-AFM x-z force maps of DEC near HOPG atomic step sites. Each gray line represents a graphene layer as part of the HOPG substrate. Step heights are observed as 0.34 nm, 0.68 nm, and 0.94 nm for mono-step ( c ), bi-step ( d ), and tri-step ( e ), res pectively. DEC layers on each side of the step edge are labeled by the corresponding layer numbers, highlighting the layer crossover effects. In c , the regions exhibiting layer crossove r and terminati on effects are also marked by the dark red and pi nk dashed rectangles, respectiv ely, to guide the eye. In d and e , the blank regions (white, triangular-shaped) adjacent to th e step edges correspond to areas the AFM tip cannot access due to its finite radius (Extended Data Fig. 3 and Supplemen tary Note 1). In e , the top layer is bent downward near the step edge , which is common for gra phitic structures (see more discussions in Supplementary Note 2). Scale bars: 3 nm for x and 0.5 nm for z ( a , c – e ). In prior reports fro m our lab, both co nformal and non-conformal inte rfacial liquid layers have been observed for aqueous solutions and lithium sa lt/ionic liquid electrolyt es on substrates with 4 heterogeneous adsorbates or clus ters 2,25 . Other labs have also reported 3D-AFM imaging of aqueous solution at non-flat substrates 26–28 . However, due to the ra ndomness of the surface corrugation and/or the limited re solution/sensitivity, the und erlying mechanism has b een elusive. Here the DEC layers are resolved at the we ll- defined atomic step sites of HOPG. The layer crossover effect is iden tified unambiguously. For example, laye r 1 on the left crosses into layer 2 on the right for the mono- and bi-ste p systems; for the tri-step syst em, layer 1 crosse s into layer 3 from left to right (Fig. 1c–e). Within the connected liquid lay ers, we also observe a shift of force values across the step edge, furt her confirming the laye r crossover effect (Supplementary Fig. 4). The finite AFM tip radius has impacts on the e xperimentally observed tr ansition region width of both the substrate and the liquid layers, but does not affect the qualitative la yer crossover features (Extended Data Fig. 3 and Supplementary Note 1). We have measured DEC/HOPG (flat area and step edge) systems through five separate sets of experiments using two different types of cant ilevers. Each experimen t was conducted with a new AFM tip having a specific (likel y different) tip radius. In total, we have obtained more than 16,000 x-z maps and >1 million force curves. The force curves are hi ghly reproducible at all the flat area sites (Extended Data Fig. 4), while the layer crossover and termination effects are consistently observed at atomic step sites in all the measurements (Fig. 1c–e, Extended Data Fig. 2, and Supplementary Figs. 1–3). The observed non-conformal effect indicates that th e interfacial liquid structure cannot be predicted by simply following the surface cor rugati ons of the substrate. To unravel the universal relationship between substrate morphology and interfacial liquid structure, a new analytical framework is needed—one that extends beyond mode ls developed for flat solid–liquid interfaces. Solid–liquid superposition model To capture the key mechanism of the interfacial liquid s tructure, we take a descrip tor-based analytical approach. We consider the exact ato mic stru cture of the solid, and para meterize the solid-modulation of adjacent liquid density via an ETCF, ℎ  󰇛 𝒓  󰇜 . Here, 𝒓  𝒓  𝒓  , with 𝒓  representing the geometric center of a liqu id molecule and 𝒓  the location of a solid site (an ato m, lattice, or other repeating unit). At th is stage, we focus on the liquid density profile, as represented by the spatial distribution of the molecules’ ge ometric centers. This simplification is universa lly applicable to liquid molecules rega rdless of thei r shape, with th e limitation that the molecular orientational informat ion is not captured. The real total correlation function (TCF) be tween the liquid molecules and a solid site represents the local real-space liquid density prof ile. At heterogeneous sites where translational symmetry is broken, TCF is subject to change at different sites. Therefore, TCF cannot serve as a universal descriptor of the inter facial liquid structure. To overcome th is limit, here we define ETCF to capture the effective contribu tion of individual solid sites to the variation in nearby liquid density. For the same solid site–liquid molecule pair, ETCF is identica l, regardless of the local solid morphology. Therefore, ETCF can serve as a powerful descriptor for deriving the interfacial liquid structure. 5 For example, in the DEC/HOPG sy stem, if we take a carbon at om (of HOPG) as a site, the ETCF between the carbon site and nearby DEC molecu les is identical everywhe re at the interface. This includes the flat basal plane areas, step edges, and other heterogeneous sites where only carbon atoms are present. If impurity atoms (e.g., oxyge n and nitrogen) exist, different ETCFs may be needed. In this work, we mainly consider syst ems where only one ETCF is requ ired to describe the key interfacial liquid profiles. As shown in the following results, such systems can be highly heterogeneous and even include mu ltiple electrolyte components. Fo r systems requir ing more than one ETCF, our developed analytical model can also be readily ge neralized to obtain their liquid density distributions. A bulk liquid in equilibrium has a uniform time-averaged density, 𝜌  , due to the mobile nature of the liquid species. For a given ETCF, we propose that the time -averaged, near-surface liquid density variations can be obt ained through a linear superposition: Δ𝜌 󰇛 𝒓  󰇜 𝜌   ℎ  󰇛 𝒓  󰇜 𝒓  , 󰇛 1 󰇜 where Δ𝜌 󰇛 𝒓  󰇜 𝜌 󰇛 𝒓  󰇜 𝜌  , and 𝜌 󰇛 𝒓  󰇜 is the liquid density at position 𝒓  . Eq. 1 is the core of our proposed SLS model. The profound impact is, as long as the solid atom/site distribution and ETCF are known, in eith er analytical or nu merical form, we can instantaneously calculate the liquid density distribution adja cent to a solid with arbitrary surface morphology. ℎ  󰇛 𝒓  󰇜 can be obtained via multiple alternative approaches, such as analytical models, MD simulations, and/or experimental inputs (e.g., 3D -AFM and X-ray scattering). In the following, we show how ETCF can be qua ntified by integrating an analytical formula with either 3D-AFM data or MD simulation results. Based on classical theory of liquids, the TCF of pure liquids and liquid mixtures are usually in the form of damp ed oscillations, due to the dense packing of the disc rete molecule s 29–31 (Supplementary Note 3). Considering the generic oscillations and the impact of s ubstrate–liquid interactions, we propos e an analytical form of ETCF: ℎ  󰇛 𝒓  󰇜 𝐴  exp 󰇛 𝛼  𝑟  󰇜 cos 󰇛 𝛼  𝑟  𝜃  󰇜 / 𝑟  𝐴  exp 󰇩  󰇛 𝑟  𝑟  󰇜  2 𝜎  󰇪 , 󰇛 2 󰇜 where the first term is hypothesized based on pred ictions of the Ornstein-Zernike (OZ) equation 32 , and the second term is an empirical corre ction of the first liquid density peak. 𝐴  , 1/ 𝛼  , 2 𝜋 / 𝛼  , and 𝜃  correspond to the amplitude, decay length, pe riodicity, and phase f actor for the damped oscillation function, respectively. 𝐴  , 𝑟  , and 𝜎 represent the amplitude, center and width of the Gaussian function for first peak correction, respectivel y (see Supplementary Note 3 for more details). Combining experimental 3D-AFM force maps of the DEC/flat HOPG interface with Eqs. 1 and 2, we extract parameters fo r the ETCF between the DEC mol ecule and carbon atom of HOPG (Supplementary Notes 4–6 and Supplementary Tabl e 1). The obtained ETCF, shown in Fig. 2a, exhibits strong damped oscillations with a periodic ity of 4.0 Å, consistent with the 3D-AFM results. We then construct the atomic c onfigurations of HOPG substrates , including both the basal plane 6 and step edge areas (Fig. 2b, bottom and Supplem entary Figs. 5–7). Substituting the carbon atomic positions ( 𝒓  ) and the DEC–carbon ETCF ( ℎ  󰇛 𝒓  󰇜 ) into Eq. 1, we calc ulate the 3D DEC density maps at the interface (Fig. 2b–i). We refer to this an alytical calculation method as “AFM-SLS”, since it uses 3D-AFM data to extract the ETC F descriptor and then employs the obtained ETCF for SLS calculations (see summary in Supplementary Fig. 8). Fig. 2 | SLS model of interfacial liquid struct ure: key features and benchmarking. System: DEC/HOPG interface. a , The ETCF for DEC/HOPG, extracted from experimental 3D-AFM maps of DEC/flat HOPG interfaces. Inset is a schematic illustrating that r sl corresponds to the distance between a carbon atom (in HOPG) and th e geometric center of a DEC molecule. b , SLS-calculated flat area density map, obtained by plugging in the at omic positions o f four graphene layers (shown in the bottom, as model for HOPG basal plane) and the ET CF in a to Eq. 1. c , Liquid density variation ( Δ𝜌 / 𝜌  ) vs z curves of DEC/flat HOPG obtai ned from different methods. AFM-SLS curve (black) is obta ined by averaging along the x direction in b . 3D-AFM experimental density curve (blue) is converted from the measured force curves (Suppl ementary Note 5) . MD curve (red) is extracted by aver aging along x in the MD-sim ulated x-z density map (Extended Data Fig. 5a). 7 d – f , j – l , Liquid density maps from AFM-SLS model ( d – f ) and MD simula tion ( j – l ) near HOPG mono-step ( d , j ), bi-step ( e , k ), and tri-step ( f , l ). Gray dots in d – f and pink lines in j – l correspond to the HOPG substrate. Red dashed lin es in d – f , j – l mark the boundaries of the tran sition region where the liquid layers are non-flat (processing method shown in Supplementary Fig. 10). Full 3D atomic structures of th e constructed solids in d – f are shown in Supplementary Figs. 5 and 7. g – i , m – o , Horizontal density profiles along the liqui d regions marked by the corresponding white dashed rectangles in d – f , j – l , respectively. The maximal Δ𝜌 / 𝜌  value within the local z range is plotted at each x. Scale bars ( b , d – f , j – l ): 0.5 nm for both x and z. At the planar basal plane area, AFM-SLS predicts flat, layered distribution of the interfacial liquid with a vertical periodicity of 4.0 Å (Fig. 2b), c onsistent with 3D-AFM imaging results (Fig. 1a). This spatial pattern directly reflects the superposition of ETCF over the substrate carbon atoms: 1) superposition from the dense, periodic packin g of carbon atoms result s in uniform in-plane density distribution of the liquid, since the interm olecular distance of DEC (4.0 Å) is much larger than the atomic bond length of HOPG (1.4 Å); and 2) the oscillating profile of ℎ  󰇛 𝒓  󰇜 translates into the oscillation of liquid dens ity along the z direction (with the same periodi city) after the linear superposition (see derivations in Supplementary Note 4). To benchmark the AFM-SLS results, we furthe r perform MD simulati on of DEC/flat HOPG (system configuration shown in Supplementary Fig. 9a). We av erage the AFM-SLS x-z density map (Fig. 2b) over x to obtain the density vs z curve, and compare with results from 3D-AFM experiments and MD simulation (Fig. 2c). The c onversion from 3D-AFM fo rce to liquid density is discussed in Supplementary Note 5. We observe remarkable consistency in both the qualitative damped oscillation pattern and the quantitative de nsity values among these three methods. These results verify the validity of the SLS method. MD simulation further re veals the anisotropic orientation of the interfacial liqu ids, particular ly those in the first layer (Extended Data Fig. 5a,c,d,e). Despite the rich internal molecular structure and anisot ropic orientation of the interfacial DEC molecules, the geometric center a pproach in the SLS model remains valid. SLS density maps of the DEC at HOPG step ed ge areas (Fig. 2d–f) also closely resemble the experimental results (Fig. 1c–e), both revealing the same layer c rossover and layer termina tion effects. Layer crossover is a manifestation of the mobile nature of liquid molec ules, which rearrange their spatial di stribution to minimiz e the overall free energy. This spatial reconfiguration likely involves a significant entrop ic contribution that disfavors c onformal coverage at the atomic step sites. Near the atomic steps, we further observe a “transition region” (marked between the red dashed lines in Fig. 2d–f) wh ere the liquid layers gradually change their z positions while remaining connected. The transition region occurs due to the competition of the carbon atoms’ contribution to the liquid density through the ETCF, as the substrate at the left vs right side of the step edge favors different liquid layer positions. For larger s tep heights, the impact of the extra substrate layers at the left si de of the step edge becomes more pronounced , causing the transition region to extend over a larger volume (wid er angle between the boundary lines). In addition to the z chan ges of layer positions , we fur ther observe step edge-modulation of liquid density along the x direction (Fig. 2g–i). Near the at omic steps, liquid dens ity profiles exhibit 8 horizontal damped oscillations superposed on top of the average density of the corresponding layer. These oscillations have an average periodicity of 2 𝜋 / 𝛼   4.1  0.2 Å and decay length of 1/ 𝛼   4.5  1.0 Å (see Supplementary Table 2), close to those of the overall vertical density oscillations. We attribute the horizontal oscillation f eatures to interferen ce e ffects originating from the sharp edge termination and propagating alon g the x direction, again demonstrating the SLS principle. To cross validate the AFM-SLS predictions on th e step edge-induced liquid density variations, we perform MD simulations on the same DEC/HOP G atomic step systems (Supplementary Fig. 9b–d and Fig. 2j–o). The MD results reveal nearly identical layer crossover effects as the SLS predictions, including the triangul ar-shaped transition regions with similar size and boundary angles (Fig. 2j–l). At the tran sition region, both the de nsity and molecular or ientation of the DEC change abruptly, as a result of the layer crossover (Supplementary Fig. 4 and Extended Data Fig. 5b,f,g). In addition, we observe si milar interference effects along th e x axis near the step edges (Fig. 2m–o). The average periodicity and decay length of these oscillations are 4.2±0.1 Å and 3.4±0.6 Å, respectively (Supplementary Table 2) , close to the AFM-SLS predictions. These quantitative agreements reveal the accuracy of SLS in predicting liquid de nsity distributions at complex interfaces. Note that horizontal liquid de nsity oscillations have not been observed in experimental force maps. Additionally, the transition regions of the liqui d layers near step edges are wider in the 3D- AFM x-z maps compared to those in AFM-SLS a nd MD results. These differences are likely due to AFM tip convolution effects that result in lo wer lateral resolution (Extended Data Fig. 3 and Supplementary Note 1). Interfacial liquid structure with more complex and larger solid morphologies A key advantage of the SLS model is the ability to quickly pre dict interfacial liquid stru cture adjacent to solid surfaces with ar bitrary morphologies. To further te st this capability, we extend the SLS simulation to another type of commonl y observed heterogeneous surface structure of HOPG, buried step edges, where one or multiple layers of graphene covers an atomic step 33 . We construct a series of these buried step systems with different ramp angles of the top graphene layer (Fig. 3a and Supplementary Fig. 11 ). Plugging in the ato mic positions of the constructed substrates and the already-obtained ETCF f unction between carbon atoms and DEC (Fig. 2a) to Eq. 1, we obtain the SLS-modeled liquid densit y distribution (Fig. 3a). At an angle of 90°, we observe layer crossover effect, similar to the results at exposed atomic steps. As the angle decreases, we find more conformal coverage of the li quid layers. At a ramp angle of 20° or lower, the crossover effect can no longer be observed and the liquid layers fully conform to the height change of the underlying substrate, due to the smoothness of the surface morphology. This conclusion is independent of the number of graphene layers coveri ng the atomic step (Supplementary Fig. 12). Experimentally, the ramp angle of buried step edges has large va riations, and can reach as low as ~0.5° 33,34 . We have imaged a number of these structur es, and observe conformal coverage of DEC above all of them, due to the low ramp a ngles (0.5°–9°, Fig. 3b,e and Supplementary Figs. 3 9 and 13). These results are consistent with the SLS predictions and further illustrate the key role of solid morphology in modulating the interfacial liquid structure. Fig. 3 | Liquid structure at more comple x and larger-scale DEC/HOPG interfaces. a , AFM- SLS liquid density x-z maps near buried step edges of HOPG with different ramp angles. Full 3D atomic structures of the constructed so lids are shown in Suppl ementary Fig. 11. b , A 3D-AFM x- z force map near a burie d step edge of HOPG. c , d , DC mode AFM x-y images of height ( c ) and friction ( d ) of a 360 nm-wide region of HOPG immersed in DEC. The cross-sectional height and friction prof iles extracted from th e dashed line regions are overlaid on top of the x-y images. e , 3D-AFM (AC mode) x-z force map taken at the same region marked by the dashed lines in c and d . f , AFM-SLS x-z liquid density map where the so lid struc ture is constructed according to e . Full 3D atomic structure of the cons tructed solid is shown in Supplem entary Fig. 14. Scale bars: 0.5 nm for x and z ( a ), 2 nm for x and 0.5 nm for z ( b ), 50 nm ( c , d ), 50 nm for x and 0.5 nm for z ( e , f ). The straightforward analytical nature of SLS al so enables direct liquid structure calcu lations of interfacial systems at arbitrar y size scales. This capability can be critical for understanding and predicting systems with a size of 100 nm and larg er, such as whole livi ng cells, microparticle electrodes in batteries, a nd ion exchange membranes 5,35,36 . As a demonstration, we perform 3D- AFM imaging of an HOPG region with a lateral le ngth of 360 nm. As shown in Fig. 3c,d, we 10 observe the co-existence of two he terogeneous sites, one exposed atomic step and the other buried, as confirmed from the x-y surface images of height and friction. Specifically, the exposed step has larger friction than the flat HOPG region, while burie d step has nearly the sa me friction as the flat area (Fig. 3d), consistent with previous reports 33,34 . In addition, the height transition is sharper at the exposed step compared to the buried step (F ig. 3c). 3D-AFM x-z map in the same region reveals horizontal layered liquid structure at the fl at regions, layer cross over at the exposed step edge, and conformal coverage at the buried step (Fig. 3e), all of which are reproduced by SLS simulation (Fig. 3f and Supplem entary Figs. 14 and 15). Such size scales are highly challenging/demanding for al l-atom MD simulations 5,35 ; in contrast, SLS offers a direct solution for the interfacial liquid structure at arbitrarily large scales, with minimal computational cost. Interfacial structure of more complex liquids To further examine the g enerality of the SLS mode l, we investigate the inte rfacial structure of a series of more complex liqui ds, including ethylene carbonate (EC) a nd DEC mixtures with a series of molar ratios (common ly used in lithium-io n batteries 18 ), 1 m (molal) lithium bis(trifluoromethanesulfonyl)imide (LiT FSI) in DEC (a battery electrolyte 37 ) at different electrode potentials, pure water, and 0.01 M K 2 SO 4 aqueous solution at different potentials (Fig. 4 and Extended Data Fig. 6). Potential is controlled using the electroche mical 3D-AFM method developed in our lab 2,21,38,39 . 3D-AFM maps of these liquids at flat HOPG areas reveal similar damped oscillation patterns (Extended Data Fig. 6 and Supplementary Fig. 16), despite the distinct and highly complex compositions of the series of liqui d systems. These oscillation featur es are consistent with previou s 3D-AFM and X-ray scattering measurements at the interface s between multi-component liqu id electrolytes and flat substrates 3,9,21,40, 41 . However, a quantitative unders tanding of these interfacial liquid structures and their exte nsion to heterogeneous solid surfaces have been lacking. Here we propose that Eq. 1 remains valid fo r multi-component liquid systems, where the ETCF can correspond to either one liquid compone nt or the sum of all the components. If individual component’s ETCF is used, Eq. 1 will produce the partial density distribution of th e corresponding component. On the other hand, if the combined ETCF is plugged into Eq. 1, we obtain the total density of the li quid (sum of the partial density of every component). As 3D-AFM force is sensitive to the total density of liquid sp ecies 2,19–21 , the ETCF extracted from 3D-AFM results will roughly correspond to the combined response of all the liquid components. Since damped oscillations have been univers ally observed in the multi-component liquid systems that we have measured, we hypothesize that Eq. 2 is also valid for these interfacial liquids. However, in these complex liquids , 3D-AFM force curves tend to ha ve large fluctuations a t the first liquid peak (Extended Data Fig. 6), likely due to variations in the solvation condition at the end of the AFM probe, as previously reported 21,42 . As a result, it is highly challenging to quantify the first peak density of these complex liquids solely based on 3D -AFM data. To ensure reliable comparison of different liqui d systems, here we set 𝐴   0 in Eq. 2 to simplify the ETCF formula. Supplementary Figs. 17 and 18 rev eal that this s implification does not resu lt in large deviation 11 from the realistic interfacial liqu id density prof ile. Combining the 3D-A FM maps at flat HOPG sites with Eq. 1 and the simplified Eq. 2, we obtain the ETCF for the series of complex liquid systems at carbon surfaces (Supplementary Note 6 and Supplementary Table 3). Since 3D-AFM data and simplified Eq. 2 are used, this ETCF ex traction and subsequent SLS calculation protocol is referred to as “AFM(simp)-SLS”. As shown in Fig. 4a and Extended Data Fig. 7, all ETCFs of the complex liquids exhibit damped oscillation patterns, similar to that of the pure DEC–carbon system. Fig. 4 | Interfacial structure of more complex liquids. a , ETCFs of a series of liquid systems at HOPG surface. ETCF is extracted using the AFM( simp)-SLS protocol (Supplementary Note 6) with parameters listed in Suppl ementary Tables 1 and 3. Mol ecular structures of the liquid components are also shown on the plots. Color code: carbon, green-cyan; oxygen, red: hydrogen, gray; sulfur, yellow; fluorine, brown; nitrog en, blue; lithium, purple; potassium, orange. b – e , x-z maps of molecular density from AFM(simp)-SLS ( b , d ) and x-z maps of 3D-AFM force ( c , e ) of EC/DEC mixture with 1:1 molar ratio ( b , c ) and 0.01 M K 2 SO 4 in water ( d , e ) near HOPG mono- steps. Full 3D atomic structures of th e constructed solids in b , d are shown in Supplementary Fig. 5. In all panels, the solvent-on ly systems are measured at open-circuit, 1 m LiTFSI in DEC at − 0.5 V vs Pt, and 0.01 M K 2 SO 4 at 1 V vs Ag/AgCl. Scale bars: 0.5 nm for x and z ( b , d ), 2 nm for x and 0.5 nm for z ( c , e ). 12 For all the EC/DEC mixtures and LiTFSI/DEC el ectroly te (at different electrod e potentials) at HOPG surface, their ETCFs exhibit a nearly identical periodicity ( 2 𝜋 / 𝛼   4.1  0.2 Å ) and decay length 1/ 𝛼   4.5  1.2 Å (Fig. 4a and Extended Data Fi g. 7). Similarly, the pure water and aqueous solution (at varyi ng electrode potentials ) at HOP G surface also exhibit mostly identical periodicity and decay len gth in ETCF, with 2 𝜋 / 𝛼   3.1  0.2 Å and 1/ 𝛼   3.9  1.5 Å , respectively. This periodicity is the same as the known in termolecular distance of bulk water 20,43,44 . Note that variations in the exact condi tions of each AFM probe have negligible impacts on the extracted 𝛼  and 𝛼  parameters (Supplementary Table 4). For bulk liquids, OZ theory has predicted that the intermolecular TCF of liquid mixtures has the same periodicity and decay length as that of a pure liquid consisting of the dominant component 30 . This effect can be intuit ively understood as the minority species permeating through the “structural framewor k” defi ned by the dominant species. Our results reveal that such bulk liquid behavior manifest s itself in the interfacial density pro file, following the SLS principle (Eq. 1). In EC/DEC mixtures, since pure EC and pure DEC have almost the same periodicity and decay length (~4 Å for both species 40,45,46 ), it is natural to expect their mixture to exhibit the s ame metrics in the ETCF. For 1 m LiTFSI in DEC, the salt to solvent molar rati o is 0.12:1; in 0.01 M K 2 SO 4 in water, the molar ratios are 0.002 :1 (K + vs H 2 O) and 0.001:1 (SO 4 2 − vs H 2 O). Therefore, it is expected that the 𝛼  and 𝛼  parameters are the same as those of the corresponding pure solvent/solid interfaces. While the 𝛼  and 𝛼  parameters tend to converge to the same values as those of the dominant species, other parameters in the ETCF can be mo re sensitive to the liq uid composition changes and/or interfacial inte ractions. In EC/DEC mixtu res, we find that the os cillation amplitude 𝐴  increases as the ratio of EC becomes higher (Fig. 4a and Extended Data Fig. 7). This is likely due to a stronger intermolecular interaction between EC molecules compared to DEC 24,46 ; for the same reason, pure EC is a solid at room temperature while DEC is a liquid. Upon the addition of 1 m LiTFSI in DEC, we observe a broadening of the first liquid layer, which manifests as an overall phase shift in the extracted ETCF (Fig. 4a and Exte nded Data Figs. 6 and 7). This layer broadening / phase shift is possibly due to the a ccumulation of TFSI − in the first lay er, as a result of their strong interaction with HOPG 47 . At different electrode potentials, we observe no significan t changes in the ETCF for both the 1 m LiTFSI in DEC and 0.01 M K 2 SO 4 in water solutions (Fig. 4a and Extended Data Figs. 6 and 7). This is consistent with our previous study of a series of aqueous solutions 21 . In the actual experimental systems, it is likely that the first liq uid density peak will change at sufficiently large electrode potentials du e to charge accumulation/depletion, although the first peak quantification by 3D-AFM is not always straightforward 21,42 . The upper layers are likely always dominated by the solvent regardless of the electrode potential, as the density of charges beyond the first liquid layer is oftentimes much lower than that o f the solvent for the dilute solu tions studied here 21 . We further examine the interfacial structure of these more co mplex liquids at heterogeneous substrate surfaces. We choose HOPG with e xposed mono-step edge as the substrate (Supplementary Fig. 5), which has common well-defi ned structural features that enable the validation of the superposition principle. Using the ETCFs shown in Fig. 4a and Extended Data 13 Fig. 7m–o, we obtain SLS-modele d liquid density maps for EC /DEC mixtures (1:1, 1.82:1, and 2.33:1 mol) and 0.01 M K 2 SO 4 in water (Fig. 4b,d, Extended Data Fig. 8, and Supplementary Fig. 19a,c). In addition, we conduct 3D-AFM measurements of the sa me interfacial systems, with results shown in Fig. 4c,e and Supplementary Fig. 19b,d. For all the systems, SLS predictions match well with the corresponding 3D-AFM maps , both revealing layer crossover effects. EC/DEC mixtures exhibit nearly id entical spatial patterns near the atomic step compared to pure DEC/HOPG, due to their sim ilar ETCFs. In the 0.01 M K 2 SO 4 aqueous solution, liquid layers bend downwards when crossing from the left (above th e upper HOPG surface) to the right (Fig. 4d,e), in contrast to the upward bending observed in pure DEC and EC/DEC mixture sy stems. This difference again manifests the supe rposition principle: a liquid laye r on the left always connects to the nearest-neighbor layer on th e right of the atomic step. Si nce the interlayer distance of aqueous solutions (~3.1 Å) is slight ly smaller than that of the HOPG step hei ght (3.4 Å), layer n (n=1, 2, 3, etc.) on the left need s to bend downwards when connec ting to the nearby layer n+1 on the right; for the same reason, DEC and EC/DEC, with slightly larger interlayer distance (near 4.0 Å) than HOPG step height, exhibit upward bending. This nearby-layer-connect ion effect is also observed for EC/DEC mixtures adjacent to e xposed bi-step edges (S upplementary Fig. 20). In addition, SLS density maps of EC/DEC also show horiz ontal density oscillations next to the step edge (Extended Data Fig. 8b,d), with average period icity of 4.3±0.3 Å and decay length of 3.6±0.1 Å, similar to those of pure DEC. The transition region boundary angles for all EC/DEC and pure DEC at HOPG mono-steps are also simila r (~6°–9°). These simila rities arise from the comparable ETCF periodicity and decay length characteristics of the EC/DEC and DEC systems. In contrast, changes in oscillation amplitude of th e ETCF do not have a direct impact o n the spatial pattern of the liquid distribution. For the aqueous solution system, the horizontal de nsity oscillation differs markedly from that of the EC/DEC systems. Even at the flat HOPG area, we observe in-plane os cillations in the first aqueous layer, with a periodicity of 2.5 Å, co mmensurate with the HOPG substrate (Extended Data Fig. 8e,f). This is because the first peak of the aqueous ETCF lies very close to the origin (2.7 Å, Fig. 4a), resulting in the partial projection of substrate atomic de nsity onto the first liquid peak after applying the SLS model. Although this in-p lane density oscillation of aqueous solutions on HOPG is beyond the sensitivity limit of 3D-AFM, prev ious studies observe d this effect on mineral and molecular crystal surfac es with larger la ttice constants using both 3D-AFM and MD simulation 26,48 , further corroborating our SL S predictions. In upper liquid layers, due to the larger separation from the flat HOPG surfa ce, the in-plane density becomes constant without oscillations. Near the HOPG mono-step site, we find the horizontal profile of the aqueous liquid density to be influenced by both the HOPG basal plane lattice and the edge termination (Extended Data Fig. 8g,h). Namely, the second aqueous la yer on the right of the step edge exhibits damped oscillation with similar periodicity and decay length as the vertical oscillations (~ 3 Å), which resembles th e behavior of (EC/)DEC (Fig. 2g–i and Extended Data Fig. 8b,d); in contrast, the first aqueous layer exhibits an initially larger oscillation period (~3 .2 Å) with large amplitude next to th e step edge, and gradually decays to the 2.5 Å periodicity with smaller amplitude afte r ~4 periods away fro m the step site. 14 Furthermore, the transition region boundary angl e for the aqueous solution (12.2°, marked on Fig. 4d) is larger than that of the (EC/)DEC systems. This is likely due to the smaller mismatch between the liquid layer periodi city (0.31 nm) and the HOPG step height (0.34 nm), which produces a smoother structural transition a nd consequently a la rger boundary angle. These series of results confirm the capabilities of SLS to accurately predict liquid structures at most of the practically relevant solid–liquid interfaces, regardless of the solvent and solute species in the liquid elec trolyte. Conclusions and outlook In summary, we have systematically investigated the structure of interfacial liqu ids at solid surfaces with a large range of complexity and size scale. Inspired by the experimentally observed real-space interfacial liquid conf igurations, we have developed the SLS model to interpre t and predict liquid structures near arbitrary solid surfaces. This model is benchmarked against experimental results and MD simulations, and bridges the co mplexity gap in our existing understanding of interfacial liquids . SLS is broadly appli cable in realistic so lid–liquid interfaces, where the liquids include a diverse set of orga nic and aqueous solvents and electrolytes. To further demonstrate the capabilities and gen e rality of SLS, we perform more SLS modeling of many other systems. Extende d Data Fig. 9 shows the SLS prediction of DEC distribution adjacent to other complex carbon sub strates. A si ngle-atom vacancy at th e HOPG surface produce s only a localized perturbation, manifested as a sl ight downward bending in the first liquid layer, while the surrounding liquid struct ure remains largely unchanged (Extended Data Fig. 9a). In contrast, a flat graphite edge plane induces liqui d layering parallel to th e substrate surface, with pronounced lateral density variati ons (Extended Data Fig. 9b,d). Th e periodicity of these lateral oscillations is 0.68 nm, the same as the substrate (vertical, ABAB-stacked graphene layers). For a rough edge plane, due to the strong substrate di sorder, the adjacen t liquid exhibits significantly disrupted ordering and rich interference patterns (E xtended Data Fig. 9c,e). Furthermore, SLS can be applicable beyond soli d–liquid interfaces, in the general systems of fluids adjacent to rigid bodies. As demonstrations of this broa der scope, we apply SLS to two systems reported in literature, li quid helium surrounding magnesium atoms 49 and 2D molecular gas near immobile sites 50 (Supplementary Note 7). As shown in Extended Data Fig. 10, both the spatial interference patterns and quantitative density profiles predicted by SLS are in good agreement with reported results produced fr om density functional theory (DFT), scanning tunneling microscopy (STM), and ki netic Monte Carlo simulations. We have not encountered any solid–liquid inte rface system where SLS predictions strongly deviate from existing validated resu lts. The exact scope of applicati ons of SLS is a topic for future studies (Supplementary Note 8). We have demonstrated the ease of integra tion of the SLS mode l with 3D-AFM experimental inputs, bridged by the ETCF descriptor. Following similar protocols, we have also successfully integrated MD inputs into SLS, and the resulti ng MD-SLS predictions ar e consistent with AFM- 15 SLS (Supplementary Figs. 17 and 18). The straight forward analytical nature of SLS will also enable its integration with many other experimental and simu lation methods, such as X-ray scattering and DFT, for multi-scale modeling of heterogeneous syst ems where solid–liquid interfaces are key components. Methods Probe and sample preparation Ethylene carbonate (99%), diethy l carbonate (>=99%), and LiTFSI (99.95%) were all purchased from Sigma-Aldrich. Ultrapure deionized water (Milli-Q water) was obtained from the Synergy UV water purification system (M illiporeSigma). K 2 SO 4 (anhydrous, >=99%) was purchased from Acros Organics. HOPG (ZYH grade, 12 mm × 12 mm × 2 mm) wa s purchase d from Struc ture Probe Inc. AFM probes were purchased fr om NanoAndMore (PPP-NCHAuD) and Asylum Research (FS-1500AuD). The tips were made of s ilicon terminated with a native oxide layer. For each AFM probe we used, we measured its detaile d mechanical, optical, and resonance metrics, and provided all of them in Supplementary Table 5. The cantilever was soaked in acetone for ~0.5 h, followed by soaking in isopropan ol overnight (~12 h). Then, it was soaked in ultrapure deionized water for ~3 h before being cleaned by UV ozone (ProCleaner Plus system, BioForce Nanosciences) right before 3D -AFM measurement. After asse mbling the AFM liquid cell, HOPG was quickly cleaved using Scotch tape, and about 120 µL of liquid was dropped on the fresh graphite surface using pipette (Eppendorf North America). The whole AFM liquid cell was purged and then sealed in argon gas. A Pt ring was used as quasi-reference electr ode for all experiments where an electrode potential was a pplied. In aqueous solutions, we further calibrated the Pt quasi- reference against Ag/AgCl, and fo und that 0 V vs Pt is equiva lent to 0.22 V vs Ag/AgCl, as reported in our prior publication 21 . AFM measurements AC, amplitude-modulation mode 3D-AFM measur ements were performe d using a Cypher ES atomic force microscope (Asylum Research, Oxfo rd Instruments). Therma l tuning and cantilever tuning were carried out in air and liquid to de termine th e spring constant, quality factor, inve rse optical lever sensitivity (InvOLS), and resonan ce frequency. Prior to and right after acquiring 3D maps, topographic images of the HOP G surface in the x-y plane were recorded to asse ss the surface cleanliness. 3D-AFM measurements were carried out in clean regions, w here the surface features include flat terraces, exposed atomic steps, and/or buried step edges. During 3D-AFM measurements, the AFM tip scanned a three-di mensional volume above the surface, with a sinusoidal z motion at a rate of 10 Hz. AFM imaging parameters are prov ided in Supplementary Table 6. Typical scan volumes for fl at surfaces, step edges, and la rger-scale structures (e.g., those containing multiple step edges) were 10 × 10 × 3 nm 3 , 40 × 2.5 × 3 nm 3 , and 400 × 2.5 × 3 nm 3 (x × y × z), respectively. Corresponding pixel resolu tions were 64 × 64 × 25,000 for flat surfaces and 256 × 16 × 25,000 for areas containing step edges. 16 3D-AFM data were acquired sequentially as x-z fram es, with y direction as the slow scan axis. For each x-z frame, four datasets we re recorded: trace and retrace (forward and backward lateral motion along x) and approach and retract (dow nward and upward motion along z). These datasets yielded nearly identical maps and were used to verify data reproduc ibility; the same protocol and reproducibility examina tion were reported in our prior publication 2 . At each pixel, the spatial position (x, y, z) as well as th e cantilever oscillation amplitude and phase were reco rded. For lateral deflection me asurements, DC, contact mode scans were perfor med. Two datasets were recorded: trace and retrace. The friction map in Fig. 3d was obtained by subtracting retrace lateral deflection map from trace lateral deflection map (Supplementary Fig. 21). The latera force map in Extended Data Fig. 1h corresponds to the trace scan. Supplementary Fig. 13a shows a lateral force map obtained along the retrace scan. All the AFM measurements reported in this work were conducted at room temperature. 3D-AFM data analysis In AC mode 3D-AFM, each data point consists of x, y, scanner extension (Ext), amplitude, an d phase. Tip–sample distance d was calculated as 𝑑 𝐸 𝑥 𝑡  𝑎 𝑚 𝑝 𝑙 𝑖 𝑡 𝑢 𝑑 𝑒  𝐶 , where C is a constant applied so that d =0 corresponds to the substrate su rface. The measured amplitude and phase data were subsequently converted to force using the me thod initially described by Payam et al 51 . Both the phase–distance and force–distance curves contain two co mponents, a long-range repulsive background and a short-ra nge oscillatory component. On ly the short-range component is directly related to the local liquid d ensity variations 19 . To remove the re pulsive background, we performed double-exponential fit and subtracted the fitted func tion from the original phase– distance and force–distance curves. In this wor k, the first minimum of th e background-subtracted force ( Δ F) curve was used as the origin ( 𝑧 0 ) for plotting the z axis. If not otherwise noted, all the 3D-AFM x-z maps and z curves in this article show the bac kground-subtracted phase or force. Supplementary Fig. 22 shows an example of the force reconstr uction and background subtraction process. SLS calculations All SLS density maps were computed using MATLA B. Graphite substrates with hexagonal lattice (honeycomb structure) and ABAB stacking were constructed, with a carbon–carbon bond length of 0.142 nm and interlayer distance of 0.34 nm. The flat basal planes and exposed step edges used for SLS modeling are illustrated in Fig. 2b a nd Supplementary Figs. 5–7. Unless otherwise specified, exposed step edges in the SLS model adopt an armchair configuration. The choice of step type (armchair vs zigzag) does not have a measurable impact on the resulting liquid density maps, as demonstrated in Supplementary Fig. 23. For buried step edges and large-scale simulations (Fig. 3a,f), the exact substrate geometries are shown in Supplementary Figs. 11 and 14. For each SLS density map, at least ~3000 substrate atoms su rrounding the target inte rfacial liquid region 17 were included and distributed in all spatial directions (x, y, a nd z) (see the example in Fig. 2b), enabling accurate calculation of the liquid dens ity throughout the targ et region using the superposition formula (Eq. 1). E ach x-z liquid density map consis ts of a 500 × 500 pixel grid. At each pixel, the distance between the pixel location and every substrate atom was calculated as 𝑟  . The liquid density variation ( Δ𝜌 󰇛 𝒓  󰇜 / 𝜌  ) at a given pixel was obtained by summing ℎ  󰇛 𝒓  󰇜 over all substrate atoms (at 𝒓  ) following Eq. 1. Repeating this procedure over the entire grid produces the full SLS density map. In typical calculations for fl at and step-edge regions, appr oximately 3,000–10,000 substrate atoms were included. For larger-scale simu lations, more substrate atoms were c onstructed. For example, the graphite substrate shown in Fig. 3f and Supplementary Fi g. 14 contains 131,328 carbon atoms. Molecular dynamics simulations MD simulations of DEC solvent adjacent to HOPG flat area and di fferent atomic steps (mono-step, bi-step and tri-step) were performed using the LAMMPS package 52 . The setups of HOPG/DEC/HOPG sandwich systems are shown in Supplementary Fig. 9. Each MD system includes two HOPG substrates, with each substrate constru cted using multiple layers of graphene sheets (ABAB stacking). The interlayer distance of graphene is 0.335 nm. The flat HOPG system contains 4 graphene layers with a surface ar ea of 5.97 nm × 6.24 nm. Th e HOPG containing step edges has 3 bottom graphene layers with an area of 9.95 nm × 6.24 nm, in addition to 1–3 upper layers with an area of 4.91 nm × 6.24 nm. The st ep edges have an armc hair structure. The separation of the two HOPG/DEC interfaces within each simulated system is at least 12 nm everywhere to ensure DEC molecules in the ch annel center exhibit bulk- like behavior. The force field parameters of HOPG were taken from Cale man et al 53 . DEC molecules were modeled using the Automated Topology Builder 54–56 . Non-electrostatic interac tions, including van der Waals (Lennard–Jones) interactions, were evaluated usi ng a truncated cutoff scheme consistent with th e employed force fields (Caleman et al. and the Automated T opology Builder). The long-range electrostatic interac tions were computed using the particle mesh Ewald (PME) method 57 . A cut- off distance of 1.2 nm was used to compute both the electr ostatic and non-electro static interactions. We applied periodic boundary conditions in the x and y direct ions. The system was non-periodic in the z direction. The simulations were pe rformed under canonical (NVT) ensemble and the temperature was maintained at 300 K by using the Nosè–Hoover thermostat 58,59 . We applied a time step of 1 fs. The simulations we re first equilibrated for 5 ns, followed by a production run for 35 ns. The MD trajectories of atom ic coordination were saved ev ery 1 ps. Each simulation was performed three times using indepen dent initial conf igurations to generate statistically independent trajectories. The resulting obs ervables were averaged to obtain the reported results. Methods for extracting horizo ntal density/force profiles To obtain horizontal density/force prof iles for a segment of x-z map (x 1 ≤ x ≤ x 2 , z 1 ≤ z ≤ z 2 , where z 1 and z 2 can vary at different x), we first extracte d the maxi mum density/force value within the 18 analyzed z interval (z 1 ≤ z ≤ z 2 ) at each x. For MD and SLS results, the maximal liquid density was directly plotted as a function of x (from x 1 to x 2 ). For 3D-AFM force maps, the extracted maximal force values were further smoothened us ing MATLAB’s “smoothdata” function with the “Gaussian” method and a window size of five adjacent points along the x direction. These smoothened maximal force values were plotted as a function of x (from x 1 to x 2 ). Data availability All data are included in the main article, Extended Data, and S upplementary Information. Source data will be provided prior to publication. Code availability SLS codes used to produce Fig. 2b are provided as a Supplement ary Information file. To produce all other SLS maps in this article, the same codes were used, except that the substrate atomic positions (provided in Supplementary Figs. 5–7, 11 and 14) and/or ETCF parameters (specified in Supplementary Tables 1 and 3) were substituted. Acknowledgement Q.A., L.K.S.B., K.S.P., S.Z., F.Z., Y.L. and Y. Z. acknowledge support from the National Science Foundation under Grant No. 2339175, the Beckman Young Investigator Award provided by the Arnold and Mabel Beckman Foundation, and the Sloan Research Fellowship from the Alfred P. Sloan Foundation. Experiments were carried out in pa rt in the Carl R. Woese Ins titute for Genomic Biology and in the Materials Re search Laboratory (MRL) Centra l Research Facilities at the University of Illinois. We are gratefu l to Minjiang Zhu for discussion on the SLS model. Author contributions Q.A. and Y.Z. conceived the overall idea. Q.A. , L.K.S.B., K.S.P., S.Z., and F.Z. conducted the experiments, guided by Y.Z. Q.A. performed th e experimental data an alysis. Q.A. and Y.Z. conceptualized the SLS model and performed calcul ations, with help from Y.L. and K.S.S. H.W. and N.R.A. performed MD simulations. Q.A. and Y.Z. wrote the manuscript, with help from a ll other authors. Y.Z. and N.A. supervised the work. Competing interests The authors declare no competing interests. 19 Additional information Supplementary information: The online version contains supplem entary material (a pdf file for Supplementary Notes, Figures, and Tables , and a zip file for SLS source code). Correspondence and requests for materials should be addressed to Yingjie Zhang. 20 Extended Data Fig. 1 | AFM x-y images of HOPG su rface immersed in liquid. Liquid systems and electrode potentials ar e: DEC at open-circuit ( a – d ), 0.01 M K 2 SO 4 in water at 0 V vs Ag/AgCl ( e ), EC/DEC 1:1 mol ( f ) and 2.33:1 mol ( g ) mixtures at open-circuit, and 1 m LiTFSI in DEC at 1 V vs Pt ( h , i ). a – g are AC mode surface height maps, while h is a DC mode lateral force image. The height profiles at positions marked by the arrow (next to the left panels in a – g ) are extracted and shown as the right panels in a – g . Across each step edge, the hei ght change usually occurs over a width of several nanometers, instead of an infini tely sharp ju mp, due to the finite size of the AFM tip (see Extended Data Fig. 3). Scale bars: 125 nm ( a ), 12.5 nm ( b ), 10 nm ( c , d ), 50 nm ( e , f ), 37.5 nm ( g ), 2 nm ( h ), 5 nm − 1 ( i ). h resolves the atomic lattice of the HOPG surface, with the Fourier transform shown in i . From h , i , the HOPG lattice constant is extracted as ~0.23 nm, consisten t with our previous observations and the known standard value 21,38,39 . 21 22 Extended Data Fig . 2 | Ga l l er y o f 3 D- AF M x -z ph as e maps of pure DEC and EC/DEC near exposed atomic steps of HOPG. The measurement date an d liquid system are specified on top of each map. All results are measured at open-circuit. These maps are ob tained from a total of 8 different experiments with multiple types of AFM tips, measured at three types of exposed atomi c steps (mono-step, bi-step, and tri-step). Th e measured tip parame ters are provided in Supplementary Table 5. The gray lines at th e bottom of each map correspond to the assigned position and structure of HOPG s ubstrate. Low phase (yellow) regions in the map roughly correspond to the liquid layers. Layer crossover eff ect is universally obser ved in all the maps. 23 24 Extended Data Fig. 3 | Ti p convolution effect. Simulated interfacial system: DEC/HOPG. a , Intrinsic liquid density map calculate d by AFM-SL S. The inset on the left shows a simu lated AFM tip that senses an averag e signal surrounding the tip ap ex (tip radius: 5 nm). b , Force vs x along a connected liquid layer acro ss a step edge, obtained from an experimentally measured fo rce map of DEC/HOPG mono-step (same as Su pplementary Fig. 4d). The transi tion region is ~1.9 nm wide. c , Horizontal maximal density curve ex tracted from the red dashed box in a . The transition region is ~0.5 nm wide, narrower than that obtained from the correspondi ng experimental system shown in b . d – f , Simulated sensitivity maps surrounding AFM tips with a radius of 5 nm ( d ), 10 nm ( e ), and 20 nm ( f ). Regions with higher sensitivity have higher weighted contribution to the corresponding convoluted density map. g – i , Simulated tip-convoluted x-z density maps after dilation operation using the intrinsic density map in a and the corresponding sensitivity maps in d – f , respectively. j – l , Horizontal maximal densit y profiles extracted from the red dashed boxes in g – i , respectively. The transition region width increases w ith tip radius. m , Schematic illustratin g the formation of an ina ccessibl e area near the step edge due to the finite size o f the tip. n , o , Simulated tip-convoluted x-z density maps illustra ting the evo lution of the inaccessible area as the step structure change s, including bi-step ( n ) and tri-step ( o ) (compare to the mono-step in h ). Full 3D atomic structures of the constructed solids in a , g – i , n , o are shown in Supplementary Figs. 5 and 7. Scale bars: 1 nm for x and 0.5 nm for z ( a , g – i , n , o ), 0.5 nm for x and 0.25 nm for z ( d – f ). Color bars for g – i , n , o are the same as that shown on the left of a . 25 Extended Data Fig. 4 | Reproducibility of 3D-AFM measurements. Results are summarized from 3D-AFM force curves of DEC at HOPG flat areas at open-circuit. a , Force curves from 5 different experiments (measuremen t dates labeled) using two type s of cantilevers (PPP-NCHAuD and FS-1500AuD). Each colored curve is an average of at least 150 single force curves. Black curve is the average of all colo red curves. Si milar force curves are observed in all five measurements. For each measurement, we fur ther extract the character istic parameters of the force curves following part of the AFM( simp)-SLS protocol, with result s summarized in Supplementary Table 4. Similar parameters are obtained regardless of the measuremen t date or cantilever type, 26 further confirming the quantitative re producibility of our 3D-AFM results. b , An example to illustrate the automated protocol for determinin g layer spacings. Local minima, iden tified using MATLAB’s “islocalmin” function, are marked by red dots. The layer spacings of interfacial liquid are calculated as the distances between adjacent minima. 𝑑  is the distance between the first and second layer, 𝑑  the distance between the second and the third layer, and so forth. c , DEC layer spacings extracted from measuremen ts at different dates. The spacings remain consistent among all experiments conducted over multiple years. E ach data point (average value and standard deviation) represents calculation fr om over 128,000 individual phase curves. d , DEC layer spacings (average value and standard deviation) m easured using different types of cantilevers. The spacings are found to be independent of the canti lever type. Over 600,000 individual phase curves are obtained and used for analysis using PPP- NCHAuD tips, and more than 300,000 curves using FS-1500AuD tips. 27 Extended Data Fig. 5 | MD simulated molecu lar orientation at the flat and mono-step transition regions. System: DEC/HOPG at open-circuit. a , MD simulated x-z density map of DEC at flat HOPG surface. Molecular orientation is calcu lated at the four high density regions marked by red boxes. b , MD simulated x-z density map of DEC near an HOPG mono-step edge (same as Fig. 2j). Molecular orientation is calculated at the 8 boxed areas across the layer crossover region. Scale bars ( a , b ): 0.5 nm for both x and z. c , Schematic of a DEC molecule and its simplified form. Green-cyan: carbon; red: oxygen; gray: hydrogen. The mo lecule is simplified as an elongated plate for the purpose of marking the a ngular orientations. The blue arro w marks the direction connecting the two oxygen atoms with single carbon–oxygen bond. The red triangle mark s the plane formed by the three oxygen atoms. d , e , Simulated orientation d istribu tion from the corresponding liquid layers marked in a . The angle is betw een the normal direction of the substrate and ( d ) the blue lin e connecting two oxygen atoms with single carbon–oxygen bond or ( e ) the normal direction of the red plane of three oxygen atoms in DEC. In th e first liquid layer, DE C molecules are mostly parallel to the substrate. In layer 2, the molecu les are largely randomly di stributed with a slight tendency toward the flat configuration. In layer 3 and above, the molecular orientation is almost fully random. f , g , Orientation distribution from the eight spots marked in b . Molecules gradually evolve from a typical layer #1 orient ation to a typical layer #2 orientat ion from left to right as layer crossover occurs. 28 Extended Data Fig. 6 | 3D-AFM force curves of more co mplex liquids at flat HOPG substrates. Red curves are individual force vs z curves . Black curves are averaged fro m red curves. The solvent-only systems are measured at open -circuit. For liquid electrolytes, the mark ed potentials are referenced vs Pt for LiT FSI in DEC, and vs Ag/AgCl for 0.01 M K 2 SO 4 in water. In our previous work, we conducted 3D-AFM measurements of 0.01 M K 2 SO 4 in water on HOPG at a series of negative potentials vs Ag/AgCl, and observed no measur able differences in the force curves 21 . 29 Extended Data Fig. 7 | Additional ETCF quant ification of more complex liquids at HOPG surface. ETCF is extracted using the AFM(simp)-SLS protocol. a – l , Extracted ETCF parameters for pure DEC and EC/DEC mixtures at open-circuit ( a – d ), 1 m LiTFSI in DEC at different potentials ( e – h ), and water at open-circuit and aque ous solution at diffe rent potentials ( i – l ). The plotted parameters include oscillation amplitude 𝐴  ( a , e , i ), periodicity 2 𝜋 / 𝛼  ( b , f , j ), decay length 1/ 𝛼  ( c , g , k ), and phase factor 𝜃  ( d , h , l ). Each plotted average va lue and error bar (standard deviation) are obtained from a set of e xtracte d ETCF functions, where each ETCF is obtained from fitting of a single force curve. Statistics for th e (EC/)DEC system are ex tracted from more than one million force curves, whereas statistics for 1 m LiTFSI in DEC and for the aqueous system are each extracted from more than 82,000 force curves. m – o , ETCF plots for EC/DEC with 1.82:1 molar ratio at open-circ uit ( m ), 1 m LiTFSI in DEC at different potentials ( n ) and 0.01 M K 2 SO 4 in water at 0 V vs Ag/AgCl ( o ) constructed using the average parameters from a – l . 30 Extended Data Fig. 8 | Horizontal density profil es of complex liquids at HOPG surfaces (flat or mono-atomic step region). a , c , e , g , AFM(simp)-SLS x-z density maps for EC/DEC mixture (1:1 mol) at mono-step ( a ), EC/DEC mixture (1.82:1 mol) at mono-step ( c ), 0.01 M K 2 SO 4 in water at flat area ( e ), and 0.01 M K 2 SO 4 in water at mono-step regions ( g ). a , c are simulated for open-circuit and e , g are simulated for 1 V vs Ag/AgCl. Panels a and g are identical to Fig. 4b,d, respectively. Scale bars: 0.5 nm for x and z ( a , c , e , g ). b , d , f , h , Maximal density vs x profiles extracted from areas marked by the dashed boxes in a , c , e , g , respectively. Red density curves correspond to the red dashed boxes, while gray cu rves to gray dashed boxes. Yellow and cyan dashed curves are fits for red and gray curves, respec tively. Fitting para meters are provided in Supplementary Table 2. Full 3D atomic struct ures of the constructed solids are shown in Supplementary Fig. 5 ( a , c , g ) and Fig. 2b ( e ). 31 Extended Data Fig. 9 | AFM-SLS maps of DEC adjacent to other complex carbon substrates. a , x-z liquid density map near HOP G containing a surface atomic va cancy. The white circle shows the position of the missing atom. b , c , x-z density maps of D EC adjacent to a flat edge plane ( b ) and a rough edge plane ( c ) of graphite. d , e , Schematics of the substrate configurations used in b , c , respectively, where the top substrate surface has zi gzag edg e structure (Supplementary Fig. 6). Scale bars: 0.5 nm for x and z ( a – c ). 32 Extended Data Fig. 10 | Application of the SLS model to reported systems beyond solid– liquid interfaces. Panels a and b are DFT calculation results from Ref. 49 . Copyright 2008 American Physical Society. a , Map of helium density surrounding two Mg ato ms in a liquid helium droplet. b , Helium density vs z at x=0 (blue solid line) and helium density vs x at z=0 (red dashed line). c , Helium density map produced from our SLS mode l with the same Mg atomic pair as in a . Scale bars: 1 nm for both x and z. d , SLS-modeled liquid helium density profile vs z at x=0 (blue solid line) and density vs x at z=0 (red dashed line). e , The interference pattern (left) and reflected sigmoid background (right) of c . They are multiplied to mi mic the helium droplet boundary condition as reported in Ref. 49 , producing c . f – k , Results for a 2D molecular gas system, highly mobile molecules of fluorina ted copper phthalocyanine (F 16 CuPc) on a Si(111)/Tl-(1×1) surface. f – i panels correspond to Fig. 4a,e,b,f of Matvija et al 50 , respectively. Copyr ight 2017 American Chemical Society. f , g , Experimental STM images of mobile F 16 CuPc molecules surrounding one ( f ) and two ( g ) fixed F 16 CuPc molecules. h , i , Molecular density maps fr om a kinetic Monte Carlo simulation, corresponding to the same systems shown in f , g . j , k , Molecular density maps calculated from our SLS model, corresponding to the same systems shown in h , i . Scale bars ( h – k ): 3 nm. In both the liquid helium and 2D molecular gas systems, SLS results are in agreement with existing results in the corresponding literatur e. Details of SLS calculation procedures for both systems are provided in Supplementary Note 7. 33 References 1. Ball, P. Water as an Active Constituent in Cell Biology. Chem. Rev. 108 , 74–108 (2008). 2. Ai, Q. et al. Nucleation at solid–liq uid interfaces is accompanied by the reconfiguration of electrical double layers. Proc. Natl. Acad. Sci. 122 , e2421635122 (2025). 3. Toney, M. F. et al. Voltage-dependent ordering of water mo lecules a t an electrode–electro lyte interface. Nature 368 , 444–446 (1994). 4. Becker, M. et al. Multiscale modeling of aqueous electric double layers. Chem. Rev. 124 , 1–26 (2024). 5. Wang, A., Kadam, S., Li, H., Shi, S. & Qi, Y. Review on modeling of the anode solid electrolyte interphase (SEI) for lithium-ion batte ries. npj Comput. Mater. 4 , 15 (2018). 6. Ilgen, A. G. et al. Bridging molecular-scale interfacial scien ce with continuum -scale models. Nat. Commun. 15 , 5326 (2024). 7. Schreier, M., Kenis, P., Che, F. & Hall, A. S. Trends in electrocata lysis: The microenv ironment moves to center stage. ACS Energy Lett. 8 , 3935–3940 (2023). 8. Zhang, Q. et al. Atomic dynamics of elec trified solid–liquid inte rfaces in liquid-cell TEM. Nature 630 , 643–647 (2024). 9. Martin-Jimenez, D., Chacon, E., Tarazona, P. & Garcia, R. Atomically resolved three- dimensional structures of electroly te aqueous solutions ne ar a solid surface. Nat. Commun. 7 , 12164 (2016). 10. Yurtsever, A. et al. Molecular insights on the crystalline cellu lose-water interfaces via three-dimensional atomic force microscopy. Sci. Adv. 8 , eabq0160 (2022). 11. Groß, A. & Sakong, S. Ab initio simula tions of water/metal interfaces. Chem. Rev. 122 , 10746–10776 (2022). 12. Stern, O. Zur theorie der el ektrolytischen doppelschicht. Z. Für Elektrochem. Angew. Phys. Chem. 30 , 508–516 (1924). 13. Limaye, A., Suvlu, D. & Will ard, A. P. Water molecules mute the dependence of the double-layer potential profile on ionic strength. Faraday Discuss. 249 , 267–288 (2024). 14. Zhang, M. et al. How charge frustration causes ion ordering and microphase separation at surfaces. Nat. Mater. 25 , 487–494 (2026). 15. Li, C.-Y. et al. In situ probing electrifie d interfacial water structur es at atomically flat surfaces. Nat. Mater. 18 , 697–701 (2019). 16. Wang, T. et al. Enhancing oxygen reduction electrocatal ysis by tuning in terfacial hydrogen bonds. Nat. Catal. 4 , 753–762 (2021). 17. Velasco-Velez, J.-J. et al. The structure of interfacial wa ter on gold electrodes studied by x-ray absorption spectroscopy. Science 346 , 831–834 (2014). 18. Xu, K. Electrolytes and interpha ses in Li-ion batteries and beyond. Chem. Rev. 114 , 11503– 11618 (2014). 34 19. Ai, Q., Bonagiri, L. K. S., Farokh Payam, A ., Aluru, N. R. & Zhang, Y. Toward quantitative interpretation of 3D ato mic force microscopy at solid–liqu id interfaces. J. Phys. Chem. C 129 , 5273–5286 (2025). 20. Li, Y., Ai, Q., Bonagiri, L. K. S. & Zhang, Y. Coexistence of two type s of liquid structures at platinum–water interfaces. J. Phys. Chem. C 129 , 21486–21495 (2025). 21. Bonagiri, L. K. S. et al. Probing the molecular s tructure at graphite–water interfaces by correlating 3D-AFM and SHINERS. Nat. Commun. 17 , 2230 (2026). 22. Hernández-Muñoz, J. et al. Subnanometer interfacial forces in three-dimensional atomic force microscopy: Water and oc tane near a mica surface. J. Phys. Chem. C 124 , 26296–26303 (2020). 23. Hernández-Muñoz, J., Chacón, E. & Tarazona, P. Density functional analysis of atomic force microscopy in a dense fluid. J. Chem. Phys. 151 , 034701 (2019). 24. García, G. et al. Structure of alkylcarbonate + n-alkane mixed fluids. J. Phys. Chem. B 118 , 11310–11322 (2014). 25. Kim, J., Zhao, F., Bonagiri, L. K. S., Ai, Q. & Zhang, Y. Electrical double layers modulate the growth of solid–electrolyte interphases. Chem. Mater. 36 , 9156–9166 (2024). 26. Umeda, K. et al. Atomic-resolution thr ee-dimensional hydration structures on a heteroge neously char ged surfa ce. Nat. Commun. 8 , 2111 (2017). 27. Miyata, K. et al. High-speed three-dimension al scanni ng force microscopy visualization of subnanoscale hydration structures on dissolving calcite step edges. Nano Lett. 24 , 10842–10849 (2024). 28. Siretanu, I., Lin, S. R. van & Mugele, F. I on adsorption and hydration forces: a comparison of crystalline mica vs. amorphous silica surfaces. Faraday Discuss. 246 , 274–295 (2023). 29. Henderson, J. R. & Sabeur, Z. A. Liquid ‐ state integral equations at high density: On the mathematical origin of infin ite ‐ range oscillatory solutions. J. Chem. Phys. 97 , 6750–6758 (1992). 30. Evans, R., Leote de Carvalho, R. J. F., Henders on, J. R. & Hoyle, D. C. Asymptotic decay of correlations in liqu i ds and their mixtures. J. Chem. Phys. 100 , 591–603 (1994). 31. Ryu, C. W. & Egami, T. Medi um-range atomic correlation in simple liquids. I. Distinction from short-range order. Phys. Rev. E 104 , 064109 (2021). 32. Ornstein, L. S. & Zernike, F. Accidental deviations of density and opalescence at the critical point of a single substance. Proc. R. Neth. Acad. Arts Sci. 17 , 793–806 (1914). 33. Lee, H., Lee, H.-B.-R., Kwon, S., Salmeron, M. & Park, J. Y. Internal and external atomic steps in graphite exhibit dramatic ally di fferent physical and chemical properties. ACS Nano 9 , 3814–3819 (2015). 34. Chen, Z., Khajeh, A., Martini, A. & Kim, S. H. Identifying phy sical and chemical contributions to friction: A comparative study of chemically inert and active graphene step edges. ACS Appl. Mater. Interfaces 12 , 30007–30015 (2020). 35. Durrant, J. D. et al. Mesoscale all-atom influenza viru s simulations suggest new substrate binding mechanism. ACS Cent. Sci. 6 , 189–196 (2020). 35 36. Foglia, F. et al. Disentangling water, ion and polymer dynamics in an anion exchange membrane. Nat. Mater. 21 , 555–563 (2022). 37. Xing, L. & Xu, K. Charging imid e anions for batteries. Nat. Chem. 17 , 1976–1976 (2025). 38. Zhou, S., Panse, K. S., Motevaselian, M. H. , Aluru, N. R. & Zhang, Y. Three-dimensional molecular mapping of ionic liquids at electrified interfaces. ACS Nano 14 , 17515–17523 (2020). 39. Bonagiri, L. K. S. et al. Real-space charge density prof iling of electrode–electrolyte interfaces with angstrom depth resolution. ACS Nano 16 , 19594–19604 (2022). 40. Steinrück, H.-G. et al. The nanoscale structure of the el ectrolyte–metal oxide interface. Energy Environ. Sci. 11 , 594–602 (2018). 41. Mezger, M. et al. Molecular layering of fluorinated io nic liquids at a charged sapphire (0001) surface. Science 322 , 424–428 (2008). 42. Fukuma, T. & Garcia, R. Atomic- and mo lecular-resolution mapp ing of solid–liquid interfaces by 3D at omic force microscopy. ACS Nano 12 , 11785–11797 (2018). 43. Narten, A. H., Thiessen, W. E. & Blum, L. At om pair distribution f unctions of liquid water at 25°C from neutron diffraction. Science 217 , 1033–1034 (1982). 44. Garcia, R. Interfacial liq uid water on graphite, gra phene, and 2D materials. ACS Nano 17 , 51–69 (2023). 45. Westphal, L . et al. Temperature-resolved crystal st ructure of ethylene carbonate. J. Phys. Chem. C 129 , 11546–11553 (2025). 46. Magd ă u, I.-B. et al. Machine learning force fields for molecular liquids: Eth ylene Carbonate/Ethyl Methyl Ca rbonate binary solvent. npj Comput. Mater. 9 , 146 (2023). 47. Vatamanu, J., Borodin, O. & Smith, G. D. Molecular insights into the potential and temperature dependences of the differential capacitan ce of a room-temperatu re ionic liquid at graphite electrodes. J. Am. Chem. Soc. 132 , 14825–14833 (2010). 48. Yurtsever, A. et al. Interplay between β -Chitin nanocrystal supra molecular archite cture and water structuring: Insights from three-dime nsional atomic force mi croscopy measurements and molecular dynamics simulations. J. Am. Chem. Soc. 147 , 39134–39150 (2025). 49. Hernando, A., Barranco, M., Mayol, R., Pi, M. & Ancilotto, F. Dens ity functional theory of the structure of magnesium-doped helium nanodroplets. Phys. Rev. B 78 , 184515 (2008). 50. Matvija, P., Rozbo ř il, F., Sobotík, P., Ošt’ádal, I. & Ko cán, P. Pair correlation function of a 2D molecular gas directly visual ized by scanning tunneling microscopy. J. Phys. Chem. Lett. 8 , 4268–4272 (2017). 51. Payam, A. F., Martin-Jimenez, D. & Garc ia, R. Force reconstruction from tapping mode force microscopy experiments. Nanotechnology 26 , 185706 (2015). 52. Plimpton, S. Fast parallel algorith ms for short-range mo lecular dynamics. J. Comput. Phys. 117 , 1–19 (1995). 53. Caleman, C. et al. Force field benchmark of organi c liquids: Density, enthalpy of vaporization, heat capacities, surface tension, iso thermal compressibil ity, volumetric expansion coefficient, and dielectric constant. J. Chem. Theory Comput. 8 , 61–74 (2012). 36 54. Malde, A. K. et al. An automated force field topol ogy builder (ATB) and repository: Version 1.0. J. Chem. Theory Comput. 7 , 4026–4037 (2011). 55. Koziara, K. B., Stroet, M., Malde, A. K. & Mark, A. E. Testing and validation of the automated topology builder (ATB ) version 2.0: Prediction of hydration free enthalpies. J. Comput. Aided Mol. Des. 28 , 221–233 (2014). 56. Stroet, M. et al. Automated topology builder versi on 3.0: Prediction of solvation free enthalpies in water and hexane. J. Chem. Theory Comput. 14 , 5834–5845 (2018). 57. Darden, T., York, D. & Pedersen, L. Particle mesh Ewald: An N ⋅ log(N) method for Ewald sums in large systems. J. Chem. Phys. 98 , 10089–10092 (1993). 58. Nosé, S. A unified formulation of the cons tant temperature molecu lar dynamics methods. J. Chem. Phys. 81 , 511–519 (1984). 59. Hoover, W. G. Canonical dynamics: E quilibrium phase-space distribu tions. Phys. Rev. A 31 , 1695–1697 (1985).