Liquid intrusion in and extrusion from non-wettable nanopores for technological applications

In this article, we review some recent theoretical results about intrusion and extrusion of non-wetting liquids in and out of cavities of nanotextured surfaces and nanoporous materials. Nanoscale confinement allows these processes to happen at conditions which significantly differ from bulk phase coexistence. In particular, the pressure at which a liquid penetrates in and exits from cavities is of interest for many technological applications such as energy storage, dissipation, and conversion, materials with negative compressibility, ion channels, liquid chromatography, and more. Notwithstanding its technological interest, intrusion/extrusion processes are difficult to understand and control solely via experiments: the missing step is often a simple theory capable of providing a microscopic interpretation of the results, e.g., of liquid porosimetry or other techniques used in the field, especially in the case of complex nanoporous media. In this context, simulations can help shedding light on the relation between the morphology of pores, the chemical composition of the solids and liquids, and the thermodynamics and kinetics of intrusion and extrusion. Indeed, the intrusion/extrusion kinetics is determined by the presence of free energy barriers and special approaches, the so-called rare event techniques, must be used to study these processes. Usually, rare event techniques are employed to investigate processes occurring in relatively simple molecular systems, while intrusion/extrusion concerns the collective dynamics of hundreds to thousands of degrees of freedom, the molecules of a liquid entering in or exiting from a cavity, which, from the methodological point of view, is itself a challenge.


Introduction
In this review, we discuss recent progress in the theory and in the simulation of phase changes of liquids confined in non-wettable cavities. In particular, we will focus on the processes of intrusion in and extrusion from non-wettable textured surfaces and porous solids, i.e., confined wetting/evaporation, and discuss their technological relevance. The conditions at which these processes occur sensitively depend on material properties and external conditions, whose understanding is crucial for the technological applicability of novel classes of materials. Here, we present a simple and yet powerful continuum theory of the intrusion and extrusion processes and recent results of advanced atomistic simulations. This combination allows to bridge scales which range from microns down to angstroms, where the continuum models fail, providing a set of theoretical tools to interpret experiments and design materia e-mail: simone.meloni@unife.it (corresponding author) als for specific applications. Theoretical achievements are described in the context of parallel experimental advances-one of the authors of this review being an experimentalist working in the field. Non-wettable porous or nanostructured materials show a non-trivial behaviour when immersed in a liquid: as the name implies, at ambient conditions, the cavities are occupied by gaseous phases. Only upon application of sufficient pressure, or other perturbing stresses (e.g., heat), capillary forces are overcome and cavities are intruded by the liquid. Thus, the nomenclature used in the field denotes as "non-wettable" those porous solids and textured surfaces that are not intruded by a specific liquid at ambient conditions. Similarly, sometimes the focus in on the liquid, which is denoted as "nonwetting" with respect to a given porous solid. Once the cavities are wet, the opposite process of evaporation of the confined liquid, or extrusion, occurs when the pressure (or other external conditions) is again decreased (changed); depending on the characteristics of confine-ment, in particular in nanoscale cavities, such pressure may be larger than the ambient one.
Intrusion/extrusion of liquids in porous systems is important for many technological applications, among which the separation of liquids [1,2], liquid-phase chromatography [3,4], energy damping and storage [5,6], porosimetry for the characterization of the porous systems [7][8][9], biological and bioinspired channels [10][11][12][13], and many more. Intrusion/extrusion is important also for textured surfaces, i.e., surfaces presenting roughness that enhance their intrinsic lyophobicity, the property to repel a liquid (see below). The macroscopic properties of these surfaces depend on the state of the liquid deposited on them, that is, whether it wets the bottom of surface roughness or remains suspended over their asperities. Thus, understanding in quantitative terms the transition from one state to the other and how this depends on the characteristics of the textures allow designing surfaces with properties tailored to the application. Atomic-scale simulations might be a necessary computational tool to engineer materials whose macroscopic behaviour depends on nanoscale processes, here intrusion and extrusion. When the characteristic length at play is small, i.e., when cavities of porous systems are of the same order of magnitude of the liquid/gas molecules, the very concept of continuum fluid passing through an aperture breaks down. This happens, for example, in the case of zeolites, metal-organic frameworks, covalent-organic framework as well as for pores of biological relevance, presenting characteristic size of pores and apertures of the order of few molecular radii. In these cases, macroscopic concepts, such as surface tension and interface energy, are insufficient at describing the energetics and dynamics of a liquid intruding/extruding a pore.
Atomic-scale simulations of intrusion/extrusion are necessary also to investigate elusive microscopic effects such as line tension, the force associated with a variation of the length of the solid-liquid-gas contact line, the Tolman length, determining the dependence of the surface tension on the curvature of the surface, orientation of liquid molecules at the interface with the gas and solid, layering of liquids in contact with solid walls, disjoining pressure and, as discovered recently, accumulation of gasses and other solutes at the solid/liquid interface, which depends on the chemical characteristics of the solid, liquid and solute [11,12].
In this work, we discuss intrusion/extrusion focusing on some fundamental aspects of these processes and on their applications. Instead of providing a broad summary of the field, and show how simulations can help its investigation, we chose to select few paradigmatic cases that allow to establish the ground for a wider range of applications. With this in mind, we decided to follow the results of a limited number of articles that can guide the reader from the basic understanding of intrusion/extrusion, which relies on a continuum theory of these processes, to the detailed analysis of the mechanism and energetics of liquids entering in cavities with different characteristics. Finally, we consider the effect of flexibility of porous/textured media on intrusion/extrusion, an aspect that has been investigated only recently. We show that one can use intrusion/extrusion in flexible materials to trigger the occurrence of unusual phenomena, like negative volumetric compressibility of unprecedented magnitude. We found convenient to consider the studies to which the authors of this article contributed, which are here often presented with a different perspective from the original papers and highlighting the contacts among the various subjects. We certainly recognize that many authors contributed to the development of this field and have inspired our research; the work of these authors is duly cited within the text to encourage the readers to consider additional results and different points of view.
The article is organized as follows: in Sect. 2, we discuss the general aspects of intrusion/extrusion and introduce a continuum theory to describe these processes. In Sect. 3, we present applications of the macroscopic theory to design superhydrophobic surfaces that can resist or recover their properties after extreme conditions are applied. In Sect. 4, we discuss the use of porous materials and non-wetting liquid intrusion/extrusion for damping vibrations or storing mechanical (or thermal) energy. In Sect. 5, we show how the intrusion/extrusion of liquids can be used to trigger negative compressibility. Finally, in Sect. 6, we draw some conclusions and present our outlook of the field. We also included three appendices in which we summarize the computational methods used in the works we analyze in this topical review, which make this article self-contained.

