Local grafting heterogeneities control water intrusion and extrusion in nanopores

Hydrophobic nanoporous materials can be intruded by water only by exerting an external action, typically increasing pressure. For some materials, water extrudes when the pressure is lowered again. Controlling intrusion/extrusion hysteresis is central in a number of technological applications, including materials for energy applications and for high performance liquid chromatography, and experimental techniques, as liquid porosimetry, but is still far from being understood. In this work, we consider water intrusion and extrusion in common mesoporous materials grafted with hydrophobic chains, showing that the macroscopic properties of the system are significantly affected by subnanometric heterogeneities in the grafting. For example, intrusion and extrusion pressures can vary more than 20 MPa depending on the chain length and density of the grafting. Coarse-grained molecular dynamics simulations reveal that local changes of radius and contact angle produced by grafting heterogeneities can pin the interface during intrusion or facilitate bubble nucleation in extrusion. These unprecedented microscopic insights can directly impact the design of energy materials and chromatography columns, as well as the interpretation of porosimetry results.


Introduction
Hydrophobic nanoporous materials combine "water fear" with confined spaces, which together strongly repel water from the pores.In these conditions, water intrudes the nanopores only if an external action is exerted, e.g., under pressure; once the intrusion process is completed, vapor bubbles must nucleate in order to give rise to the opposite phenomenon of extrusion, which typically occurs at lower pressures.These processes can be exploited in a number of applications, which range all the way between energy storage and damping 1,2 , depending on the hysteresis between intrusion and extrusion pressures.There are other contexts in which water extrusion is undesired: in high performance liquid chromatography (HPLC), dewetting is at the origin of retention losses, 3 1 arXiv:2305.15250v3[cond-mat.soft]14 Jun 2024 which makes the stationary phase unavailable to analytes, while in nanopore sensing it might cause undesirable noise in the ion current 4 .Finally, intrusion isotherms are used in water porosimetry 5,6 to measure the pore size distribution, interpreting them in terms of the classical Laplace's equation 7 .In all these cases, a connection between the microscopic characteristics of the pores (radius, contact angle, local heterogeneities, etc.) and the macroscopic properties of the system (intrusion and extrusion pressures, hysteresis, dissipated energy, etc.) is highly needed to design specific applications, control undesired phenomena, and interpret experimental results.
Nevertheless, the reference theoretical framework relies on the definition of average physical parameters describing the geometry of the pores and the hydrophobicity of the grafted surface, such as average radius and contact angle.A fundamental question still awaits an answer: how do the microscopic heterogeneities in the grafted pore surface affect global properties of interest for applications?
Here, we employ coarse-grained molecular dynamics simulations to reveal the microscopic mechanisms of water intrusion and extrusion in silanized mesoporous materials.Interestingly, even though the substrate is the same, small differences in the grafting -different chain lengths and grafting densities -cause significant changes in the intrusion and extrusion processes, affecting crucial technological parameters, such as the intrusion and extrusion pressures, the shape of the pressure-volume isotherms, and the dissipated energy.Simulations allowed us to reveal the microscopic origin of such variability, which is rooted in the local properties arising from chemical and topographical heterogeneity within the pores.For intrusion, the grafted layer's thickness is insufficient to justify the increase in the intrusion pressure via Laplace's equation; only by including local changes in the pore radius and in the hydrophobicity is it possible to account quantitatively for its variability.For extrusion, on the other hand, simulations show that local accumulations of hydrophobic material and constrictions act as nucleation seeds for vapor bubbles.These effects were not accessible by previous simulation efforts on idealized hydrophobic nanopores without explicit functionalization 18,23 .
The present results shed light on the multiscale nature of intrusion and extrusion phenomena, revealing how molecular details control the macroscopic behavior of nanoporous materials, determining their technological applicability.In particular, physical and chemical heterogeneities in the functionalization of mesoporous materials are shown to play a crucial role, being able to change the experimental intrusion/extrusion pressures by up to 60 MPa and to almost double the dissipated energy.Controlling grafting can thus be favorably exploited to rationally design materials for various applications including energy storage, vibration damping, and liquid chromatography.In addition, microscopic insights facilitate the interpretation of liquid porosimetry measurements, revealing significant deviations from the expected macroscopic behavior due to grafting heterogeneities.