Macroscopic model of intrusion/extrusion
The objective of this section is to discuss the intrusion and extrusion processes according to a classical continuum model of a solid-liquid-gas system. This simple model is useful both because it allows to understand the basic physics of intrusion/extrusion and for a quantitative prediction of the energetics of the process in macroscopic and microscopic cavities down to few nanometers.
It is instructive to start from the simple case of the intrusion/extrusion from a 2D square cavity excavated from a flat surface (see Fig. 1a). The intrusion of a liquid into this cavity can be described resorting to a sharp interface model, i.e., a model in which the free energy of the system, or any other suitable thermodynamic potential, is the sum of bulk terms, the bulk free energy of the liquid, solid and gas present in the system, plus interface energies associated with the three pairs of phases, solid-liquid (sl), solid-gas (sg), and liquid-gas (lg). In the present derivation, we neglect additional terms, for example the line tension(s) [15], because it is not crucial to understand the general concepts of intrusion/extrusion we intend to discuss here. Nevertheless, we remark that other authors, including some of us, have investigated the effect of line tension on intrusion/extrusion [9,16] is convenient to develop the sharp interface model of intrusion/extrusion in the grand canonical ensemble, in which the thermodynamic potential is the grand potential Ω = Ω bulk +Ω interf ace . Within the grand canonical ensemble, the bulk grand potential of a single phase is the negative of the product of its pressure, P x , times its volume, V x , where x stands for any of the three phases: Ω bulk Interface terms are proportional to the area between pairs of phases, A xy , through the corresponding surface tension γ xy , where x and y denote two phases: Ω interf ace xy = γ xy A xy . Considering that the surface of the solid does not change during intrusion, the configuration of the system is defined solely by specifying the liquid-gas dividing surface, the interface Σ lg . Indeed, Σ lg and its intersection with the surface of the solid implicitly determine Σ sl , Σ sg , and the volumes of the three phases (Fig. 2). Thus, Ω is a function(al) of Σ lg : Ω(Σ lg ). Neglecting the contribution of the solid bulk term, which does not change during intrusion/extrusion, Ω(Σ lg ) reads: where ΔP = P l − P g and cos θ Y = (γ sg − γ sl )/γ lg is the definition of the (cosine of the) Young contact angle.
The last equality stems from the observation that the volume available to the liquid and vapor is constant, V = V l + V g as is the area of the exposed solid surface, A s = A sl + A sg . A solid is lyophobic when θ Y > 90 • , that is, when wetting a solid, i.e., increasing A ls , is associated with an energy penalty.
The model of the free energy of a three phase system, Eq. (1), will be used in the following to illustrate the fundamental mechanism, energetics, and kinetics of intrusion and extrusion. We will initially digress from this objective to illustrate the basic principles that allow one to identify the transition path, i.e., the most probable sequence of configurations of the system during the intrusion/extrusion process. We will start by considering the nucleation of an emerging phase in the bulk of a pre-existing one, e.g., the bulk nucleation of vapor within a liquid. In classical nucleation theory (CNT) [17], bulk nucleation is described in terms of the volume of the growing phase, V g in the paradigmatic example of nucleation of vapor in a bulk liquid. In the modern language of transition events [18,19], V g is the (putative) reaction coordinate of the process. The reaction coordinate of a process is a collective variable (CV), i.e., a function of the entire set of variables describing the system--in the case discussed above Σ lg -that is sufficient to characterize the sample along the transition-the nucleation and growth of a vapor bubble in the present case. For the time being, V g is just a putative reaction coordinate as it must be proved shows that Σ lg , together with its intersection with the solid surface, determines Vs, V l , Vg, Σ sl and Σsg, hence the configuration of the system that this variable is sufficient to characterize the process. Indeed, one immediately realizes that there is an infinite number of bubbles of volume V g = V * having different shapes Σ lg ; if homogeneous nucleation takes place following paths in which the bubble has a complex shape, V g is an insufficient CV at describing nucleation. We will discuss this problem in detail in Sect. 2.1.
Within CNT, originally developed for homogeneous nucleation, the growing bubble is assumed to be spherical, which can also be deduced from the following probabilistic considerations. Consider a bubble of volume V g ; this bubble can be bound by an infinite number of liquid-vapor interfaces, Σ lg (V g ). The probability to find a bubble in the configuration Hence, the configuration of maximum probability for a bubble of volume V g is the one with surface Σ * lg (V g ) minimizing the grand potential. Since all bulk terms are determined by the value of V g , minimizing the grand potential Ω(Σ lg (V g )) for bulk nucleation amounts to minimize the liquid-vapor interface, i.e., to find the interface Σ * lg of minimum area enclosing V g . In other words, an energetic problem is turned into a geometrical one. From basic geometry it is known that the 3D figure of minimal surface/volume ratio is the sphere, which is thus also the most probable shape of a vapor bubble of prescribed volume. It is worth remarking that this geometrical/energetic problem does not have additional solutions; so, for the case of vapor/gas bubble nucleation in the bulk, there are no local minima of the free energymetastable states-in which the system can be trapped. This has important consequences, considerably simplifying the bulk vapor nucleation path.
Let us now consider the probability to observe departures from a spherical bubble. The relative probability of observing two bubbles with the same volume but enclosed by surfaces Σ lg , say the spherical bubble, and Σ lg is where . Now consider, for example, a sinusoidal perturbation of Σ lg of suitable frequency, i.e., such that the frequency is consistent with the periodicity imposed by the radius of a sphere, and prescribed amplitude. Consider two spherical bubbles of different radii, both consistent with the frequency of the sinusoidal perturbation; clearly, ΔA is larger for the bigger bubble, as can be easily understood for the corresponding 2D problem considering that the length of a sinusoidal wave, L = 2πR 0 dx 1 + (A sin(νx)) 2 , grows with the radius of the circle. This observation implies that perturbations of prescribed amplitude (and frequency) are less probable in larger bubbles. Since any perturbation consistent with the fact that Σ lg is closed can be expanded in terms of sine/cosine waves, the general conclusion is that the departure from the spherical shape decreases with increasing size of vapor bubbles. This analysis is consistent with the empirical observation for the complementary problem of condensation that small droplets containing as few as 10-25 Lennard-Jones particles are already of spherical shape [20].
The above considerations allow one to conclude that the most probable vapor nucleation path at relatively low temperature and in quasi-static conditions, i.e., such that the probability to observe a bubble of given shape is p(Σ lg ; V g ) ∝ exp [−Ω(Σ lg (V g ))/k B T ], is a sequence of (almost) spherical bubbles of growing volume, as obtained from the minimization procedure discussed above. Importantly, there are no local minima of the free energy that can trap the bubble in other attractive basins, i.e., bubbles of underneath different geometries, and that the departure from this path decreases with growing bubble volume, as the energy penalty to produce perturbations also grows with size. A consequence of the above analysis is that V g is sufficient to characterize the nucleation path of vapor/gas bubbles, i.e., V g is the reaction coordinate of this process.
In confinement, we use a variational approach analogous to the one illustrated for CNT to identify the intrusion/extrusion path and the associated grand potential profile within the (continuum) sharp interface model: 1. The process, which in this quasi-static framework is the same for both intrusion and extrusion, is described in terms of the volume of the vapor/gas phase, V g , within the pores of the solid; V g decreases during intrusion and increases along extrusion. 2. For a given value of V g , one seeks the surface Σ * lg minimizing the grand potential Ω(Σ lg ) of Eq. (1). This procedure represents the extension of the classical nucleation theory to confined environments, which we denominated the confined Classical Nucleation Theory (cCNT) [21,22].
The application of cCNT to a simple 2D square cavity reveals that intrusion/extrusion is a process of unexpected complexity [21]. Hypothetical intrusion mechanisms were proposed in pioneering works, based on physical intuition: one in which the liquid enters in the cavity with an almost flat liquid-gas interface (Fig. 3a [23]) and another in which the liquid sags into the pore while being pinned at the top of corrugations (Fig. 3b [24]). On the other hand, cCNT shows that the process follows a more complex path: after an initial phase in which the liquid intrudes the cavity with a flat meniscus, a bubble is formed in one of the corners at the bottom of the pore (Fig. 1a [21]). This more complex behavior is related to the ruggedness of the grand potential landscape induced by the confining environment, which is characterized by the presence of several absolute/local minima of the grand potential at fixed V g ; importantly, the relative probability of these minima changes with V g . In other words, V g is not a good reaction coordinate in confinement: there exist several minimal bubble morphologies with the same V g and different Σ lg ; thus, one certainly needs additional collective variables to represent the state of the system along the process. We will come back on this problem in Sect. 2.1.
The complex mechanism identified by cCNT for the wetting of a 2D cavity was found to be consistent with recent confocal microscopy experiments [25], which showed the formation of asymmetric air bubbles in the corners of cylindrical pores during intrusion. Indeed, the experiments showed that the bubblein-the-corner was a possible mechanism by which a liquid wets a cavity; a second mechanism was observed in which the liquid enters into the cavity with a symmetric meniscus. This difference might be due to intrinsic differences between the experimental setup and theoretical hypotheses, first of all that insoluble gases (air) are not considered in cCNT. Moreover, forward flux sampling simulations [26,27] showed that the symmetric meniscus mechanism is more likely due to inertial effects: at the typical experimental pressures, when the intrusion barrier is small, and in absence of insoluble gasses, the liquid enters into cavities at a high speed and the meniscus does not have time to relax to the constrained equilibrium configuration. Of course, cCNT, which assumes local equilibrium, cannot model well enough highly out of equilibrium conditions. The observation that the intrusion mechanism is nontrivial might be relevant also in different contexts. For example, Patel and coworkers [28] identified nonclassical vapor nucleation paths for the formation of a bubble in a liquid between two hydrophobic plates. This nonclassical path, which results in a reduced nucleation barrier, corresponds to the abrupt transition from a spherical-cap bubble attached to one wall into a cylindrical bubble joining the two plates. Indeed, this mechanism can also be predicted by classical continuum models, provided that one uses a suitable theory of confined nucleation, such as cCNT described above.
In Fig. 1b, the grand potential along intrusion is shown at several values of the pressure. One notices that, at low pressures, the stable state is the extruded one, which is associated with a lower grand potential. By increasing the liquid pressure, i.e. for ΔP 0, the bulk term in Eq. (1) favors intrusion, i.e., it reduces the grand potential of configurations with large V l . Thus, for sufficiently large liquid pressures the intruded state becomes the stable one. However, a barrier ΔΩ † separates the two states over a broad range of pressures; thus, if the system was initially in the extruded state, intrusion can be prevented for kinetic reasons, i.e., the system cannot overcome the barrier on the relevant timescale, e.g., the experimental one (from microseconds to hours). However, by further increasing the liquid pressure the barrier is reduced, until it reaches a value comparable with the thermal energy available to the system, k B T . Since the intrusion time τ grows exponentially with the barrier, τ = τ 0 exp[ΔΩ † /k B T ], it is only when ΔΩ † ≈ k B T that τ becomes short enough that intrusion takes place on the experimentally accessible timescales (Fig. 1c). This is the statistical mechanics picture of intrusion. This representation is apparently in contrast with the classical picture, deriving from force balance arguments, in which intrusion takes place at a well-defined value of pressure, P int . However, given the exponential dependence of τ on ΔΩ † , and the dependence of the latter on the liquid pressure, τ is found to sharply decrease around a specific value of the pressure. Thus, although one cannot identify a single value of the pressure at which intrusion takes place, the pressure interval within which τ goes from very long times to experimentally relevant timescales is narrow enough that, in practice, one can identify an intrusion pressure also within the statistical mechanics picture of intrusion. This reconciles the two pictures of the intrusion process, the classical and the statistical mechanics ones.
For the square cavity, it is interesting to illustrate the origin of the intrusion barrier, ΔΩ † , by an intuitive argument. By intruding the liquid in the cavity one increases the solid-liquid interface area at the expenses of the solid-gas one. If the solid is lyophobic, θ Y > 90 • and cos θ Y < 0, intrusion brings to an increase of the grand potential. This is the only relevant surface contribution during the initial phase of intrusion before the meniscus forms a bubble in the corner, because in this first phase the liquid-gas interface area does not change. However, to complete intrusion, also the bottom wall of the cavity needs to be wet, which requires the cancellation of a portion of the liquid-gas and gassolid interfaces (see Figs. 1a and 3a). Thus, the energetic cost of further wetting of the solid is compensated by the gain of removing two interfaces. Depending on the properties of the system, namely the values of the three surface energies γ xy , the final part of the wetting may be associated with a reduction of the grand potential, giving rise to the formation of a barrier. This mechanism at the basis of the intrusion barrier was recognized by Patankar in his seminal work on the wetting of surface textures [23,24], in which he imagined two mechanisms depending on the depth corrugations: flat meniscus intrusion and sagging (Fig. 3).
The arguments developed above are valid also for extrusion. The physical process in this case is the nucleation of a vapor bubble favored by lyophobic confinement [29]. Thus, also for extrusion one can identify a characteristic pressure, P ext , at which the extrusion barrier is of the order of k B T . Also in this case, a cavity initially intruded by a liquid can be kinetically trapped in this (meta)stable state. Depending on the morphology and size of the cavity and on the characteristics of the liquid, the extrusion pressure may be negative, i.e., extrusion can be achieved only under tensile conditions [30]. However, liquid rupture or formation of bubbles in the bulk of the liquid might occur before extrusion, preventing the observation of the latter process [31]. This is why in the literature it is often reported that wetting of cavities is irreversible [32]. One should remark that the extrusion barrier decreases with the characteristic length of the cavity; thus for nanometric lyophobic cavities (contact angle θ Y > 90 • ) extrusion can take place at positive pressures [33][34][35][36][37][38].
The presence of intrusion and extrusion barriers is also responsible for the hysteresis observed in liquid porosimetry during the complete cycle. In liquid porosimetry, a porous material, usually in granular form, is dispersed in a non-wetting liquid and the mixture is placed in a cylinder-piston system. During the experiment, the piston is displaced and the pressure of the liquid is measured, thus defining a P-V curve for intrusion; a similar process can be repeated for extrusion. Here, for the sake of simplicity, we illustrate the experiment in which the pressure is controlled while the volume of the mixture is measured. Before and after intrusion, one measures a moderate reduction of the volume of the mixture due to the liquid compressibility. On the contrary, at or close to P int the volume changes abruptly due to intrusion, shrinking by a value corresponding to the total pore volume. Correspondingly, by reducing the pressure starting from the intruded state, the volume initially increases because of liquid expansion. Subsequently, when the pressure approaches P ext , the liquid extrudes and the volume increases suddenly by an amount corresponding to the porous volume. Thus, liquid porosimetry allows for the experimental determination of the intrusion and extrusion pressures; intrusion-extrusion hysteresis means that P ext < P int and, as we discussed above, these two critical pressures are determined by intrusion/extrusion barriers. Thus, a key question in porosimetry is determining the chemical and morphological characteristics of the pore which determine the intrusion/extrusion pressures; more generally, such knowledge allows to establish the structurefunction relation for non-wetting porous materials, which is crucial both in basic science and in technological applications. In the following, we will present selected continuum calculations and atomistic simulations that clarify the relation between physical properties and intrusion/extrusion conditions; we will further apply this knowledge to design novel materials for technological applications.

Descriptors of intrusion and extrusion
Before moving to applications, we discuss some additional fundamental aspects of intrusion/extrusion. In the section above we have explained that the most detailed descriptor of intrusion and extrusion within the sharp interface model is the liquid/gas interface Σ lg : for a rigid solid and an incompressible liquid Σ lg completely describes the state of the system. Nevertheless, in experiments one can only (indirectly) measure the volume of liquid in the cavities, and theories and simulations have been developed on the basis of this descriptor of intrusion/extrusion (see, e.g., Refs. [16,28,29,39]). Departure from macroscopic predictions (of Eq. (1)) of the atomistic free energy profile for the confined nucleation of vapor bubbles have been interpreted as nanoscale effects not captured by the sharp interface model, i.e., effect that manifest themselves only for nanoscopic bubbles. In particular, the mismatch between atomistic and sharp interface (Eq. (1)) results has been used to quantify the so-called line tension, a term that is proportional to the length of the solid-liquid-gas triple line [15,16,28], which gives a significant contribution only for nanoscopic bubbles.
The presence of hysteresis in liquid porosimetry experiments, in which one controls the volume of the cylinder-piston system, suggests that the volume of the liquid in the cavity is an insufficient descriptor of the process. Indeed, hysteresis is connected to multiple attractive basins in the space orthogonal to the variable controlled [18,40,41], here the volume of the computational sample. In this section, we describe continuum calculations and atomistic simulations that brought us to identify a suitable descriptor for intrusion/extrusion.
Let us first analyze the results shown in Fig. 4. One notices that the transition state corresponds to the morphological transition of the meniscus from a symmetric morphology to the bubble-in-a-corner one. Within cCNT, where the descriptor of wetting is the volume of the liquid in the cavity, this transition is discontinuous. Let us define Φ † l the volume fraction of the cavity occupied by the liquid at the transition state, Φ † l = V l /V ; for Φ l < Φ † l the symmetric meniscus configuration has the lowest grand potential, while for Φ l > Φ † l the bubblein-the-corner configuration is the most stable one. Of course, the path in which at Φ † l the system abruptly changes configuration from symmetric to bubble-in-thecorner along intrusion (or vice versa during extrusion) in unphysical as no continuous dynamics can bring to such a discontinuous passage from one morphology to another.
To overcome this problem and identify the real intrusion/extrusion path of a simple 2D cavity, we have implemented the string method for a discretized representation of the sharp interface and a corresponding atomistic model [45] (See App. C). The string method [42][43][44] is a technique capable of identifying the most likely transition path when the process is characterized by large free-energy barriers separating the initial and final states. In this case, one cannot use standard continuum/atomistic techniques as the timescale of the process largely exceeds that achievable by brute force approaches. An additional advantage of the string method is its good scalability with the number of (collective) variables, which allows to use multiple CVs. A similar approach has been used also by Panter and Kusumaatmaja [46], Ren [47], Pashos et al. [48] to investigate the intrusion of non-wetting liquids into surface textures in combination with a diffuse interface model of the fluid. We studied intrusion as a function of (the discretized) Σ lg obtaining the transition path (a) (b) Fig. 4 a Comparison between the continuum (top) and atomistic (bottom) density fields of a liquid intruding a 2D cavity as obtained from string calculations [42][43][44]. One notices the similarities between continuum and atomistic results, suggesting that up to the ∼ 10 nm length scale macroscopic modeling gives a reliable (qualitative) picture of the process. b Free energy profiles as obtained from continuum and atomistic simulations for (i) a single CV measuring the amount of liquid in the cavity and (ii) for the vectorial CV representing the density field of the liquid. The close correspondence between continuum and atomistic results confirm the reliability of simple continuum models down to the ∼ 10 nm length scale. Reprinted from J. Chem. Phys. 142, 104701, 2015, with the permission of AIP Publishing shown in Fig 4a [45]. It is worth remarking that the string method, which is a minimization algorithm in the transition path space, allows to find the most likely path in the attractive basin of the first guess. In other words, the string method does not guarantee that the one identified is the only possible intrusion/extrusion path. Indeed, our analysis of the the intrusion in the 2D cavity by cCNT [21] revealed that several intrusion paths are possible including the nucleation and growth of a liquid droplet in the middle of the cavity or at its walls that will eventually bring to the complete filling of the texture. Nevertheless, studies based on a continuum model of the three phase system [21,46] and the corresponding atomistic representation [26,29,30] confirm that the one shown in Fig 4a is the most likely intrusion path among all possible paths.
One notices that the morphologies of the meniscus obtained from the string calculation and cCNT coincide for a large portion of the transition path, essentially departing from each other only close to the tran-sition state, in correspondence of the nonphysical morphological transition of the cCNT. Similarly, the grand potential profiles obtained for Σ lg and V l are very close to each other except in correspondence of the transition state. Knowing more precisely the configuration of the meniscus at the transition states and the value of the corresponding barrier might provide not only better estimates of the intrusion/extrusion kinetics but might help to devise strategies to enhance/delay these processes.
Atomistic simulations of the intrusion/extrusion process performed using i) the total number of fluid particles and ii) the fluid density field in the cavity (Fig 4a) as proxies of V l and Σ lg in the macroscopic representation, respectively, reveal a substantial accord between the macroscopic and microscopic intrusion/extrusion paths for cavities down to ca. 10 nm (Fig 4). The agreement between atomistic free energy and macroscopic grand potential of intrusion/extrusion is also remarkable (Fig 4b). A similar behavior has been observed for intrusion/extrusion in more complex cavities, such as those corresponding to nanopillars decorating a surface [41,49] and hydrophobic nanopores [50]. The analysis of the continuum vs atomistic intrusion paths serves also to test the accuracy of macroscopic models to describe intrusion down to the nanoscale. Indeed, if one takes into account that, due to the interatomic potential, a portion of the cavity of length ∼ 1σ, where σ is the characteristic size of atoms, e.g., the Lennard-Jones characteristic length, is inaccessible by the liquid, there is a surprising agreement between the macroscopic and microscopic energetics. Here, one must remark that this agreement has been obtained without any fitting [45,49,51]; indeed, parameters of the macroscopic model, Eq. (1), have been determined by independently measuring the surface tension γ lg by the Kirkwood-Buff method [52] and the contact angle θ Y by depositing a droplet on a flat surface. V l is determined from the average density ρ within the cavity of volume V assuming that the system is two phase, with liquid and vapor having bulk densities ρ l and ρ g , respectively: Finally, A lg and A ls are determined by the marching cube method [53], a technique that allows to reconstruct interfaces from a scalar field, here the density field of the fluid.
In principle, for atomistic simulations one can perform the committor analysis to test whether the density field is a suitable descriptor of intrusion/extrusion. The committor analysis consists in determining what is the probability to reach the final (initial) state first starting from a given state. At the transition state the committor must be distributed around probability 0.5 [18]. However, the collective variable corresponding to the (discrete) density field has a high computational cost, because it requires to count of the number of particles within (tens to hundreds of) cells of edge ∼ 3Å(∼ 1 σ) and by applying some smoothing to prevent impulsive forces on particles located at the borders the cells [54,55]. The computational cost forces one to discretize the string in relatively few images (typically 32-64), to limit the number of MD simulations to be performed per string iteration. This limitation, combined with the rapid change of the free energy around the putative transition state (the free energy barrier at the thermodynamic conditions of the simulations we have performed is several hundreds of k B T [41]), prevented us to identify an individual image along the transition path with a committor distribution centered around 0.5, which would be the final evidence that the density field is a good CV. Also refining the string near the putative maximum of the free energy and performing the committor analysis to check whether the distribution of the committor improved is computationally very challenging as the duration of simulations bringing the liquid from the transition to the initial or final states is very long, one 1ns or more per trajectory, requiring, say, 10 × 1 ns ≈ 100 ns to estimate the committor for each initial configuration at the transition state, to be repeated for 10-100 initial configurations, for a total of 1-10µs of MD simulations. This is a much longer simulation time than that required for typical molecular problems (e.g., conformational isomerization), for which a careful committor analysis has been performed. Nevertheless, the distribution of the committor around 0 and 1 for the two images before and after the putative transition state reassured us that this CV is adequate [45].
It must be remarked that the discrete density field CV is not uniquely defined, as it depends on the size of the cells used to discretize the space. In Refs. [49] and [41], we proposed a qualitative and computationally inexpensive test based on the idea discussed above that for a CV to be suitable it must produce a continuous transition path. Thus, (i) one can compute the transition path using a density field with a prescribed discretization of the space and, after convergence of the string calculation, (ii) compute the density field along the string using a finer discretization of the space based on atomistic trajectories at the last string iteration; (iii) compute the Cartesian distance between the density field at successive images: d i = j ρ i j − ρ i−1 j , with i the index of images along the string and j the index of the cells used to discretize the space. If the transition is continuous and the string is parametrized with equally spaced images, one expects the distance between density fields to be (almost) constant (see Fig. 5). For nonuniform discretization of the string the distance must be properly re-weighted. Figure 5 also illustrates the difference in the intrusion path as predicted via a single CV, V g , or with the density field for two discretization levels of the density field [41]. For the intrusion path determined with these three CVs, we report Σ lg obtained from the density field on the fine grid. Σ lg is represented in color code: it is blue when it is at the top of the pillars decorating the surface and red when it is at the bottom of textures, passing through green and yellow at intermediate heights. One notices that with the single collective variables (first row of Fig. 5a), the liquid intrudes the surface textured in blocks, e.g., the space in between two entire rows of pillars, then the space between four pillars and so on. With the discretized density field CV, on the contrary, the process is continuous: the liquid intrudes all the textures, eventually touching the bottom of the textures between two pillars. The general aspects of the intrusion path are analogous for the coarser and finer discrete density fields, with some subtler differences still present during the contact with the lower wall, that appears continuous only with the finer density field CV.
We conclude this section by summarizing some results: -The liquid/gas interface, the meniscus shape Σ lg , is a suitable descriptor of the intrusion process. -A suitable microscopic (atomistic) analogue of this descriptor is the discretized density field. Typically, the adequacy of degree of discretization can hardly be tested by a proper committor analysis. Nevertheless, an heuristic test based on the Cartesian distance between density fields at successive images along a string has been proposed. -We could not identify a universal principle for choosing a suitable discretization of the density field. This must be sufficiently fine to capture the curvature of the meniscus during liquid intrusion, especially high curvatures that occur when the meniscus transforms from one morphology to another (see Fig 4a/b), which, in turn, depends on the size and morphology of the cavity. -In some cases, these descriptors are difficult to use in macroscopic calculations and in atomistic simulations. In particular, the cost of computing the density field in atomistic simulations adds a significant overhead to simulations that makes this approach unsuitable, especially when one wants to gain a qualitative understanding of the process. -Fortunately, the simpler V l descriptor is sufficient to characterize the process, both qualitatively and quantitatively, for a large part of the transition path, and to give reasonable estimates of the freeenergy barrier.