Coarse-grained model of grafted mesopores
Our goal is to build a sufficiently detailed model of a mesopore grafted with different hydrophobic chains yet with a low enough computational cost to simulate experimentally relevant time and length scales.To this aim, we constructed a coarse-grained (CG) model of a silica nanopore a) Figure 1: Setup of simulation system and TEM micrographs of MCM-41 material.a) Cross-section of the simulated system, which consists of a cylindrical pore confined between two water reservoirs whose pressure is controlled by movable pistons.The inset below shows the front view of the pore with three graftings, obtained by hydrophobic chains of different lengths (C8, C12, and C18).The molecular structures of these silanes are shown along with their respective coarse-grained scheme (Si in yellow, C in black, and H in white).TEM micrographs of MCM-41 from Liu et al. 24 are reported in panel b (section perpendicular to the pores' axis) and c (parallel).functionalized with a hydrophobic silane grafting. 16The simulated systems have geometrical similarities with the nanoporous material MCM-41, made of essentially independent, parallel cylinders with diameters ranging from 2 to 11 nm 16,25 .We used the most recent MARTINI force field (MARTINI 3 26 ) in which, on average, two to four heavy atoms and associated hydrogens are mapped into one CG bead.This framework enables the computation of systems larger than those accessible to all-atom models and for longer times, with an overall speedup of at least two orders of magnitude 27 , while maintaining sufficient chemical and spatial resolution.
We modeled the silica matrix with a slightly hydrophilic cubic fcc lattice, from which a cylindrical pore of diameter 5.2 nm and length 20 nm was excavated.The pore was immersed in water, and the liquid pressure was controlled mechanically by the force acting on two pistons far away from the pore 28 , allowing to perform in silico intrusion/extrusion (I/E) cycles analogous to experimental ones.Periodic boundary conditions were applied in the directions orthogonal to the pore axis.
To functionalize the pore, we used three alkyl (CH 3 )C m H 2m+1 linear chains having different lengths: C8 (m = 8), C12 (m = 12), and C18 (m = 18).The chains were grafted to the surface at one end at random locations and oriented perpendicularly to the normal vector of the local internal surface (Supplementary Figure S1).Each bead on the solid surface was bonded to at most one silane.For each of the three chain lengths (C8, C12, C18), we produced systems with a variety of grafting densities spanning the range 0.6 to 1.2 groups per square nanometer (gps nm −2 ).In the following section, we characterize the degree of hydrophobicity achieved by different graftings.We remark that, albeit some features of MARTINI water do not precisely match those of actual water, the phenomenology investigated here is related to solvophobic phenomena, which are expected to be independent of the particular liquid 29 .Figure 2: Contact angle θ Y computed on a flat surface with different grafting densities and chain lengths.Each point results from an average of 5 independent simulations of cylindrical droplets deposited on the surface; error bars denote the standard deviation.A linear interpolation between data points (solid lines) is shown as a guide to the eye.The dashed lines are the average SASA of the silica surface, computed using a test particle with the same van der Waals radius as the smallest water bead.Blue, green, and orange represent C8, C12, and C18, respectively.

Hydrophobicity of grafted surfaces
In order to assess the effect of the grafting on the hydrophobicity of the pore, we computed the contact angle θ Y of cylindrical water droplets deposited on a flat silica surface grafted with different chain lengths (C8, C12, and C18) within a broad range of grafting densities (0 to 3.0 gps nm −2 ).The average contact angles are reported in Figure 2 as a function of the grafting density.
In accordance with the experimental wisdom of silica being hydrophilic at ambient conditions 30 , we modeled the substrate to obtain a Young contact angle θ Y ≈ 75 • .Although the actual contact angle of silica has been reported within a wide range of values, its dependency with the percentage of accessible silanol groups is well known 31 .For a 50% fraction of accessible silanol groups, which would correspond to an MCM-41 material 32 , the contact angle is in the vicinity of 75 • .For completeness, in the Supplementary section "Effect of silica hydrophobicity" we discuss I/E when more hydrophilic substrates are considered.Figure 2 shows that, for grafting densities below 0.5 gps nm −2 , the interaction of water with the hydrophilic substrate dominates, showing a progressive increase of θ Y with the grafting density; effects of the chain length are hardly discernible in this regime.At intermediate grafting densities, between 0.5 and 1.5 gps nm −2 , longer chains result in higher θ Y for the same grafting density.Finally, at large grafting densities, the contact angle reaches a plateau at θ Y ≈ 117 • , showing that the effect of additional chains on the hydrophobicity saturates.The grafted pores we used in the following sections have grafting densities 0.6, 0.8, 1.0, and 1.2 gps nm −2 , which correspond to somewhat higher effective densities because of curvature.
We further computed the Solvent Accessible Surface Area (SASA) of the silica substrate for the different graftings (dashed lines in Figure 2).This analysis clearly shows that the grafting diminishes the contact between water and the hydrophilic substrate, causing a progressive increase in the hydrophobicity of the grafted surface.Longer chains provide better coverage for the same grafting density, which corresponds to lower SASA; this affects the contact angle more prominently at intermediate grafting densities.As the grafting density further increases, the chains become more packed and the hydrophilic patches vanish altogether.Accordingly, the grafted layer approaches a complete coverage and the hydrophobicity becomes independent of the chain length. 2 gps nm −2 for the three chain lengths; insets show the dry and wet states of the pore.These cycles mimic the typical I/E experiments that can be performed on similar materials, 16 except that, in simulations, negative pressures can be reached.e) Experimental I/E curves 33 for an MCM-41 matrix grafted with C8 at two different grafting densities: 1.24 (dashed line) and 1.35 gps nm −2 (solid line).We remark that experimental materials have a much smaller radius than in our simulations, R C8−1.24 = 1.75 nm and R C8−1.35 = 2.1 nm before grafting, and that the compressibility of the system was subtracted by the authors from the I/E curves.

In silico intrusion/extrusion experiments
We performed in silico I/E cycles for all grafted pores by progressively increasing the water pressure until complete filling of the pore was achieved and subsequently decreasing it until emptying.These cycles are analogous to typical I/E experiments performed on similar materials, 16 except that, in simulations, negative pressures can be reached, and the cycles are necessarily faster, although this is partially compensated by the shorter pores.Figure 3 illustrates the I/E curves obtained for different grafting densities (panels a-c) and chain lengths (d).The most apparent finding is that there are significant qualitative and quantitative changes in the intrusion and extrusion curves for a relatively narrow change in chain lengths (C8-C18) and grafting densities (0.6-1.2 gps nm −2 ).Intrusion and extrusion pressures vary by > 20 MPa and the shapes can account for rather abrupt I/E processes, as expected for monodisperse cylindrical pores, to rather progressive ones (C18, 1.2 gps nm −2 ).
Figure 3a-c shows that higher grafting densities shift both the intrusion and the extrusion curves to higher pressures, for all considered chain lengths.Panel d reveals that, for the same grafting density, the I/E curves move to higher pressures as the chain length is increased.Our data show the same qualitative behavior -growth of the I/E pressures with grafting density-as the experimental data (Figure 3e) for an MCM-41 matrix with different grafting densities and radii 33 .It is also worth noticing that our simulations capture the change in slope with grafting density observed in such experiments.
In the attempt to rationalize the effect of grafting on I/E, we note that both the chain length and the grafting density are expected to affect the size of the pore available to the water.Indeed, Figure 4c shows that the mean thickness of the grafting layer grows with chain length and Figure 4: On the left, relation between intrusion/extrusion pressures and dissipated energy with ⟨R⟩ (mean layer thickness).On the right, ⟨R⟩ (mean layer thickness) variation with pressure.a) I/E pressures, defined as the inflection points of the respective curves, against mean layer thickness and ⟨R⟩.b) Dissipated energies for each system, as computed by multiplying the area between extrusion and intrusion curves times the typical pore density of MCM-41 34 (330 mm 3 /g ).c) Mean layer thickness and mean radius of the grafted pores as a function of pressure.Notice that the pressure increases until 30 MPa, shown with a vertical line, then starts decreasing.⟨R⟩ is calculated as R pore subtracted of the mean layer thickness; to calculate distances, we use the point position of the beads, in consistency with that used for computing the contact angles in Figure 2. The cross/circular markers represent dry/wet configurations, delimited by a linear fit to the I/E pressures.The color legend refers to both figures.More information on how to compute the mean radius, mean layer thickness, and dissipated energy can be found in the "Supplementary Method" section in SI.
grafting density, modifying the mean radius of the grafted pore ⟨R⟩ between ca.1.4 and 2 nm.The variance in ⟨R⟩ increases with grafting density, with the maximum attained for C18.The compression of the grafted layer upon intrusion is seen to be negligible: the softness of the ligands does not seem to play a critical role for the considered systems.In fact, Supplementary Figure S6 shows that the mean radius increases/decreases no more than 0.02 nm during intrusion/extrusion.In addition, performing I/E cycles with artificially immobilized silanes does not alter the results, as shown in Supplementary Figure S8.
The difference in ⟨R⟩ across systems suggests that the amount of hydrophobic material added inside the hydrophilic channel is key to understanding the I/E mechanisms.Figure 4a indeed highlights that the intrusion and extrusion pressures P int/ext grow monotonically with ⟨R⟩, in agreement with experimental data 5 .Similarly, the dissipated energy E d per mass of porous material depends on ⟨R⟩, Figure 4b.One can notice that E d for the lowest grafting density and shortest chain is ca.1.5 times the energy of the highest grafting density and longest chain, consistent with previous results which reported an increase of the dissipated energy with pore radius 18 .
Overall, these results suggest that the type and quality of grafting can significantly change the macroscopic behavior of hydrophobic mesoporous materials, affecting their technological applicability.For instance, not all the considered materials display extrusion at positive pressures, which allows for a single-time use only (e.g., for "bumpers" 1 ) .In the following, we dissect the microscopic origin of such unexpected variability, showing that the effect of ⟨R⟩ considered in Figure 4 just accounts for a general, qualitative trend, while local variations in pore radius and hydrophobicity due to grafting heterogeneities are the quantitative key to understand intrusion and extrusion.S1.Bare nanopores details can be seen in Supplementary Figure S9.a) Intrusion curves normalized according to Laplace's law (1) using the mean radii and contact angles.The points show the inflection point of the curve, which we define as P int .The vertical line shows where P int is exactly equal to Laplace's equation prediction.In b) the curves are rescaled using the minimum radii and mean contact angle, while in c) we used mean radii and maximum contact angles.Details on the calculation of the maximum contact angle can be found in the Supplementary Information.Panel d) shows the rescaling with both the local parameters.i, ii, and iii labels data refer to C8, C12, and C18 chains, respectively.