Wetting and recovery of superhydrophobic surfaces
Materials possess an intrinsic level of lyophobicity, the property to repel a liquid, that depends on their surface chemistry. Lyophobicity, which is named hydrophobicity when the liquid is water, is usually measured by the contact angle of the liquid, i.e., the angle formed by the tangent to a droplet deposited on a surface at the triple solid-liquid-gas line. Lyophobic surfaces have a contact angle greater than 90 • , lyophilic ones have a contact angle smaller than the same number. The larger the contact angle, the more lyophobic the surface. For an (atomically) flat surface, the contact angle corresponds to the Young one, which depends on the inter-face energies between the three phases (see Eq. (1)). As a matter of fact, no surface of simple, atomistically flat materials exist with a Young contact angle larger than ∼ 110 • with the most common liquid, water. However, it is common experience to observe larger values of contact angles for corrugated (non-atomistically flat) surfaces. This is because, owing to capillary forces, a droplet can remain suspended above the asperities present on the surface in the so-called fakir state, also know as the Cassie-Baxter state [56], after the name of the researchers who first described this state. In the Cassie-Baxter state, the solid-liquid contact area is reduced; this fact, in addition to an enhanced contact angle, confers to the liquid a series of properties that we will discuss shortly. Indeed, the contact angle of an intrinsically lyophobic surface is enhanced by surface corrugations even when the liquid wets their bottom, i.e., when the liquid is in the so-called Wenzel state [57], also named the impaled state. There are several other properties which make the Cassie-Baxter and Wenzel states very different between them, among others, contact angle hysteresis, i.e., the difference between the advancing and receding contact angles, θ a and θ r , the contact angles measured when the size of the droplet increases or decreases, respectively, or the front and rear contact angle for a sliding droplet. Empirically, the larger the contact angle hysteresis, the stronger the liquid adhesion to the surface. In the Cassie-Baxter state, in which the liquid-solid contact is reduced, the adhesion of a liquid to a solid is lower than on the corresponding smooth surface. On the contrary, in the Wenzel state, in which the actual liquid-solid contact area is larger than for a smooth surface, liquid adhesion is increased. Thus, for example, a liquid in the Cassie-Baxter state can easily roll-off of a surface, requiring a minimal force or a small tilting of the surface. This property, together with other technologically relevant features that we refrain from discussing in detail, goes under the name of superhydrophobicity. Superhydrophobicity is often observed in nature; for example, lotus leaves exploit this property to keep their surface clean: rain droplets carry the dirt with them as they roll off the surface of this plant [58]. Another natural example in which surperhydrophobicity is seen in action are little "bugs" skating on water: water striders move on the surface of water thanks to the action of capillary forces and the peculiar hairy surface of their legs, which prevent water from wetting them, leaving the liquid in the fakir state [59,60].
Superhydrophobic surfaces have numerous technological applications: they can be used, for example, to prevent metal corrosion [61], to enhance water repellency, to realize self-cleaning optical windows and lenses, to reduce viscous drag [62], to impart anti-icing and antibio-fouling properties [63], for water desalination and many more applications [61,[63][64][65][66]. What we want to stress here is that the ideal superhydrophobic surface must present (regular or irregular) surface corrugations which reduce the solid-liquid contact area. Moreover, these corrugations must be designed so that under operative conditions, with possible variations of pressure and temperature, they resist wetting, i.e., the liquid cannot intrude the surface cavities within the experimental timescale or, if intrusion takes place, the liquid is able to spontaneously extrude upon restoring normal conditions (self-recovery superhydrophobic surfaces [33]).
An intense experimental and theoretical research activity has been devoted to investigate the fundamental mechanism of wetting and recovery, and its dependence on the characteristics of surface textures. Li et al. reported that a self-recovery textured surface can be obtained by suitably optimizing the chemistry of the material and, in particular, adopting hierarchical surface textures [67], similar to Verho et al. who considered dual-scale superhydrophobic surfaces [68]. Checco et al. have shown that textured surfaces with partial self-recovery can be fabricated on a large-scale exploiting the self-assembling of block-copolymer followed by plasma etching [69]. Vrancken et al. associated self-recovery properties of selected superhydrophobic surfaces to the directionality of their patterning, which allowed the triple solid-liquid-gas to dewet the grooves with negligible or no energy barrier to be overcome [70]. Recovery can also be achieved by applying external stimuli, e.g., electric [71,72] or magnetic pulses [73], if the liquid is superparamagnetic in the latter case. Of course, recovery by external stimuli is more complex and expensive in applications than the passive means mentioned before. In parallel to these (and other) experimental studies, an intense theoretical research has tried to establish the fundamental principles of (self-)recovery. Among the others, Prakash et al. have shown that inserting additional hydrophobic features in the cavities of surface textures reduces the barrier the liquid has to overcome along the extrusion path [74]. This can be easily understood noticing that these additional features increase the contact area between water and the hydrophobic solid for a given volume of surface corrugations. Using diffuse interface methods, Kusumaatmaja et al. came to the conclusion that a liquid droplet can penetrate into pillars decorating a surface by two mechanisms: by depinning from the top of the pillars, overcoming a free energy barrier, or by intruding in the interpillar space as a result of the increase of the droplet curvature due to liquid evaporation [75]. The authors concluded that wetting can be prevented by decorating surfaces with tall pillars or, by extension, tall textures.
To show how the theory presented in the previous sections can be instrumental for designing superhydropho-bic surfaces, we report cCNT investigations of the path, energetics, and kinetics of intrusion and extrusion of water for typical surface textures. We aim at finding the geometry of corrugations with the highest possible barrier opposing intrusion, so that if the thermodynamic conditions favor the Wenzel state, the system can still remain in the original Cassie-Baxter state for kinetic reasons until the unfavorable (pressure/temperature) fluctuation is finished [33]. At the same time, or alternatively, we looked for textures minimizing the extrusion barrier to obtain values lower than that allowing recovery on a timescale short as compared to experimentally relevant timescales. The self-recovery strategy is that, if liquid intrusion in the corrugations accidentally takes place, the system can spontaneously and quickly return to the Cassie-Baxter state once the fluctuation that induced the transition has ceased. In practice, we try to provide a theoretical foundation to the experimental [76][77][78][79] and heuristic [80,81] design of superlyophobic (hydrophobic, oleophobic, omniphobic) surfaces, additionally requiring that they remain permanently dry [34,82]. In particular, in these previous studies it was found that a judicious use of morphology of the corrugations, namely the geometry of their top in contact with bulk liquid, may allow to fabricate corrugated surfaces resisting to the intrusion even of liquids with low surface tension. Within a specific morphology of surface textures, e.g., pillars with re-entrant square head, genetics algorithms have been employed to heuristically optimize their geometrical characteristics, e.g., spacing between pillars, size of the reentrant head [80]. We aim at rationalizing these findings within a unifying theory. With this aim in mind, we considered a set of surfaces decorated with regular corrugations among the ones more commonly used in experiments: surfaces decorated with pillars, ridges, cylindrical pores and square pores (see Fig. 6a). Assuming that the operative conditions of these surfaces are ambient conditions, that the liquid is water and that the intrinsic (Young) contact angle of the material used to fabricate the surfaces is θ Y = 125 • 1 , we computed the free energy profile of intrusion/extrusion using the cCNT presented in Sect. 2. Given the complexity of these geometries, this calculations is solved numerically using the soft-1 θY = 125 • used in the present analysis is only seemingly inconsistent with the maximum value of the contact Young declared above for water, ca. 110 • . We remark that this value is about a simple solid, whose surface is made by a compact layer of atoms/molecules. Here, however, we consider the more general case of either (i) hierarchical rough materials, in which the material used to fabricate surface texture contains a second level of cavities of nanometric or subnanometric textures, of (ii) crystalline porous materials, such as MOFs, and zeolites, presenting surface cavities revealed by the cutting or pores intrinsic to the crystallographic structure. In both cases,the contact between solid and liquid atoms/molecules is reduced, which can increase the contact angle. Indeed, preliminary, unpublished simulations revealed for ZIF8, the zinc zeolitic imidazolate framework, a contact angle of ∼ 130 • . ware SurfaceEvolver [83]. We are aware of the limitations of this approach, which we discussed in detail in Sect. 2.1. In particular, we know that this approach might lead to underestimating the free-energy barrier. Nevertheless, here we are more interested in the qualitative aspects of the free energy profile and its connection to the morphology and size of the surface textures than in the quantitative details. Moreover, by comparison with brute force molecular dynamics simulations of liquid extrusion from selected cavities, we show that the limitations implicit in the present approach still allow to predict the correct order of magnitude of the kinetics of the self-recovery process.
A key finding of this research is that the extrusion barrier dramatically depends on the morphology of surface textures. In particular, only the pores, especially the square ones, present a sufficiently low extrusion barrier to allow self-recovery of wet surfaces at ambient conditions (see Fig. 6a). One must remark that even for a seemingly simple pore of this type the extrusion path is rather complex, with several changes in the morphology of the liquid meniscus (see Fig. 7). Initially, extrusion takes place by forming a gas bubble at one of the bottom corners; the bubble subsequently expands until it reaches a critical size beyond which it develops along one of the edges at the bottom of the pore, joining two corners. Upon further growth, it develops first along two edges, joining three corners, hence along all edges and joining all corners, with the liquid touching the bottom of the cavity only in the central part. Then, the liquid detaches from the bottom and forms a regular meniscus in contact with the four lateral walls of the cavity with a contact angle θ Y = 125 • .
Another key finding of this research is that, at ambient conditions, the transition state occurs at the begin- ning of the extrusion process, in correspondence of the formation of a tiny bubble at a corner at the bottom of the pore, for further theoretical considerations on the recovery process see also Ref. [82]. This apparently minor detail has, indeed, important consequences on applications: the functional characteristics of superhydrophobic surfaces are generally connected to the morphology of their top, i.e., to the part of the corrugations which is in contact with the bulk liquid when this is in the Cassie-Baxter state. Thus, the domains relevant for the recovery and for the functional characteristics are placed at distinct locations, at the bottom and top of the cavity, respectively. This feature can be exploited to design the two domains separately: the bottom can be tailored to enhance extrusion and the top to minimize the solid-liquid contact. This led us to propose a modular design, implanting atop of a square pore pillars or ridges, which reduce the contact area between corruga-tions and bulk liquid (Fig. 6a). To prove this concept we computed the free energy profile of extrusion for these two modular morphologies and have found that, while the modular textures with square pore+ridges do not significantly alter the free energy of the process, in the case of pillars implanted on the squared pore a new intermediate minimum emerges. This new metastable state, which we denominated internal-Cassie-Baxter state, is separated from the proper, fully extruded, Cassie-Baxter state by a sizable barrier which prevents the self-recovery of the superhydrophobic state (complete extrusion). Thus, although modular surfaces allow to reconcile the seemingly incompatible characteristics of surface textures with small liquid-solid contact and self-recovery capability, care must be paid to the emergence of intermediate metastable states introduced by joining two corrugation morphologies, such as pores and pillars. The analysis reported above is based on several hypotheses, i.e., that a sharp interface model can adequately describe a liquid intruding in and extruding from a nanometric cavity (∼ 10 nm), that using the volume of liquid in the cavity as the CV does not lead to a severe underestimation of the actual barrier, the one associated with Σ lg (see Sect. 2.1), that the estimation of the pre-exponential factor τ 0 for the determination of the extrusion time from the Blander-Katz formula [14] is accurate enough, at least of the correct order of magnitude. To validate these hypotheses we have computed the free energy profile of extrusion of water from the squared pore+ridges morphology from atomistic simulations and compared it with the results stemming from the sharp interface model, obtaining a surprisingly good agreement between results obtained with the two approaches (Fig. 8). We also performed brute force molecular dynamics simulations starting from the fully intruded state and measured the time it takes to recover the fully extruded Cassie-Baxter state when the barriers are low, confirming that few nanoseconds are sufficient to accomplish this task and that modular surfaces are good candidates for a new generation of multifunctional textures solving the limitation of standard ones.
Concerning the Blander-Katz formula [14] for the pre-exponential factor, Menzl et al. [84] have validated its accuracy for the specific case of liquid/vapor transition; this is further confirmed by the consistency between brute force MD and the extrusion timescale obtained from present free energy calculations for the corresponding process under confinement.

Lyophobic porous systems for absorbing vibrations and storing energy
Another field of applications connected to intrusion/extrusion of liquids within cavities is the damping of vibrations and storage of energy [5,6,9,38,[85][86][87][88]. These applications may seem antithetical, but they are strongly related. Damping of mechanical vibration or absorption of shocks via liquid intrusion/extrusion is possible owing to the large hysteresis between the two processes, associated with the presence of free-energy barriers [89,90], indicating that the mechanical energy spent to trigger intrusion-the energy of the vibration or shock in this case-is dissipated in other forms. On the contrary, when hysteresis is negligible or small, the energy spent to intrude the liquid in the porous material is returned upon extrusion, making the heterogeneous system of porous solid and non-wetting liquid a sort of mechanical battery or molecular spring. This short description of the two processes makes it clear that to develop devices implementing these technologies one must be able to design materials with (i) controlled intrusion/extrusion hysteresis and (ii) such that intrusion and extrusion take place at operative conditions, e.g., in the range of pressures accessed during the operations. Since during intrusion the energy spent for triggering the process is (mostly) converted into the interface energy between the liquid and solid, it is also important that conditions (i) and (ii) are met for materials having very large (internal) surface area density, the maximum liquid/solid interface area per gram of the solid. Porous materials with several hundreds to thousands of square meters per gram are suitable for the applications mentioned above. Suitably coated porous silica [31,[91][92][93], hydrophobic metal or covalent organic frameworks (MOF and COF, respectively) typically satisfy these conditions, although the complexity of the structure of internal surfaces and the nanometric size of these cavities prevented, so far, to identify the characteristics of the pores determining the intrusion/extrusion behavior.
An intriguing case is that of two porous silica materials, WC8 and RPB. These systems contain cavities that can be coated to make the materials hydrophobic [31]. RPB and WC8 present cavities of similar sizes, which can be intruded by water through surface apertures with similar characteristics. Thus, one expects similar intrusion/extrusion cycles between RPB and WC8. However, while intrusion takes place at similar pressures, ∼ 15 MPa, WC8 extrudes at ∼ 2.5 MPa and RPB does not extrude at all (Fig. 9c). Thus, CH3-coated WC8 is suitable for damping vibrations thanks to the large intrusion/extrusion hysteresis, which, however, closes at positive pressures allowing repeated intrusion/extrusion cycles. On the other hand, RPB cannot be used for this purpose because, once intruded, water remains in the pores; hence, RPB acts as a bumper suitable only for the mitigation of a single impact. Thus, given the analogous nominal characteristics of these two materials, a question arise: what determines the differ-ent extrusion behavior between WC8 and RPB? Transmission Electron Microscopy (TEM) images revealed that RPB cavities are independent, cylinder-like, with a smooth internal surface (Fig. 9a). On the contrary, WC8 pores consists of interconnected cavities, resulting in nanoscale apertures on the internal surface of the main pores (Fig. 9b). Can these different characteristics be at the basis of the different extrusion behavior?
To address the above question we investigated the formation of a gas bubble within a spherical cap cavity with and without nanoscale apertures on the internal surface (Fig. 9d), mimicking the main characteristic differentiating RPB and WC8. A preliminary continuum analysis (Sect. 2) showed that extrusion from this geometry follows a simple path, corresponding to the formation and growth of a spherical cap bubble at the internal liquid/solid interface. We thus performed atomistic free energy calculations, using a single collective variable measuring the size of the bubble, via, e.g., the total number of water molecules in the cavity, which is sufficient to describe the process. At ambient conditions, the free energy profile of the smooth spherical cap cavity presents two minima, one in correspondence of the fully extruded state and the other in correspondence of the fully intruded one, denoted by Φ = 0 and Φ = 1, respectively. Values of Φ larger than 1 correspond to compression of the liquid within the cavity. At ambient conditions these two minima are separated by a barrier ΔΩ † of several hundreds of k B T (Fig. 9d), which results in an extrusion time (τ = τ 0 exp[ΔΩ † /kBT ] -with τ 0 determined by the Blander and Katz formula [14])longer than the age of the universe. Thus, if the system is intruded by applying a suitable pressure, ∼ 15 MPa in the two experimental cases, water within the porous material with smooth internal cavities (as RPB) cannot extrude, remaining kinetically trapped inside the pores.
Let us consider now the same cavity decorated with nanoscale apertures. These apertures do not get intruded by water, neither at ambient pressure nor at ∼ 15 MPa, the intrusion pressure, because of the extreme confinement and hydrophobicity [34]. The effect of these apertures on the free energy profile is (i) to destabilize the intruded state with respect to the extruded one (ΔΩ = Ω int − Ω ext is larger in presence of the nanoscale apertures) and (ii) to decrease the extrusion barrier ΔΩ † (Fig. 9d). The higher the internal surface area decorated with nanoscale apertures the larger ΔΩ and the smaller ΔΩ † . Detailed atomistic simulations showed that, for nanoscale apertures covering ∼ 50 % of the internal pore surface, the extrusion barrier vanishes at ambient conditions. This result shows that, in addition to pore size and surface chemistry, one can tune intrusion/extrusion characteristics of a porous material for either vibration damping or energy storage by engineering the internal surface of porous materials, in particular, the pore connectivity.
The effect of empty nanoscale apertures, which reduce the liquid/solid contact area, can be rationalized considering that their presence increases the effective hydrophobicity of the internal surface of the main pores according to the so-called Cassie-Baxter effect [56]. This is shown in Fig. 9e, presenting a trend of the extrusion barrier obtained from cCNT consistent with the atomistic predictions.

Intrusion/extrusion of liquids in flexible materials for exceptional negative compressibility
In the previous section, the discussion of applications of heterogeneous systems composed of a non-wetting liquid and a nanoporous material suggested that it is convenient to reduce the size and control the geometry of pores, to facilitate and engineer extrusion and increase the pore surface area. Among several crystalline porous materials, Metal Organic Frameworks (MOF) are an emerging class of materials for energy and other technological applications [94][95][96]. Intrusion/extrusion of liquids in MOFs, which are highly flexible materials, add some complexity to the experimental and theoretical investigation of this process. For example, Debenedetti and co-workers have shown that thermodynamic stability and evaporation kinetics of water in hydrophobic confinement is extremely sensitive to the flexibility of the confining material [97,98]. However, the flexibility of MOFs opens new perspectives for technological applications of liquid intrusion/extrusion. In particular, we have recently discovered that intrusion of liquids inside selected flexible MOFs brings to exceptionally large values of negative volumetric compressibility [99]. Here, we will focus on this latter subject, which gives us the opportunity to discuss intrusion/extrusion at an atomic scale for systems in which the macroscopic modeling is expected to fail.
Let us begin by illustrating experimental results obtained on a sample of ZIF8 dispersed in water. ZIF8 is the acronym of Zeolitic Imidazolate Framework, a hydrophobic MOF (apparent contact angle 147 • -see inset of Fig. 11a) containing large cavities of the shape of truncated octahedrons (Fig. 10a,b) whose macroscopic form is that of a granular material. When dispersed in water at ambient pressure the internal cavities of ZIF8 remain dry. By applying an increasing pressure on the water-ZIF8 mixture the overall volume of the sample initially gently reduces due to the compressibility of the liquid and solid (pre-intrusion branch of the red curve of Fig. 11a). At ∼ 25 MPa the volume suddenly shrinks; this event is associated with the intrusion of water in ZIF8 cavities. Once ZIF8 is fully intruded, any further compression results in a reduction of the overall volume with a volume/pressure slope similar to that observed in the pre-intrusion stage (postintrusion branch of the red curve of Fig. 11a). In this post-intrusion stage, the liquid and porous materials are compressed. One can then proceed along the reverse path, decreasing the pressure starting from the final stage of the intrusion process (blue curve in Fig. 11a). Here, the sample first expands gently, then extrudes, though the extrusion branch of the curve is different  In Ref. [99], an intrusion/extrusion experiment was performed while collecting neutron scattering data, from which it was possible to determine the trend of the lattice parameter of the cubic cell of ZIF8 vs pressure (Fig. 11b). One notices that in the pre-intrusion and post-intrusion stages of the process, and in the corresponding stages along extrusion, the lattice parameter changes according the usual positive compressibility of conventional materials: the lattice parameter reduces when the pressure increases and expands when pressure decreases. However, in correspondence of intrusion (extrusion) the materials expands (contracts) while the pressure increases (decreases) as a consequence of the liquid entering (leaving) ZIF8 cavities. During intrusion the coefficient of volumetric compressibility of the porous material, κ = −1/V (∂V /∂P ) T , is negative with an absolute value an order of magnitude larger than previously reported materials with the largest negative compressibility. Thus, this research inspired by the investigation of materials for energy storage and vibration damping led us to identify novel strategies for producing giant negative compressibility based on liquid intrusion/extrusion, which can be used, for example, to regulate the liquid flow in nanofluidic devices, i.e., design and fabricate nanoapparatuses playing the role of a nanovalve: grains of ZIF8 are introduced in a nanochannel, if the pressure rises above the intrusion pressure the liquid intrudes in the ZIF8 cavities, the grains expand and reduces or block the liquid flow (Fig.12).
Negative linear compressibility has been recently observed in another highly flexible MOF, Cu 2 L [100], whose structure significantly differs from ZIF8. Thus  [29,30,101] of the intrusion/extrusion process of water in ZIF8 (see also App. B). We recognize that the process is very complex and, as discussed above, a single CV counting the number of water molecules within the porous material is insufficient at determining intrusion/extrusion path and the corresponding free energy profile. Many questions need to be addressed before one can imagine possible CVs for this purpose. In fact, as shown in Fig. 10b, ZIF8 contains a series of cavi-ties connected by hexagon-like and square-like windows. Thus, a first question to address is whether intrusion occurs with the fluid penetrating into all cages with an almost uniform density that grows during intrusion (vapor condensation-like mechanism) or it fills the cavities one after the other (percolation-like mechanism). If intrusion occurs cage-by-cage, the next question is whether liquid percolation follows a specific path, with the liquid passing from a cavity to another of the 8 con- tiguous ones on the basis of some well defined pattern or if the percolation path is random. Moreover, as we have already observed for simple cavities [9,21,29,45], describing the intrusion or extrusion process of pores might require complex CVs, such as the density field of the fluid. Last but not least, some elements of the ZIF8 framework, namely imidazolate linkers, can take different orientations, which alters the internal morphology of ZIF8 cavities and of the windows connecting neighboring cavities, which represent further degrees of freedom to be taken into account in the CVs to study the intrusion/extrusion path in detail. This makes the detailed intrusion/extrusion mechanism of liquids in ZIF8, and possibly other MOFs and COFs, especially challenging as these processes likely involve also the degrees of freedom of the solid medium, that so far we did not consider in studying the intrusion/extrusion mechanism. However, the objective of this preliminary work is mainly exploratory and we used a single CV, the number of H 2 O molecules within a slab of ZIF8 embedded within two water slabs, N H2O (Fig. 10c).
The free energy of N H2O presents several minima, corresponding to the filling of individual ZIF8 cavities, separated by sizable barriers. One notices that the profile changes with pressure, with the absolute minimum passing from the state corresponding to the empty cages to the one of filled cages in 1 MP a, when going from 50 to 51 MPa (Fig. 11c). The difference in the simulated intrusion pressure as compared to the experimental value is probably due to limitations in the force field [102], which is not designed for ZIF8 in contact with water. Along the intrusion process, we measured the effective lattice parameter in our slab by computing the average distance between equivalent crystal sites of different unit cells forming the ZIF8 slab at the value of N H2O corresponding to the minimum of the potential of the mean force at the current pressure. In other words, we computed this average distance out of restrained MD trajectories with a value of N H2O corresponding to the minimum of the mean force at the current pressure. For the pre-and post-intrusion slab, the average distance between equivalent crystal sites is determined in standard MD simulations of the sample containing either the empty or filled slabs, respectively. One notices (Fig. 11d) that the trend of the lattice parameter with pressure is qualitatively consistent with experimental results. Moreover, the variation of the lattice parameter during intrusion is in quantitative agreement with the experimental value.
Having confirmed that the experimental phenomenology is reproduced by MD simulations, we analyzed our trajectory to investigate the origin of negative compressibility. We considered the variations of ZIF8 characteristics that can affect the lattice parameters. The most likely variations are those corresponding to soft degrees of freedom, namely the distance between Zn and nitrogen atoms of imidazolate, and rotations of the Zn(Im) 4 tetrahedrons (Im stands for imidazolate). Concerning this latter degree of freedom, at ambient conditions, before intrusion, Zn(Im) 4 tetrahedrons have an orientation such that the hexagon-like windows allowing percolation of water from cavity to cavity are starshaped (see Fig. 11e,i). The water intruding the cavities may induce a rotation of Zn(Im) 4 tetrahedrons, resulting in an expansion of the lattice. Indeed, with some surprise we observed that water resides in these windows, of the size of ∼ 6Å, once the main cavities are filled (see Fig. 10h). It is this ultraconfined water that induces Zn(Im) 4 rotation (Fig. 11e) and the stretching of Zn-N bonds (Fig. 10f) resulting in the expansion of the lattice upon extrusion.
The results described above allow us to confirm that negative compressibility is connected to the flexibility of the material, but the material must possess some structural characteristics allowing expansion, which in ZIF8 is given by the presence of star -shaped windows connecting the framework cavities. More research is needed to identify the typical characteristics of MOFs and other flexible porous material leading to negative linear and volumetric compressibility. Inspiration can be drawn from metamaterials with negative compressibility [103][104][105], which are also typically made of interconnected nodes that can expand or compress due to rotations of some of the units forming the network. In a positive feedback, some of the characteristics of flexible porous materials might led to propose metamaterials with novel structures.

Conclusions and outlook
In this review, we described some of the results of our research in the field of intrusion/extrusion of liquids in/from textured surfaces and porous materials. Here we mostly focused on the mechanical aspects of the process: intrusion/extrusion mechanism and its (free) ener-getics under the action of a mechanical perturbation, typically the increase or decrease of pressure. We have introduced a model, the confined classical nucleation theory, which constitutes a simple and yet comprehensive framework to understand liquid intrusion and extrusion processes. We have shown that this theory provides a (semi-)quantitative description of the process, providing free energetics in good agreement with atomistic simulations down to scales of tens of nanometers. On smaller length scales, or when the atomistic structure of porous materials is more complex, like in crystalline porous materials, the continuum description of the process is likely insufficient (work is in progress about this issue).
Many aspects of intrusion/extrusion have been intentionally left out of this review; some of them have not been discussed to keep the focus on the main aspects of liquid intrusion/extrusion, others because there are open challenges to be addressed in the future. For example, here we discussed intrusion only as a quasi-static process, so that the system, whether it is modeled at the atomistic or continuum level, can reach the (local) equilibrium at the current value of the collective coordinates. However, liquids have inertia and the comparison between simulations and experimental results revealed this can play a key role in the intrusion and extrusion of cavities, including those of nanoscale size [26]. However, these results require further investigation on a broad set of systems to ascertain the relevance of inertia in the intrusion and its dependence on the nature of the liquid and of the porous material. Concerning further applications of liquid intrusion/extrusion, we focused on the narrow field related to energy and materials with negative compressibility. Actually, intrusion/extrusion of liquids into/from porous materials is of great importance for a much broader range of applications, e.g., liquid-phase chromatography, liquid separation, artificial and biological nanopores [13], and many more. Here, we refrain from discussing too many applications and aspects at the same time with the objective to give the theoretical and simulation fundamentals of intrusion/extrusion phenomena for researchers interested in the field.
Let us conclude this article listing some challenges to be addressed in the close future, some of which concern well-established experimental evidences that are still waiting for a theoretical explanation and other that arise from very recent results. Typically, intrusion and extrusion are accompanied by a heat flux, the process can be either endo-or exothermic. The heat released/adsorbed in the processes is very important as, for example, one can use liquid intrusion/extrusion to collect thermal energy from the environment or to fabricate actuators reacting to thermal stimuli [106]. Moreover, it has been recently shown that some electric current is generated during liquid intrusion/extrusion in hydrophobic porous material [6,88], which allows to exploit intrusion/extrusion to produce electric current, a process at the basis of a project recently funded by the European Commission, ELECTRO-INTRUSION. The triboelectric process at the basis of generation of electric current is still a rather elusive phenomenon and even more obscure is the correlation between the efficiency of this form of energy conversion and the characteristics of the liquid and porous solid. Changing the pure liquids considered in this review to a solution adds another important dimension to intrusion/extrusion problem: for instance, interesting new phenomena and applications may arise in the case of salt solutions [107][108][109][110][111]. Dissolved gases also play an important role in confined extrusion [3,11,12], accelerating the process with direct relevance for the resolution of liquid-phase chromatography [4] and, potentially, to explain the general anaesthetic effect [10].
Funding Open access funding provided by Universitá degli Studi di Ferrara within the CRUI-CARE Agreement. This project has received funding from the European Union's Horizon 2020 research and innovation programme under grant agreement No. 101017858.

Data Availability Statement
This manuscript has no associated data or the data will not be deposited [Authors' comment: Data sharing not applicable to this article as no datasets were generated or analysed during the current study.] Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/.

Appendix A: Landau free energy
In this topical review, we extensively referred to the probability density of a state of a system, the corresponding free energy and techniques to compute them. Since the focus of this article is not the detailed discussion of these topics, we refrained from presenting them in the main text. However, with the aim of offering to the readers a self-contained reference text, especially for those who are new to these themes, we added a couple of appendices on the atomistic techniques mentioned in the work. The reader is referred to other review works [112][113][114] and textbooks [115,116] for a thorough presentation and discussion of these techniques and their theoretical foundations.
In statistical mechanics, any thermodynamic potential is related to the logarithm of a suitable probability density function (PDF) in the relevant ensemble. In this topical review article, we extensively considered the PDF of having a prescribed number of water molecules within the cavities 163 Page 20 of 24 Eur. Phys. J. B (2021) 94 :163 exist (see Fig. 13). In this case, one first determines the most likely path connecting the initial and final state starting from a first guess and then compute the free energy along this path. This approach returns both i) the mechanism of the process and ii) its energetics, free energy difference between the initial and final state and the barrier. Several techniques have been proposed to achieve this objective; here, we shortly describe the string method in collective variable [44], that we used in our research on liquid intrusion/extrusion in porous systems. We remark that in this appendix the focus is on the algorithm, readers are invited to read the original paper for the theory beneath this method.
The objective of the string method in collective variables is to identify the minimum free energy path, MFEP, the continuous curve in the space of collective variables satisfying the equation Qij(z(α)) ∂F (z(α)) ∂zj , where zi is the value of the ith collective variable, or of the ith component of a vectorial CV, α is a parameterization of the MFEP (usually, α ∈ [0, 1]), means "parallel to", F (z) = −kBT log(M (z)) is the free energy of more than one CV, with In other words, the free energy depends on the multidimensional PDF M (z), which the probability density that the ith CV, ξi(r) takes a prescribed value: ξi(r) = zi. In Eq. 10, Qij(z) is a metric matrix due to projection of the phase space onto the collective variable space, and is defined as [44] Qij(z) = ∇zξi · ∇rξj z (11) ≡ d z ∇zξi · ∇zξj e −βV (r) δ(ξ k (r) − z k ) d z e −βV (r) δ(ξ k (r) − z k ) ,.
The hypothetical case Qij = δij helps understanding the intuitive meaning of Eq. (9): MFEP if the path connecting the initial and final states, two minima of the free energy, along which the mean force has only parallel components. Within the hypothesis Qij = δij, the MFEP is like the path a hiker would take to go from one valley to the next one passing through a mountain pass: always against or along the gradient of the force. Of course, when Q is non-trivial, i.e., when Qij = δij, the path deviates from this example taken from normal life.
The string method is an algorithm that allows one to identify the MFEP. The string itself is a discretization of the path connecting two (meta)stable states and is usually parameterized according to its relative arc-length, α = zα z b |dz|/ ze z b |dz|, with z b and ze beginning and end of the string. The discrete points along the string are called images and are labeled with their position on the string α l , where l is the index of the images. Starting from a first guess, the simplified and improved version of the string method [118] used in our works consists of three steps: 1. Calculation of the free-energy gradient and of the metric matrix, at the current position of the images; Qij(z(α l , t)) ∂F (z(α l , t)) ∂zj ; (12) 3. Parameterization of the string to enforce equal arclength parameterization among contiguous images.
For sufficiently large t, the algorithm guarantees that z(α l , t) converges to the MFEP of Eq. (9) [44,118]. We remark that t is not a realistic time, it is only a pseudo-time introduced to define the pseudo-dynamics of z. Summarizing, the string method is a first order dynamics with the constraint of constant value of the parametrization α l following the generalized force −Q(z(α l , t)) ∇F (z(α l , t)). The reparametrization step, enforcing the constant value of α l , removes the effect of the parallel component of the mean force on images, hence avoiding the need to remove this component before evolving string images by a numerically inaccurate projection operation. Overall, the string method is a constrained minimization of the system at a set of states at fixed α l in a space equipped with the metricQ local to the basin in the path space to which the first guess path belongs to.