Microscopic origin of intrusion
Previous results showed that the characteristics of grafting, i.e., chain length and grafting density, affect both the overall degree of hydrophobicity (θ Y , Figure 2) and the mean radius (⟨R⟩, Figure 4) of the pores.Are these average quantities sufficient to explain the changes observed in the intrusion pressure across different systems?To answer this question, we employ the macroscopic Laplace equation, which describes the pressure P int after which the meniscus depins from the cavity mouth, thus giving rise to intrusion: 35 where γ lv is the liquid/vapor surface tension and R the pore radius.According to eq. ( 1), the intrusion pressure should increase with higher contact angles and lower radii, which is in agreement with the general trends shown in Figure 3.For a quantitative comparison, in Figure 5a.i-iiiwe report the intrusion curves rescaled using eq.( 1) adopting the average contact angles θ Y , the mean radius ⟨R⟩ of all the replicas for each pore configuration (different grafting density and chain length), and γ lv = 21 ± 1 mN m −1 which was calculated using the test-area method. 36f only the mean values of θ Y and ⟨R⟩ were relevant for the intrusion process, the rescaled curves relative to all different graftings would overlap.However, a true collapse of the curves is obtained only for the bare nanopores made of different materials and with different pore radii, while the grafted nanopores get closer to each other but do not overlap.The systems that exhibit higher deviations from the predicted behavior are those with lower amounts of grafting; this trend is unexpected since lower hydrophobicity and higher radius should result in lower intrusion pressures.Therefore, Figure 5a reveals that the average radius and contact angle do not fully account for the variability of the I/E curves.
The intrusion curves in the plot are the results of an average between cycles obtained starting from different realizations of the same system, i.e., with the same macroscopic parameters (grafting density and chain length) but different microscopic configurations, with random distributions of the ligands.The intrusion curves in Figures 6a and b are relative to individual intrusion events of C18-0.8 and C18-1.2 systems, respectively.These curves present larger slopes than those of bare channels, as shown in Supplementary Figure S10.The steps in the intrusion curves reflect a microscopic stick-and-slip mechanism that water undergoes during intrusion.The free-energy profile for water intrusion (Figure S11 in Supplementary "Free energy calculations" section) indeed reveals the presence of intermediate free energy barriers between the dry and wet states.Microscopically, the pore surface is characterized by alternations of constrictions, namely local reductions of the pore radius due to grafted silanes, and hydrophilic patches, namely (sub)nanometer areas of exposed bare silica: when water meets hydrophilic patches, it is locally attracted and thus slips.When water meets constrictions, it is locally pinned and thus sticks.These microscopic heterogeneities can be quantified in terms of local variations of pore radius and contact angle.
In order to verify which of these two local parameters is the most relevant, we rescaled the intrusion curves using the minimum pore radius and the average contact angle or the maximum contact angle and the average radius of each individual trajectory, as shown in Figure 5b and c, respectively.Figure 5b shows that the rescaling with the minimum radius of the nanopore and the mean contact angle generally improves the agreement of the pores with higher grafting density with the bare ones.For C18, a collapse of all curves is obtained except C18-0.6.Overall, this analysis shows that, at higher grafting densities, intrusion is affected by topographical heterogeneities, which are capable of pinning the meniscus within the pore.
On the other hand, rescaling intrusion curves with the maximum contact angle and the mean pore radius (Figure 5c) improves the curves' overlap with the bare pores for most systems except C18-1.0 and 1.2, accounting for the importance of chemical heterogeneities (the appearance of hydrophilic patches), which are indeed expected at lower grafting densities and for shorter chains, due to the higher SASA of the substrate (Figure 2).More quantitatively, as shown in Figure 5c, the number of hydrophilic patches for C18-0.6 and C18-0.8 is more than twice than in C18-1.2.On the other hand the typical size of hydrophilic patches, including the largest one, is larger for C18-0.6 resulting in a total amount of hydrophilic silica exposed which is significantly higher than in the other two systems.This analysis confirms the presence of two kinds of heterogeneities (topographical and chemical), but the distinction is not necessarily sharp, as hydrophobic chains can affect both the local radius and contact angle leading to interface pinning; this scenario seems to apply to intermediate cases, which are indeed rescaled correctly using either local quantity.Finally, if both local quantities are taken into account, all curves tend to collapse into a master curve (Figure 5d).This analysis clarifies that the intrusion mechanism is severely dependent on local heterogeneities which need to be taken into account to predict changes to the intrusion process.For more compact graftings, topographical patchiness tends to pin the meniscus by creating bottlenecks inside the channel (inset of Figure 6b, displaying C18-1.2).The systems with lower amounts of hydrophobic grafting are more affected by chemical patchiness, with the interplay of hydrophobic and hydrophilic areas controlling water intrusion (inset of Figure 6b, showing C18-0.8).To complete wetting in grafted pores, at variance with bare pores, water must overcome the most severe pinning conditions, which can be topographical (minimal radius), chemical (most hydrophobic patch), or combinations thereof 37 .

Microscopic origin of extrusion
Although it can be described within the same framework, 18 the conditions at which extrusion occurs are not described by Laplace's equation (1).Actually, extrusion originates in the nucleation of a vapor phase from the confined liquid (cavitation), which involves the nontrivial competition of bulk, surface, 35,38 and line [16][17][18] terms.Thus, thermal fluctuations can stochastically lead to the emptying of an individual pore by overcoming a free energy barrier ∆Ω † .
By looking at individual extrusion trajectories, we noticed that extrusion starts at specific locations within the pore characterized by the smallest radius (Figure 7a and b, Supplementary 2) grafting densities.The match between high contact angle and small radii favors nucleation, as shown by the critical bubbles in the reported MD snapshots.c) Position of the nucleation bubble along the pore (z nuc ) as a function of the position of the minimum radius along the pore (z min ).d-f) Extrusion pressure P ext for individual MD simulations as a function of the minimal radius R min (filled circles) and of the average radius ⟨R⟩ (filled x).The predictions of the nucleation theory (2) are reported in dashed lines for the average contact angles θ Y of each grafting density and chain length combination.The shaded areas are used to demarcate the limiting curves for each chain length: blue for C8, green for C12 and red for C18.
Figure S12-S13).Figure 7c shows a fair correlation between the position of the nucleation bubble along the channel and the position of the minimum radius.This region usually also corresponds to the most hydrophobic region, because at those locations more hydrophobic material is accumulated.If along the pore there are multiple constrictions, multiple bubbles can nucleate at the same time, as shown in Supplementary Figure S14.Together, constrictions and enhanced hydrophobicity, which are maximum at the highest grafting densities, conspire in favoring the formation of a cavitation nucleus.
More formally, expressing the nucleation free energy as a sum of volume, surface, and line terms, one can build a simple classical nucleation theory for the formation of a bubble in a cylindrical pore 16,39 .In particular, the nucleation barrier ∆Ω † can be expressed as where V * , A * lv , A * sv , and l * slv are the volume, liquid-vapor area, solid-vapor area, and triple line length of the critical bubble, respectively; τ is an effective line tension that takes into account nanoscale effects, see, e.g., 18,40,41 .We used the estimates of Lefevre et al. 16 for the geometrical quantities.For the thermodynamic parameters, we kept a constant value of ∆Ω † = 10 k B T, which is compatible with the simulated times, and considered a range of τ close to the ones previously reported in literature 18,42,43 , −10 pN < τ < −5 pN).
Figures 7d-f compare P ext of individual simulations with the nucleation theory (2).We use a definition of P ext which is slightly different from Figure 4: we take the point in which the water volume starts a steep decrease instead of using the inflection point of the extrusion curve.In this way, we select with more precision the pressure at which a vapor phase nucleation starts, which is the process described by the theory.We further consider the minimal radius R min of each pore (local quantity) and the average contact angle for each chain length/grafting density pair.With these quantities, the theory seems to capture well the trend of P ext , showing that higher values are connected to larger θ Y .In contrast, the average value ⟨R⟩ predicts extrusion pressures which are far off, especially for C18.
It is worth remarking that both τ and ∆Ω † are expected to vary for each individual system, reflecting the variability in the grafting.Hence, a more precise quantitative agreement could be expected only if one had a careful measurement of these quantities and a nucleation theory that takes into account the non-ideal geometry of the pore.
In summary, for extrusion, heterogeneities in the grafting gives rise to variability in the local radius and hydrophobicity that determine where the most probable sites for nucleation are located.The extrusion pressures can be estimated using a nucleation theory, which however needs to be informed with local quantities that are hard to characterize τ , ∆Ω † , and R min , in addition to θ Y of the specific grafting.
Our Molecular Dynamics results, demonstrating the pivotal importance of the molecular features of the grafting on the I/E process, were found to be in agreement with recent numerical simulations 44 which were performed by means of a minimal lattice model, along with specialized experiments 45 which were conducted on straight pores functionalized with extremely controlled and reproducible grafting patterns.

Discussion
Before discussing the impact of our work on three different technological applications, it is worth discussing in more detail the possibility to achieve a quantitative comparison between the simulated I/E cycles and the experimental ones.There are three main sources of discrepancies between simulated and experimental data which need to be considered: pore geometry, forcefield used in simulations and time and length scale differences between experiments and simulations.
The experimental samples obviously exhibit a more heterogeneous pore geometry, including pore interconnections, structural defects, and polydispersity in lengths and radii, which are not captured by the simple cylindrical pore geometry used in our simulations.These effects, whose understanding is still in its infancy 46 , will add up to the effect of heterogeneities in the grafting which is the focus of this work.
The I/E process is ruled by interfacial energies.Atomistic forcefields can reasonably approximate experimental surface and line tensions.Coarse grained models, and the MARTINI forcefield in particular, have the tendency to underestimate them.While our coarse grained approach is expected to capture qualitatively the physico-chemical driving forces for I/E processes, a comparison between experimental and simulated intrusion pressures requires the application of a correction factor.As the liquid-vapor surface tension of water in the MARTINI forcefield is calculated as γ sim = 21 mN m −1 while the experimental value is γ exp = 71.7 mN m −1 47 , simulated pressures (as well as dissipated energies) should be multiplied by a factor of γ exp /γ sim = 3.41 to be compared to the experimental ones.This rescaling for instance predicts that the maximum variability of P int and P ext due to heterogeneities in the experimental case is 60 MPa.In Supplementary Figure S15, we show a quantitative agreement between simulated and experimental 5 intrusion and extrusion pressures as a function of grafting density, once this rescaling is taken into account.
Typical lengths of experimental MCM-41 pores, which are the archetype of our simulations, can be up to several µm long, 48 while typical experimental time scales are of the order of seconds 34 .On the other hand, the simulated pore length in this work are 20 nm, while constant-pressure simulation windows had a duration of 10 ns.We have addressed in detail the quantification of the maximum discrepancy between experimental and simulated intrusion and extrusion pressures that may emerge from these different time and length scales.The details can be found in the Supplementary section ""Bridging simulations and experiments: effect of time and length scales on experimental and simulated I/E pressures", but we remark that for P int the expected discrepancy is tens of kPa while for P ext is of the order of few MPa, which we consider fully satisfying, given the overall goal of elucidating the physics of I/E in grafted nanopores.
We now discuss the results from the perspective of three typical technological applications of water intrusion/extrusion in grafted silica nanopores, namely energy materials, high performance liquid chromatography, and porosimetry.
I/E of liquids within hydrophobic pores is at the heart of a broad family of energy materials: 21 depending on the I/E hysteresis, such systems can be used for vibration damping or energy storage applications.In both cases, the design for the specific operative conditions hinges on the ability to control P int , P ext , and their range.Until now, the complexity of the pores and their nanometric size hindered a satisfying characterization of the I/E behavior of promising materials, like suitably coated porous silica.Herein, our coarse-grained approach allowed us to derive important structureproperty trends for energy applications, showing that grafting details cannot be ignored to control I/E cycles precisely.
As shown in Supplementary Table S1, in simulations I/E pressures can vary more than 20 MPa for the same bare pore radius, depending on the length and density of the functionalizing chains (this corresponds to 68 MPa in experiments).For instance, using the longest ligands (C18), one can double the grafting density to more than double the intrusion pressure, moving at the same time the extrusion from negative pressures to almost 14 MPa (48 MPa in experiments).Our results also suggest that small grafting densities with short ligands are to be preferred to increase the dissipated energy because these coatings preserve some hydrophilic patches.However, in these conditions, often no extrusion is possible at positive pressures, and thus the material becomes suitable for single-use energy dissipation (so-called bumpers 1 ).
In summary, the information gathered via coarse-grained simulations (summarized in Table S1) is vital to learn how to tune the operative conditions and the stored and dissipated energy by proper nanopore functionalization.In future studies, our coarse-grained approach will be extended to include more intricate pore geometries, accounting for the full complexity of porous materials for energy applications.
Results from our simulation campaign are of interest also for High-Performance Liquid Chromatography (HPLC), an analytical technique used to separate, identify, and quantify components in a mixture.HPLC is common in several industrial fields, such as pharmaceutical, biotech, and food processing and safety.Reversed-phase liquid chromatography (RPLC) is a type of HPLC that uses a non-polar stationary phase made of a hydrophobic material, such as C8 and C18-functionalized silica, and a polar mobile phase to analyze mixtures containing apolar and moderately polar compounds.The analysis of very polar compounds may require 100% aqueous mobile phases, which are associated with irreversible retention losses, as explained by Gritti 3 .Both experiments 49 and Monte Carlo simulations 50 revealed that dewetting of water from the hydrophobic stationary phase (also wrongly referred to as "phase collapse") is the reason behind retention loss, which makes part of the RPLC column unavailable to retain analytes.Several parameters can affect the dewetting kinetics: temperature, pore size, and water contact angle.However, the internal microstructure of the material, namely the pore size distribution, the surface coverage and chemistry, the pore connectivity 22,46 , and the presence of dissolved gases 51 affect the retention loss as well.In practice, to minimize the loss of retention, the hydrostatic pressure in the column should never go below the extrusion pressure; ideally P ext < 0. Large intrusion pressures, on the other hand, are indicative of the difficulty of filling the column for the first time or of regenerating one that underwent dewetting.
The present results show that the grafting can significantly affect both the intrusion and the extrusion pressures (Supplementary Table S1) and should be carefully considered in the design of HPLC columns.For instance, longer chain lengths and higher grafting densities are shown to exacerbate the dewetting problem and make column regeneration more difficult.The grafting details are known to impact the selectivity towards certain analytes 52 and the grafting densities are related to the retention of analytes in chromatography 53 .For this reason, our coarse-grained approach can be a predictive tool to guide the challenging research of a balance between minimal retention losses and maximum selectivity.
Finally, we proceed to discuss the importance of the present results for the field of porosimetry.Indeed intrusion experiments that are conceptually similar to those performed here in silico are routinely used in porosimetry to infer geometrical information about the pores themselves.Depending on the intruding liquid, we distinguish between water and mercury porosimetry, with the latter being a standard technique due to the very high surface tension of mercury and its typical high contact angles on most surfaces, even non-functionalized ones. 16,54n porosimetry, it is generally assumed that intrusion is a deterministic process that follows Laplace's equation ( 1); more precisely, it is assumed that all the pores of a given radius are filled as soon as the pressure exceeds their corresponding Laplace's pressure.Based on this assumption, the range of pressures at which intrusion takes place is interpreted as resulting from an underlying distribution of pore radii.As a result, in porosimetry, the experimental intrusion curve for a given material is used to provide information about the pore size distribution.The limitation of this procedure is that it does not take into account the variability of the radius and of the contact angle along each pore, which we just proved to be important to shape the intrusion curves.In our simulations, the actual distribution of local radii along the pore is a directly measurable quantity; indeed Figure 8 show that grafted pores can have a rather broad pore size distribution with the variance, which is related to the pore roughness, increasing with grafting density.In Figure 8 we compare the actual pore size distribution with "in silico porosimetry", i.e., the radius distribution obtained by using Laplace's law to interpret the simulated intrusion curves in Figure 3. Results show that pore size distributions obtained via porosimetry vastly differ from the actual geometrical distribution.
It is particularly striking that the most probable values can differ by almost 1.5 nm for the lowest grafting densities and the variance can be either smaller or larger than the actual one depending on the chain length.In the C18 case, at 0.6 and 1.2 gps nm −2 (Figure 8c), porosimetry results could imply a bimodal distribution of the radii, which is not reflected in the actual radii population.Finally, long tails corresponding to large radii apparent in the in silico porosimetry are not reflected in the actual pores, where the distribution has a well defined range.Overall, our analysis shows that porosimetry, in particular water intrusion porosimetry on functionalized mesopores, can significantly underestimate the mean pore size and be poorly quantitative on the overall distribution.
In principle, also extrusion curves could be used to calculate the radius distribution through a nucleation theory; however, the stochastic nature of extrusion, together with local quantities that need to be taken into account, and the difficult independent estimates of line tension and nucleation kinetics, makes this approach more involved.On the other hand, Figure 5 shows that Laplace's equation describes quantitatively all the intrusion curves when local quantities (R, θ Y ) are used; this result suggests that a detailed characterization of the grafting, although experimentally challenging, could significantly improve the accuracy of porosimetry.

Conclusions
We performed coarse-grained MD simulations of water intrusion/extrusion in silica nanopores with different graftings.Our simulations show that grafting plays an active role in the intrusion and extrusion processes, even when the bare nanopores have an almost ideal geometry.Depending on the chain length and the grafting density, intrusion was observed in the approximate range 3 < P int < 23 MPa while extrusion was observed for −10 < P ext < 12 MPa.
As a consequence, related quantities, such as dissipated energies, were also observed to be strongly dependent on the grafting details.A distinctive behavior was observed for intrusion and extrusion in the pore grafted with the longest hydrophobic chains at the highest grafting density considered in this work (C18 with 1.2 gps nm −2 ): very pronounced slopes were observed in the intrusion/extrusion branch of the PV isotherms.
Using data collected in the simulation campaign, we were able to propose a microscopic physical explanation of the mechanisms ruling the intrusion and extrusion processes.For intrusion, heterogeneities change the local pore radius and contact angle affecting the intrusion pressure according to a microscopically informed Laplace's law, deviating from the average Laplace's law which describes pinning at the cavity mouth only.On the other hand, extrusion is essentially a nucleation process, which is not expected to follow Laplace's law.Indeed, we find that extrusion starts by nucleation at sites within the pore where local features favor it: local constrictions were observed to play a major role in this process because they both decrease the volume of the critical bubble and they increase the local hydrophobicity.The interplay between hydrophobic and hydrophilic regions emerges as a promising aspect to be tailored for specific responses.
We discussed the specific importance of simple grafting heterogeneities on a number of relevant technologies, which include nanoporous systems for energy applications, HPLC stationary phases, and standard porosimetry protocols.In all applications of hydrophobized mesoporous materials, grafting emerges as a lead character whose fine control could provide an optimized functional design and improve characterization methods.Table S1: Summary of the main quantities characterizing the grafted pores considered in this simulation work, including the mean radius ⟨R⟩, the minimal one R min , the mean contact angle θ Y , and the maximum one θ Y,max , the intrusion pressure P int , the extrusion one P ext , and the dissipated energy E d .We remind that pressures and energies should be multiplied by 3.41 to translate these quantities to experimental values.
in equal number percentage.We performed 1 µs CG simulations in the NVT ensemble at 298 K to compute the contact angle of water on silica surfaces with different functionalizations.For each system, we performed five runs with different starting configurations.In every run, the contact angle was calculated from the last 500 ns of the trajectory.To compute the contact angle, we used an in-house Python script that performs a circular cap fit to the isodensity contour of the border of the droplet with the vacuum interface.The same simulation and analysis protocols were used to compute the contact angle on the bare hydrophobic surfaces, made of C2 and C4 beads.The Solvent Accessible Surface Area (SASA) was calculated using the GROMACS tool sasa.We used the smallest bead radius in our water model (0.191 nm) as the radius of the probe.For each grafting density and chain length, the results are averaged over all the realizations of the pore.

Intrusion/extrusion cycles
To compute intrusion/extrusion cycles, two movable pistons were used.The pistons consist of a thin fcc lattice with dimensions (20x20x2.5)nm 3 .The length of the nanopore is 20nm, and the radius is 5.2nm previous to any grafting.The starting dimensions of the two water cubic boxes are (20x20x20)nm 3 .We used an in-house Python script, exploiting the MDAnalysis library, ?? to analyze our simulations.Firstly, we computed the density of the water reservoirs during the I/E cycle to monitor the pressure applied by the pistons on the system.We calculated the densities of two water boxes located at a fixed distance from the pistons at each pressure step.Then we compared the results with those obtained for a water box with the same bead size composition, simulated in the NPT ensemble using a Parrinello-Rahman barostat.As shown in Fig. S2, the three sets of values match, proving the proper functioning of our pulling protocol.The same script was used to compute the filling of the pore, calculating the average number of intruded water beads at each pressure for every run.We then derived the filling fraction by dividing the number of intruded beads by the maximum filling of each run.The I/E cycles for all the single realizations of the pores are reported in Fig. S3.Supplementary Figure S2: a) Comparison between water densities in the two reservoirs (reservoirs 1 and 2) at different pressures during the intrusion and the density in the NPT water box (reference).b) Comparison betweeen the nominal pressure (in green) and the pressures derived from water densities in the two reservoirs.Supplementary Figure S12: Local radius of the grafted pore and local contact angle θ Y as a function of the axial coordinate for two representative systems with low (C18-0.6)and high (C18-1.2) grafting densities and 4 different realizations of the pore.The MD snapshots show the water inside the pore when the pressure is nearby the extrusion pressure.While, at low grafting densities, the nucleation site is not always evident, at high grafting densities the coordinates where the vapor phase nucleation occurs match the ones of the minimum radius for every single run.

Figure 3 :
Figure 3: In silico intrusion/extrusion cycles.The effect of different grafting densities is shown for C8 (a), C12 (b), and C18 (c).d) I/E cycles at grafting density 1.2 gps nm −2 for the three chain lengths; insets show the dry and wet states of the pore.These cycles mimic the typical I/Eexperiments that can be performed on similar materials,16 except that, in simulations, negative pressures can be reached.e) Experimental I/E curves33 for an MCM-41 matrix grafted with C8 at two different grafting densities: 1.24 (dashed line) and 1.35 gps nm −2 (solid line).We remark that experimental materials have a much smaller radius than in our simulations, R C8−1.24 = 1.75 nm and R C8−1.35 = 2.1 nm before grafting, and that the compressibility of the system was subtracted by the authors from the I/E curves.

Figure 5 :
Figure 5: Intrusion curves for all the grafted and for three bare pores, i.e., without any grafting and with different substrate hydrophobicities and radii (Bare 1: θ Y = 100 • and R = 1.3 nm, Bare 2: θ Y = 100 • and R = 1.8 nm, and Bare 3: θ Y = 105 • and R = 1.3 nm).A summary of the grafted pore measured quantities can be found in Supplementary TableS1.Bare nanopores details can be seen in Supplementary FigureS9.a) Intrusion curves normalized according to Laplace's law (1) using the mean radii and contact angles.The points show the inflection point of the curve, which we define as P int .The vertical line shows where P int is exactly equal to Laplace's equation prediction.In b) the curves are rescaled using the minimum radii and mean contact angle, while in c) we used mean radii and maximum contact angles.Details on the calculation of the maximum contact angle can be found in the Supplementary Information.Panel d) shows the rescaling with both the local parameters.i, ii, and iii labels data refer to C8, C12, and C18 chains, respectively.

Figure 6 :
Figure 6: Details of pore heterogeneities: patches and constrictions affect water filling rates.Panels a) and b) show the water filling fraction versus time relative to single runs for C18-0.8 and C18-1.2,respectively.Vertical lines show the times at which pressure steps are applied to the system.The insets show water pinning inside the channel.In panel a) the hydrophilic patches are highlighted in pink; panel b) shows a local constriction, in which the radius of the pore has a minimum due to the presence of the silanes (in red).c) Quantitative analysis of the hydrophilic patches in C18 systems with different graftings: the first panel shows the area of the largest hydrophilic patch, the second the total number of patches, and the third panel reports the total area of bare silica exposed.The error bars refer to the standard deviation of the values for the different realizations of each system.

Figure 7 :
Figure 7: Local variations of radius and contact angle determine preferential nucleation sites.Nucleation theory can better describe simulation results when minimum radii are provided.The local radius of the grafted pore and local contact angle θ Y is shown as a function of the axial coordinate for two representative systems with a) low (C18-0.6)and b) high (C18-1.2) grafting densities.The match between high contact angle and small radii favors nucleation, as shown by the critical bubbles in the reported MD snapshots.c) Position of the nucleation bubble along the pore (z nuc ) as a function of the position of the minimum radius along the pore (z min ).d-f) Extrusion pressure P ext for individual MD simulations as a function of the minimal radius R min (filled circles) and of the average radius ⟨R⟩ (filled x).The predictions of the nucleation theory (2) are reported in dashed lines for the average contact angles θ Y of each grafting density and chain length combination.The shaded areas are used to demarcate the limiting curves for each chain length: blue for C8, green for C12 and red for C18.

Figure 8 :
Figure 8: Comparison of the geometrical radii and radii derived from Laplace's law.Radii distributions are shown for different chain lengths, C8 (a), C12 (b), and (c) C18, as obtained from simulated porosimetry analysis based on Laplace's law (filled regions with dashed lines), compared to the actual distribution of radii as obtained by analyzing simulations (solid lines).Vertical dotted lines denote the bare pore radius.
Tailored surface graftings with complex chemistry, polymerization, and polydispersity could enable new functions, pushing towards increased efforts in the development of functionalization protocols, characterization techniques, as well as predictive theoretical tools.experimental consultancy.G.A. T.A., and G.R. supervised research.All authors critically revised the manuscript.